Atmospheric Chemistry and Physics Discussions

The evolution of two-dimensional drop distribu- tions is simulated in this study using a Monte Carlo method. The stochastic algorithm of Gillespie (1976) for chemical re- actions in the formulation proposed by Laurenzi et al. (2002) was used to simulate the kinetic behavior of the drop popula- tion. Within this framework, species are defined as droplets of specific size and aerosol composition. The performance of the algorithm was checked by a comparison with the analyti- cal solutions found by Lushnikov (1975) and Golovin (1963) and with finite difference solutions of the two-component ki- netic collection equation obtained for the Golovin (sum) and hydrodynamic kernels. Very good agreement was observed between the Monte Carlo simulations and the analytical and numerical solutions. A simulation for realistic initial condi- tions is presented for the hydrodynamic kernel. As expected, the aerosol mass is shifted from small to large particles due to collection process. This algorithm could be extended to incorporate various properties of clouds such several crystals habits, different types of soluble CCN, particle charging and drop breakup.


Introduction
The understanding of aerosol-cloud interactions contains large uncertainties that must be reduced to accurately estimate the impact of aerosols on weather and climate.One of the most problematic aspects of aerosol-cloud interactions is the collision-coalescence process that is a mechanism that modifies the aerosol distribution, i.e. the aerosol particles that are the nuclei for individual droplets are combined during the Correspondence to: L. Alfonso (lesterson@yahoo.com)coalescence process in the same way as the mass of the individual water droplets are merged.After the evaporation of the drop formed by coalescence, the residual aerosol particle has the mass of the original two nuclei.
The aerosol distribution becomes important as the cloud drops evaporate and the solutes are recycled into aerosols that can serve as CCN: the larger the mass of a hygroscopic aerosol, the lower the supersaturation needed to form a cloud droplet.In the marine environment, the aerosol recycling process is believed to be the major mechanism responsible for the bimodal shape of the aerosol size distributions (Flossmann, 1994;Feingold et al., 1996).The heterogeneous chemical reactions, which add nonvolatile solute to each cloud droplet, strongly depend on the salt content and pH of the droplet (Alfonso and Raga, 2004).Since aerosols also have a significant influence on cloud microphysics and cloud radiative properties, it is necessary to simulate aerosol processes realistically and with adequate accuracy.
In general cloud models with detailed microphysics describe the aerosol and cloud droplets with two separate onedimensional size distributions.With this approach only the average aerosol mass contained in cloud droplets of a particular size is predicted by the model and it is not possible to keep track of the spectral aerosol mass distribution within the cloud droplets.For the deterministic case, the aerosol processing due to collision-coalescence was addressed by Liu (1998) and Bott (2000) by extending the flux method to two-dimensional distributions.Within this framework each particle is characterized both by the mass of its dry aerosol nucleus and by its water mass.Nevertheless, an extension of the exact stochastic framework developed by Gillespie (1976) for a two parameter droplet spectrum has never been reported in the cloud physics literature.
Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Alfonso et al.: Monte Carlo simulations of two-component drop growth The main advantage of the stochastic approach, described in this paper, over deterministic methods is that it can be easily extended to include not only the solute mass, but other particle properties such as crystal habit, different populations of CCN, chemical composition and the breakup of droplets (Alfonso et al., 2006).
Here we apply the general multi-component algorithm described by Laurenzi et al. (2002) to the solution of the kinetic collection equation (KCE) in cloud models dealing with twodimensional microphysics.
The discrete, two-component KCE, which is an extension of the discrete one-component kinetic collection equation, is given as: Where N (m, n; t) is the average number of particles consisting of m and n monomers of the first and second kind respectively (with water mass from size bin m and aerosol mass from size bin n).The water mass in size bin m equals the volume of a droplet in the smallest (monomer droplet) bin multiplied by m, the aerosol mass in size bin n equals the volume of an aerosol in the smallest bin (monomer aerosol) multiplied by n.Species are the type of particles with a given m, n composition.For example N(2, 1; 0)=1000 means that as an initial condition we have 1000 particles, all of them with a water mass that is twice the mass of the monomer droplet (in our simulations the monomer droplet is 4.189×10 −9 g) and with aerosol mass that equals the mass of the monomer aerosol (7.414×10 −15 g).The continuous version of this equation is more familiar: In Eqs.
(1) and ( 2) K(m, n; m , n ) is the collection kernel, now dependent on the composition of coagulating particles.The discrete KCE (1) gives the time rate of change of the average number of species with water mass from bin m and aerosols from bin n as the difference of two terms, the first term gives the gain in the number of particles whose water mass is in size bin m, and aerosol mass is in size bin n.It is calculated as a sum of binary collections between drops: one with water mass from size bin m , and aerosol mass from size bin n , and another with drop mass from bin m-m and aerosol mass from bin n-n and the second term describes the average rate of depletion of (m, n) particles due to their coalescence with particles from other species.To solve Eqs. ( 1) and ( 2) initial conditions are needed: For the discrete case, we also put N (0, 0; t)=0 for every t.
The numerical solution of the KCE (1) and ( 2) is difficult due to the double integral and nonlinear behavior of the equation and several numerical techniques can be found in the literature.In cloud physics modeling, Eq. ( 2) was numerically integrated by the flux method developed by Bott (2000) and independently by Liu (1998), both assuming that the probability for the collision of two cloud droplets depends only on the water mass of each one and not on the mass of the aerosol nuclei.
Other methods are computationally more expensive, such as the previously mentioned Monte Carlo (MC) algorithm developed by Laurenzi et al. (2002).This method has the advantage that it can be employed to determine both the expectations and fluctuations for multi-component aggregation.On the other hand, the KCE may not be valid at longer time periods, when a single drop acquires a mass much larger than the rest of the population and becomes separated from the continuous mass spectrum.In such a situation, the statistical fluctuations at the high-mass end of the spectrum must be taken into account.The Monte Carlo method is also very useful while investigating the role of coalescence in redistributing the aerosol mass in early warm rain stages when the artificial, numerical broadening of the drop distribution must be avoided.

The stochastic algorithm
A detailed description of the stochastic algorithm for multicomponent aggregation of particles can be found in Gillespie (1976) and Laurenzi et al. (2002), and we briefly summarize it here.Consider a well-mixed and spatially homogeneous volume V in which particles belonging to N s distinct species are present.Each species is characterized both by its water mass and by the mass of its dry aerosol nucleus, ūµ = (u m , u n ), such that, a droplet with composition ūµ is a member of the µth species.After time t=0 the species will randomly coalesce according to: where A m,n and B m ,n are droplets with compositions ūµ = (u m , u n ) and ūν = (u m , u n ), respectively.The transition probabilities for coalescence events follow Laurenzi et al. (2002) and are given by: a(i, j )=V −1 K(i, j )n i n j dt≡ (5) Pr {Probability that two particles of species i and j (for i = j ) with populations (number of particles) n i and n j will collide within the imminent time interval Pr {Probability that two particles of the same species i (with population number of particles) n i collide within the imminent time interval} In Eqs. ( 5) and ( 6), K(i, j ) is the collection kernel, and V is the cloud volume.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.The transition probabilities Eqs. ( 5) and ( 6) are then represented by one index (a ν ).This stochastic model is solved using the algorithm introduced by Gillespie (1976) for chemical kinetics and modified by Laurenzi et al. (2002).The expected behavior of the system can be evaluated by averaging over many realizations of the stochastic process, described by the following steps: 1.At t=0, the event counter is set to zero and the initial number of species n 1 , n 2 , . .., n N is defined 2. The quantity α is calculated as: A random number r 1 is generated from a uniform distribution in the interval [0,1] and considering that 1−r 1 =r * 1 is also a uniformly distributed random number is calculated 3. A random number r 2 is generated from a uniform distribution in the interval [0,1].and a collision ("chemical reaction") chosen with index µ from the inequality The number of species is changed to reflect the execution of collision.
3 Model results

Comparison of the Monte Carlo algorithm with analytical solutions of the two-component KCE
In order to check the performance of the Monte Carlo algorithm, a simulation with a constant kernel was performed and compared with the analytical solution found by Lushnikov (1975).Although the same analysis was performed by Laurenzi et al. (2002), there are some differences with the results outlined here since in Laurenzi et al. (2002) the simulations were performed for different values of the constant kernel and with a large number of particles (10 000 and 20 000 particles in the initial species), and only a single stochastic experiment was run (infinite system approximation).In the simulations presented in this section, the initial number of particles was much smaller (60) and the average was calculated over 1000 realizations.Solutions to Eqs. ( 1) and ( 2) can be obtained for an important class of collection kernels, such as when the kernel depends only on the total number of monomers (droplets and aerosols) in each colliding particle.In this case: Lushnikov constructed an explicit form for the composition distribution for this type of kernel, which corresponds to coagulation of initially monomeric particles.In this case N (1, 0; 0)=c 1 and N (0, 1; 0)=c 2 , corresponding to the situation with initially c 1 droplets and c 2 aerosols.The composition distribution may be expressed as (Lushnikov, 1975): Where m + n n are the binomial coefficients, and N (m + n, t) is the number of particles composed of (m + n) monomers (m monomer droplets and n monomer aerosols).Lushnikov (1975) showed that N (m + n, t), for the type of kernels ( 10) is a solution of the one-component kinetic collection equation: In this case, N (i, t)= m+n=i N (m, n; t).The initial condition for Eq. ( 12) is N (i, t)=N 0 δ i,1 .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 K(i, j )=A, B(x i +x j ) and C(x i x j ) where x i  and x j are the masses of the droplets from bins i and j .For the constant kernel K(x i , x j )=A and a monodisperse initial distribution with concentration c 0 , the analytical size distribution of the discrete KCE has the form: The analytical solution of Eq. ( 1), calculated according to the expression Eq. ( 11) for the constant kernel K(x i , x j )=A, is compared with true stochastic averages over N r realizations of the stochastic process (in our simulations N r =1000) : where N r (m, n; t) is the number of particles for species with droplet mass from bin number m and dry aerosol mass from bin n in the r-realization of the stochastic algorithm at time t.
The Monte Carlo simulation was conducted for initially monomeric particles (droplets and aerosols) with concentrations c 1 =30 cm −3 and c 2 =30 cm −3 (N (1, 0; 0)=30 cm −3 and N (0, 1; 0)=30 cm −3 ).Long (1974) calculated the coefficients for the polynomials K(x, y)=A, B(x+y) and C(xy) approximating the one dimensional collection kernel when the largest of the colliding drops is smaller than 50 µm.For the constant kernel, he found a value of A=1.20×10 −4 cm 3 s −1 .We used the same value for the constant discrete two-component collection kernel: In our simulations, the monomer droplet is 10 µm in radius (droplet mass 4.189×10 −9 g) and the monomer aerosol is an ammonium sulfate aerosol, 0.1 µm in radius (aerosol mass 7.414×10 −15 g).The aerosol-water mass grid was chosen according to Droplet mass(i)=i×m 0 (i=1, .., N droplets ) and Aerosol mass(j ) = j ×n 0 (j =1, .., N aerosols ).Here m 0 and n 0 are the masses of the monomer droplet (4.189×10 −9 g) and the monomer aerosol (7.414×10 −15 g), respectively.In all the simulations, the cloud volume was set equal to 1 cm 3 .We have defined 30 bins for the water mass grid and 30 bins for the aerosol grid.The pure monomeric species are also considered (those containing pure droplets and pure aerosols).Then, the total number of species in our numerical experimental can be calculated as: Where N droplets and N aerosols are the number of bins for the water mass grid and the aerosol grid respectively.The last two terms in Eq. ( 16) account for the monomeric species (droplets and aerosols).In our case the maximum number of species that can be generated during the simulations is 960.
The solutions obtained from the Monte Carlo calculations (averaged over 1000 realizations) for the species N(1, 1; t), N(1, 0; t) and N (0, 1; t) are shown in Fig. 1.The analytical solution is also shown in Fig. 1 (represented by the solid curve), and indicates the good agreement between these solutions of the KCE (1).
The two dimensional discrete size distributions for a) the analytical solution given by Eq. ( 11) and b) the average over 1000 realizations after 100 s are displayed in Figs. 2 and 3.The same comparison is shown after 200 s (Figs. 4 and 5) and 400 s (Figs. 6 and 7), respectively.Note that the differences between the Monte Carlo averages and the analytical solution of the KCE are again negligible.
The one-dimensional distribution, which is a solution of the one-dimensional kinetic collection Eq. ( 12), can be obtained from the two-dimensional spectrum by integrating over the aerosol grid for any point in time, as: In Eq. ( 17) N aerosols , and N droplets are the number of bins (grid points) in the aerosol and water grid, respectively.Two other simulations were performed with different initial conditions: N (1, 1; 0)=100 cm −3 and N(1, 2; 0)=150 cm −3 , corresponding initially to 100 and 150 particles per cubic centimeter from species (1, 1) and (1, 2), respectively.From Eq. ( 17), a monodisperse initial condition for the onecomponent KCE can be obtained from the two-component initial condition as: For this particular case (constant kernel and monodisperse initial conditions) we can use the analytical solution Eq. ( 13) of the KCE in order to compare with the two-component Monte Carlo.
The drop size distributions calculated from the Monte Carlo (after t=50 and 200 s), which are obtained by integrating the particle distribution over the aerosol grid according to Eq. ( 17), and the analytical solutions of the KCE with constant kernel (A=1.20×10−4 (cm 3 s −1 )) from a monodisperse initial condition N (1; 0)=250 cm −3 are displayed in Figs. 8  and 9. Again, a good agreement between the two approaches is found.
As was remarked in detail by Laurenzi et al. (2002), the species accounting formalism outlined in Sect. 2 reduces both computer storage and simulation time.For example, for a simulation experiment with 300 particles per cubic centimeter at initial time (which is considerably larger than the 60 particles used in our comparison study) only 82 species    were created during a 3000 s. simulation.This process is handled by dynamic allocation of memory permitting calculations with thousands of droplets in the initial distribution.a comparison with numerical solutions of the discrete twocomponent KCE.The hydrodynamic kernel (Pruppacher and Klett, 1997) has the form: In Eq. ( 19), V (m) and V (m ) are the terminal velocities of the falling particles with masses m and m , respectively, which are calculated following Beard (1976).Values of the collision efficiencies E(r, r ) were taken from Hall (1980).
The numerical integration of Eq. ( 1) was performed using the fourth order Adams-Moulton predictor-corrector scheme (Gerald and Wheatley, 2004), following the approach adopted by Valioulis and List (1984) for the KCE (12).The use of the predictor-corrector method severely restricts the number of particle sizes, because of the computational cost, but this is not a limitation for the purposes of comparison between the two methods.For both the Monte Carlo and the finite difference scheme, droplet and aerosol masses are expressed as multiples of the mass of the initial particles (monomer droplet and monomer aerosol).
At present, only the extension of the flux method (Bott, 2000) to solve the two-component KCE is actually available, but in Bott's approach, the water and the aerosol mass doubles after 4 grid cells.Then it is difficult to compare the results with the Monte Carlo algorithm which has more refined droplet and aerosol grid.Therefore, the predictor-corrector method provides the resolution in droplet and aerosol sizes to detect small differences between the two algorithms.Constant Kernel (1.2x10 -4 cm 3 s -1 ) t=200 sec Analytical Solution From 1D KCE Monte Carlo  The performance of the predictor-corrector finite difference method for the two-component KCE (2) was checked by integrating the two-component spectrum for the Golovin (sum) kernel on the assumption that the probability for the collision only depends on the water mass of each droplet and not on the mass of the aerosol nuclei: Where b=8.83×10 2 cm 3 g −1 s −1 following Long (1974), and x m , x m1 are the water masses from bin numbers m and m 1 , respectively.As an initial condition: N (1, 2; 0)=50 cm −3 and N (1, 8; 0)=50 cm −3 (which corresponds to 50 droplets with water mass 4.189×10 −9 g and aerosol mass 1.483×10 −14 g, and 50 droplets with water and aerosol masses equal to 4.189×10 −9 g and 5.931×10 −14 g, respectively).The twodimensional distribution was integrated at each time over the aerosol grid (see Eq. 17) and compared with the analytical solution for the one-component KCE (see Eq. 12) presented by Golovin (1963) for the sum kernel and monodispersed initial condition N (1; 0)=N (1, 2; 0) + N (1, 8; 0)=N 0 =100 cm −3 .In Eq. (21a,b), is the gamma function, N 0 the initial droplet concentration, b the constant for the Golovin kernel, v 0 is the mass of the monomer droplet (4.189×10 −9 g) and t is the time in seconds.The results of this comparison are displayed in Fig. 10 where the size  distribution is specified by the fraction N(m, n; t)/N 0 , where N 0 is the total number of particles in the initial distribution.An additional test of the finite difference method was performed by utilizing for the solution of both the twocomponent and one-component kinetic collection equations the hydrodynamic kernel Eq. ( 19).Again, the integration was made on the assumption that the probability for the collision only depends on the water mass of each droplet.In this case, the two-dimensional initial distribution for the two-component KCE was bidispersed in the water mass grid with N (1, 1; 0)=50 cm −3 and N(2, 2; 0)=50 cm −3 , which corresponds to the initial condition for the one-component caseN (1; 0)=50 cm −3 and N(2; 0)=50 cm −3 .Figure 11 depicts the distributions after 600 and 900 s, respectively, for this case.
Finally, the numerical solution of Eq. ( 1) obtained with the finite difference method was compared with the true stochastic averages over 1000 realizations of the Monte Carlo algorithm (see Eq. 14) for the hydrodynamic kernel Eq. ( 19).The simulations were run with realistic initial particle distributions.Figure 12 shows the initial two-component spectrum for the simulations.The spectrum has a droplet concentration of 158 cm −3 .This distribution was obtained (following Liu, 1998) by assuming a gamma distribution function for the drop coordinate and an exponential distribution for the aerosol size coordinate.
The two-dimensional discrete size distributions for a) the numerical solution given by Eq. ( 1) and b) the average over 1000 realizations after 300, 900 and 1200 s are displayed in   Figs. 13, 14 and 15, respectively.As the purpose of this simulation is to detect possible differences between the two numerical approaches, the number of bins displayed in the figures was restricted.As can be observed, the differences between the Monte Carlo averages and the numerical solution of the two-component KCE are negligible.The physical analog to this integration is to completely evaporate all drops until the dry aerosol spectrum is left.As can be observed, there is a net loss for small particles and net gain for large particles.Due to the collection process, the aerosol mass is shifted from small to large particles, with particles larger than 0.8 µm after 1500 s.

Discussion and conclusions
The multi-component MC algorithm proposed by Laurenzi et al. (2002) and based upon Gillespie's (1976) stochastic approach to chemical reactions was implemented to simulate two-component droplet growth by stochastic coalescence.Within this framework, all assumptions included in the stochastic collection equation are avoided.Additionally it permits calculation of statistical fluctuations for twocomponent droplet aggregation.On the other hand, the continuous KCE may not be valid when a single drop acquires a mass much larger than the rest of the system and becomes separated from the smooth mass spectrum (Alfonso et al.,  ) as a function of aerosol radius for the hydrodynamic kernel for t=150 s (dashed line) and t=1500 s (solid line with diamonds) obtained after averaging over 1000 realizations of the Monte Carlo algorithm.

2008
).In such a situation, the statistical fluctuations at the high-mass end of the spectrum must be taken into account.
For the two-dimensional case each species is characterized both by its water mass and by the mass of its dry aerosol nucleus.Very good agreement was observed between analytical solutions and numerical solutions of the KCE and MC simulations.
Moreover, the above described algorithm could be easily extended to the multi-component case in order to include various other properties of clouds.In a more general implementation of this approach species can be defined as types of particles with several attributes (droplet radius, CCN composition, chemical composition, electric charge, etc.) as well as the breakup of droplets (Alfonso et al., 2006).For this case, the state of a k component system is defined by a set of drops with properties or compositions ūi = (u 1,i , u 2,i , u 3,i , ..., u k,i ) where u k,i denotes the amount of the component or the property k in species i.For example, for the ice phase, it may represent the crystal habit, or the ice crystal mass.Then, the transition probability Eq. ( 5) may be defined as the probability that a specific pair of particles (drops, ice crystals, aerosols) with set of properties ūi = (u 1,i , u 2,i , u 3,i , ..., u k,i ) and ūj = (u 1,j , u 2,j , u 3,j , ..., u k,j ) will aggregate in the next time interval.
The stochastic approach should make more feasible the modeling of highly complicated microphysical processes and offers a method to evaluate these processes in much greater detail than has been previously possible.

Figure 1 .
Figure 1.Simulated time evolution of species a) N(1,0), b) N(0,1) and c) N(1,1) for a system modeled by the constant kernel.The solid lines are the analytical solutions of the two-component KCE.

Fig. 1 .
Fig. 1.Simulated time evolution of species (a) N (1, 0), (b) N(0, 1) and (c) N(1, 1) for a system modeled by the constant kernel.The solid lines are the analytical solutions of the two-component KCE.

Figure 3 .Fig. 3 .
Figure 3. Discrete two dimensional droplet distribution ( , ) N m n obtained by averaging over 1000 realizations of the Monte Carlo algorithm for the constant kernel at t=100 s.Simulations were conducted with initial conditions (1,0;0) 30 N = cm -3 and

3. 2 Figure 8 .
Figure 8.The number of particles averaged over 1000 realizations and normalized to initial number of particles N0=250 cm -3 (represented by the line with crosses) and the analytical solution of the one dimensional kinetic collection equation (KCE) (represented by the dark solid line) as a function of size for t=50 s.

Fig. 8 .
Fig. 8.The number of particles averaged over 1000 realizations and normalized to initial number of particles N 0 =250 cm −3 (represented by the line with crosses) and the analytical solution of the one dimensional kinetic collection equation (KCE) (represented by the dark solid line) as a function of size for t=50 s.

Fig. 10 .
Fig. 10.Normalized size distributions obtained from analytical solution of one-component KCE for Golovin kernel, versus size distributions from numerical solution of two-component KCE for two times (t=600, 900 s).

Fig. 11 .
Fig. 11.Normalized size distributions obtained from numerical solution of one-component KCE for hydrodynamic kernel, versus size distributions from solution of two-component KCE for two times (t=600, 900 s).
Two-dimensional size distributions N(m,n;t) resulting from a) averaging over 1000 realizations of the Monte Carlo algorithm for the hydrodynamic kernel and b) numerical solution of two-component KCE at t=300 s.

Fig. 13 .Figure 14 .
Fig. 13.Two-dimensional size distributions N (m, n; t) resulting from (a) averaging over 1000 realizations of the Monte Carlo algorithm for the hydrodynamic kernel and (b) numerical solution of two-component KCE at t=300 s.

Fig. 16 .
Fig.16.Liquid water content (LWC, in g/cm 3 ) as a function of drop radius for the hydrodynamic coalescence kernel for t=150 s (dashed line) and t=1500 s (solid line with diamonds) obtained after averaging over 1000 realizations of the Monte Carlo algorithm.

Figure 17 .
Figure17.Aerosol mass concentration (in g/cm 3 ) as a function of aerosol radius for the hydrodynamic kernel for t=150 s (dashed line) and t=1500 s (solid line with diamonds) obtained after averaging over 1000 realizations of the Monte Carlo algorithm.

Fig. 17 .
Fig.17.Aerosol mass concentration (in g/cm 3 ) as a function of aerosol radius for the hydrodynamic kernel for t=150 s (dashed line) and t=1500 s (solid line with diamonds) obtained after averaging over 1000 realizations of the Monte Carlo algorithm.
Figure 10.Normalized size distributions obtained from analytical solution of onecomponent KCE for Golovin kernel, versus size distributions from numerical solution of two-component KCE for two times (t=600, 900 s).