A statistical – numerical aerosol parameterization scheme

A new modal aerosol parameterization scheme, the statistical–numerical aerosol parameterization (SNAP), was developed for studying aerosol processes and aerosol– cloud interactions in regional or global models. SNAP applies statistical fitting on numerical results to generate accurate parameterization formulas without sacrificing details of the growth kernel. Processes considered in SNAP include fundamental aerosol processes as well as processes related to aerosol–cloud interactions. Comparison of SNAP with numerical solutions, analytical solutions, and binned aerosol model simulations showed that the new method performs well, with accuracy higher than that of the high-order numerical quadrature technique, and with much less computation time. The SNAP scheme has been implemented in regional air quality models, producing results very close to those using binned-size schemes or numerical quadrature schemes.


Introduction
Aerosol particles may strongly influence air pollution, cloud and precipitation formation, and climate and environment changes.Key factors that determine the influence of aerosols are their size spectrum and chemical compositions.However, these factors are highly variable and thus can be difficult to simulate in either regional-or global-scale atmospheric chemistry or air pollution models.Moreover, different particulate chemicals may coexist in a specific air parcel by external or internal mixing.These mixing states have additional influence on the physical and optical properties of particles (Chylek and Wong, 1995;Jacobson, 2000).The large number of possible combinations between aerosols of different origins further complicates their roles in atmospheric processes (Jacobson, 2001;Nenes et al., 2002).Therefore, increasingly sophisticated analytical methods are required to fully understand the roles of aerosols in the atmosphere.
Earlier regional models for studying aerosol processes, such as RADM2 or CAMx, keep track of only the aerosol mass concentration.Such bulk methods are insufficient in resolving size-sensitive processes, such as dry and wet deposition, cloud drop activation, light scattering and absorption, and impacts on health.Therefore, an increasing number of models are adapting size-spectrum schemes.Size-spectrum schemes can be incorporated into regional or global aerosol models in different ways.One approach is to use sectionalsize models that categorize the particles into a manageable number of bins according to their sizes (e.g., Gelbard et al., 1980;Wexler et al., 1994;Jacobson, 1997;Russell and Seinfeld, 1998;Yu and Luo, 2009;Bergman et al., 2012).The accuracy of sectional models very much depends on the number of bins applied.Having fewer bins inevitably leads to higher levels of error (Landgrebe and Pratsinis, 1990;Kumar and Ramkrishna, 1996a).Numerical diffusion is a fundamentally challenging problem for the sectional methods when solving the mass transfer among bins.The problem is more serious for the collision-coagulation processes, which need to be handled with advanced numerical techniques (e.g., Drake, 1972;Tzivion et al., 1987;Landgrebe and Pratsinis, 1990;Chen and Lamb, 1994;Kumar and Ramkrishna, 1996b).Also, the growth kernel in each bin is often assumed to be constant; in reality, however, the growth kernel usually is very sensitive to aerosol size and thus may vary significantly between bin limits.Using a large number of bins can reduce the numerical diffusion; at the same time, however, it results in an increase of the computational burden.In particular, the computational time required for particle coagulation processes is proportional to the square of the bin number.Therefore, when computational resources are Published by Copernicus Publications on behalf of the European Geosciences Union.
J.-P.Chen et al.: A statistical-numerical aerosol parameterization scheme limited, sectional schemes are not suitable for regional-or large-scale models.
Another frequently used approach for aerosol simulations is the so-called modal scheme.In typical modal schemes, a complete aerosol size distribution is composed of several modes, and each mode is represented by a relatively simple mathematical function.The evolution of the size distribution is solved by deriving analytical solutions for an integral of the size distribution multiplied by the growth kernel.Computation is less intensive for such modal approaches because the number of variables that need to be tracked is significantly reduced.Zhang et al. (1999) evaluated several air quality models and showed that the modal approach is within reasonable agreement of the sectional model, and requires only about 1% of the CPU time when calculating coagulation.A similar conclusion was reached by Mann et al. (2012), who compared sectional and modal aerosol modules in a global chemical transport model.Because of this, the modal approach has been widely adopted in many current aerosol models (e.g., Seigneur et al., 1986;Whitby et al., 1991;Binkowski and Shankar, 1995;Whitby and Mcmurry, 1997;Ackermann et al., 1998;Harrington and Kreidenweis, 1998;Schell et al., 2001;Wilson et al., 2001;Vignati et al., 2004;Mann et al., 2010;Pringle et al., 2010;Liu et al., 2012).
The main weakness of modal parameterization is that analytical solutions are needed for calculating the evolution of size distribution, but the exact solutions are not always available due to complicated mathematical forms of the growth equations.In such a situation, the growth equation must be simplified to get an analytical solution; however, this simplification can lead to large uncertainties.Therefore, we developed in this study a set of aerosol parameterization methods to provide better accuracy and computation efficiency for aerosol simulations.These methods are applied to parameterize microphysical processes -such as ice nucleation, condensation, coagulation, and sedimentation; they are also used to provide diagnostic equations, such as the Kelvin effect on aerosol wet size.

Methodology
The basic concept behind our new approach is to perform offline numerical integration over the aerosol size spectrum for each aerosol process.The numerical integration for each individual process is performed under specified conditions that cover all possible variations in atmospheric states and aerosol size modal parameters.Properties and conversion rates for each aerosol mode obtained from the numerical integration are then analyzed statistically and fitted into so-called modal (or bulk) formulas.

Size distribution function and its moments
The first step of our modal approach is to select a mathematical function that best represents the number density distribution of each modal population.Observational results showed that aerosol size distribution can generally be represented well by the multimode lognormal function (Whitby, 1978); several studies have indicated that such a distribution is selfpreserving (Friedlander, 1960;Hidy, 1965;Liu and Whitby, 1968;Lai et al., 1972).Therefore, we select the lognormal function to represent each modal distribution: where n is the number density distribution function, r is the particle size, N is the total number of particles, σ is the standard deviation (in the ln r coordinate), and µ is the modal radius.The whole aerosol size distribution may be composed of several such modal functions.The lognormal distribution requires three parameters for description: N, σ , and µ.However, these modal parameters are not extensive properties and thus cannot be used as prognostic variables in atmospheric models.In practice, the desirable tracking variables are the moments of the size distribution, such as the zeroth moment (i.e., number concentration) and third moment (i.e., volume or mass concentration).The kth moment is defined as (2) For n(r) in the lognormal form, an analytical solution for Eq. ( 2) can be solved as The zeroth and third moments are logical choices for tracking variables because of their direct relevance to many physical properties.Yet the selection of the next moment is optional.For example, in cloud microphysical parameterization, Milbrandt and Yau (2005) used the zeroth, third, and sixth moments.The sixth moment represents the radar reflectivity, which is an important characteristic of large precipitation particles.Binkowski and Shankar (1995) (hereafter BS95) also selected the sixth moment for their aerosol parameterization because it allows for easier derivation of analytical solutions.However, the cross-sectional area, represented by the second moment, is important to light scattering and atmospheric radiation and is consequently more relevant to aerosol studies.Thus, we select the second moment as the third tracking variable for this study.Note that the current modal aerosol module in USEPA Models-3 Community Multiscale Air Quality (CMAQ) model, although based on BS95, does not track the sixth moment but instead considers the second moment (Binkowski and Roselle, 2003).Also note that this CMAQ model will be used to test our scheme in Sect.4.2.
To further reduce computation time, some of the modal aerosol models (such as in NCAR CAM5) actually applied only two prognostic variables.These models typically keep track of the changes in number and mass moments, but use fixed spectral widths (σ ) for the lognormal size distribution.Mann et al. (2012) found that such two-moment modal module may produce strong bias in aerosol.They also showed that the choice of σ can have significant impacts on the model results.Thus, including a third variable is important in achieving high model accuracy.
The size distribution parameters in Eq. ( 1) can be diagnosed from the three moments as which can then be used to calculate the modal size: . (5) Note that the methodology shown in the next section is not restricted to the lognormal size distribution.It can also be applied to the gamma-type distribution functions, which are mathematically and numerically attractive for the representation of particle size spectrums.But in this study we focus on the lognormal distribution.

Parameterization methods
After the mathematical form and the key parameters of the size distribution are determined, the evolution of size distribution can be described in terms of the rate change of the moments: where K k is the growth kernel for the kth moment.This growth kernel represents the fundamental growth equation for each process.A few examples of the growth kernel will be discussed in detail in Sect.3. When the growth kernel is not in a simple form, solving such integrals requires computationally intensive numerical techniques, such as Gauss-Legendre or Gauss-Hermite numerical quadrature.Therefore, parameterization of Eq. ( 6), which enables the efficient and accurate calculation of aerosol and cloud microphysical processes, is desirable for many meteorological and air pollution models.
Common treatments of Eq. ( 6) include the use of lookup tables and kernel simplification.The lookup table approach calculates the kernel or the whole integral as a function of their key parameters and then arranges the results in tables that, when applied in models, can be searched according to the current values of those parameters.This method has the advantage of fast calculation, as it primarily involves searching, and has high accuracy when the tables are large enough.Some sectional models also applied lookup tables to reduce computational costs (e.g., Yu and Luo, 2009).However, the method may become cumbersome to use when the process involves too many parameters that require large table dimensions.In addition, the lookup table method usually cannot be used directly for physical interpretation or analysis of the functional dependence on key parameters.Alternatively, the kernel simplification approach is commonly applied in the parameterization of both aerosol and cloud microphysics.Its specific purpose is to allow for easy evaluation of Eq. ( 6) into analytical solutions.However, such simplifications are often too rough and can result in large errors.
We investigated four methods of parameterization: (A) mean-size approximation, (B) kernel transformation, (C) integral transformation, and (D) optimal-size approximation.The mean-size approximation approach can be considered as a no-skill method.We will show that the other three methods are significantly more accurate and will be further selected for our final parameterization based on the accuracy of the analyses.Since the last three methods apply statistical fitting on numerically integrated results, our overall method is named the statistical-numerical aerosol parameterization (SNAP).

Mean-size approximation method
Mean-size approximation (hereafter called MSA) is achieved by replacing all or some of the size variable r in the growth kernel with a constant size r so that the kernel, or part of the kernel, can be taken out of the integral in Eq. ( 6).It is mathematically possible to approximate the growth kernel K(of any moment) by a polynomial function of r with sufficient number of terms, i.e., K = i a i r i .We apply such a polynomial function here just to demonstrate the error associated with MSA.The corresponding growth rate for each term of order i (neglecting the coefficient a i ) can be written as This equation has an exact solution M i as given earlier in Eq. ( 3).On the other hand, the mean size approximation is Several forms of the mean size r can be used for MSA.A group of these forms is called the moment-weighted mean size r n ≡ M n M 0 1/n .For example, r 2 and r 3 are the surface-and volume-weighted mean sizes, respectively.According to Eq. ( 3), r n can be converted to Let us use ˜ i,r n to represent the approximate solution using these nth-moment-weighted sizes.Its ratio to the exact solution i can be derived as Other forms of the mean size include the modal size µ in Eq. ( 1) and the effective radius r e ≡ M 3 /M 2 , which is commonly used for radiation budget calculation.Ratios of the solution using these two mean-size approximations to the exact solution can be derived as The approximations using µ and r e are special cases of Eq. ( 10), with n = 0 and n = 5, respectively.Thus, µ and r e may be called the zeroth-and fifth-moment-weighted sizes, respectively.Figure 1 shows the errors associated with these mean-size approximations, which exhibit the following features: (1) the error increase with the width of the size spectrum (i.e., σ ), the order of the kernel (i.e., i), and the difference between n and i (i.e., |n − i|) in Eq. ( 10).Therefore, the error can be minimized if n is set as equal to i.
(2) The error is positive for n > i and negative for n < i.This indicates that the signs of error may be opposite for the growth of different moments.For growth kernels containing several polynomial terms, it would be best to select n that lies between the orders of all dominating terms, such that their errors may cancel each other.

Kernel transformation
A complicated growth kernel prohibits the derivation of an analytical solution for Eq. ( 6).However, it is possible to transform such kernels into manageable mathematical forms.
We call this approach SNAP-KT.For a lognormal n(r), useful mathematical forms include the power-law function r a , the exponential function exp(b ln 2 r), or their combinations.The conversion of growth kernels into such functional forms is done by statistical fitting of the numerically solved results.Some examples will be given in the next section.These fitting functions can be generalized as r a exp(b ln 2 r), which can also be expressed as exp(aln r + b ln 2 r).This allows Eq. ( 6) to be expressed as Its solution can be derived by introducing the variable exchange x≡αln r−γ , where α≡ 1 2σ 2 − b, β ≡ a + k+ lnµ σ 2 , and γ ≡ β 2α .We then have One can verify that Eq. ( 14) reduces to Eq. ( 3) when a = b = 0.In other words, Eq. ( 3) is the special case of F (k, 0).

Integral transformation
SNAP-KT formulations, such as Eq. ( 14), are computationally efficient.Yet satisfactory fitting of the growth kernel, as discussed above, is not always available.When this is the case, we can turn to the integral transformation method (hereafter called SNAP-IT), which involves two steps: (1) solving Eq. ( 6) numerically by discretizing the size spectrum into fine bins (as fine as possible) for a wide range of ambient conditions and size spectrum parameters (e.g., µ and σ ); and (2) analyzing the results by statistical fitting to obtain a transformed formula.However, a technical problem may arise while performing the fitting.Besides the three moments, the growth equation often contains other dependent variables, such as air temperature and pressure.Few statistical software packages can handle nonlinear fitting on multiple variables.For example, the commercial software we are using can handle only two variables at a time.Processing all of the variables may require intensive trial by error or iteration before a satisfactory parameterization formula can be acquired.Consequently, a conversion of the growth kernel for the purpose of variable separation before performing the numerical integration may be necessary.However, such variable separation is not always easy, and this greatly limits the application of this approach.We overcome this deficiency by taking advantage of the MSA method in which the dependence on ambient parameters is largely retained in the simplified kernel.We obtain SNAP-IT first by rewriting Eq. ( 6) as where Ĩk is the modal-value approximation (cf.Eq. 11) of I k , and g 1 is a correction factor that brings Ĩk closer to I k .The corrector g 1 should depend strongly on the spectral width σ because Ĩk is calculated by assuming a monodisperse size distribution (and thus σ = 0).We derive g 1 by integrating Eq. ( 6) numerically for a range of σ , as well as other size distribution parameters and ambient parameters to obtain the "true" value of I k .Each I k value is then divided by the MSA value Ĩk , and their ratios are fitted to obtain g 1 as a function of σ and other parameters.In this way, the ambient-parameter dependence is largely retained in Ĩk , while the dependence on the spectral width σ is largely contained in g 1 .Note that some computational efficiency is lost by keeping the details 1  of the growth kernel in Ĩk , as compared with a direct integral transformation (i.e., without utilizing MSA).This loss in computational efficiency is well compensated by the accuracy that is gained.

Optimal-size approximation
In the MSA approach, we assume that I k ∼ = Ĩk (µ), and in SNAP-IT, we find a correction factor to improve this approximation.The deviation of Ĩk (µ) from I k indicates that the modal value µ (or any other mean size) may not be the best representative size.In fact, we showed in Eq. ( 10) that this "best size" is actually a function of the order of the kernel and spectral width σ , and potentially some ambient parameters as well.Thus, instead of using a specific mean size (i.e., µ) and then correcting the whole integral with g 1 , as done in SNAP-IT, it may be possible to find in advance an optimal mean size, which can be adjusted with the imposed conditions to provide an accurate value of Ĩk directly according to the following relationship: To determine the formula for the optimal size µ for this SNAP-OS method, we first calculate I k for a range of relevant parameters.For each I k value, we search by iteration for a value of µ that, when placed into Ĩk , gives an exact value of I k .Afterward, the ratios of µ to µ (i.e., g 2 ) under various conditions are analyzed statistically to fit into a function of the key parameters, such as σ or µ.
The SNAP methods can be summarized as follows: (1) SNAP-KT: kernel transformation to obtain a semianalytical solution for the integral; (2) SNAP-IT: integral transformation that provides a modification factor to the MSA method; and (3) SNAP-OS: parameterization of optimal size that replaces the constant size in MSA.The MSA method is taken as a benchmark, and we will demonstrate that the SNAP parameterization methods are all significantly more accurate and thus have high skills.

Parameterization of microphysical processes
In this section, we apply the above methods to various aerosol microphysical processes and analyze the parameterization accuracy by comparison with the numerical solutions.The numerical solutions for I k are obtained by discretizing the size spectrum with 10 bins per decade, and then summing the rates from individual bins.Higher bin resolutions are also tested.Figure 2 shows the dependence of precision on bin resolution using 100 bins per decade as a reference.The example given is for intramodal coagulation, which will be mentioned in Section 3.4.One can see that the error decreases by over two orders of magnitude for an order increase in bin number.The difference between 10-bin and 100-bin calculations is less than 0.5 %, which can be regarded as the precision of the numerical solutions.For noncollisional processes the error is generally smaller as their fundamental equations contain only a single integral.In this study, the error is defined as abs exp{ j i=1 [abs(ln Ĩk /I k )]/j } − 1 , where j is the number of conditions selected for evaluation.

Ice nucleation
Heterogeneous ice nucleation from insoluble aerosol particles (which are thus called ice nuclei) such as mineral dust, soot or bio-aerosols is an important factor in the glaciation of clouds.This process is usually not considered in traditional aerosol models, which do not emphasize aerosol-cloud interactions.On the other hand, current cloud models generally do not consider the emission and production of aerosol particles, so the ice nucleation process is highly parameterized due to the lack of realistic ice nuclei.Because of the importance in climate and hydrological cycle, detailed aerosol-cloud interactions have become an essential component in advanced regional and global models, for which ice nucleation is a critical mechanism that badly needs improvement (cf.Tao et al., 2012).According to the classical theory, the heterogeneous ice nucleation rate can be generalized into the following form for several pathways of nucleation (cf.Chen et al., 2008): where r is the radius of the ice nuclei, A is a parameter that depends on the ambient conditions only, f is a sizedependent geometric factor, g a is the activation energy, g g is the homogeneous germ formation energy, and k B is the Boltzmann constant.The overall nucleation rate for a population of ice nuclei is then expressed as which represents the rate of decrease in ice nuclei concentration due to conversion into cloud ice.This integral cannot be solved analytically, as the geometric factor f , which appears twice in the kernel J H N , has a very complicated form: where m ≡ cos (θ), θ is the contact angle, q ≡ r r g is the ratio of the nuclei size to the nucleation-germ size, and φ ≡ 1−2mq + q 2 .There are several pathways of heterogeneous ice nucleation.Here, we take the immersion freezing nucleation as an example.Its key parameters include temperature and saturation vapor pressure over water (with solute and curvature effects) of the supercooled droplet wherein the ice nuclei are immersed.
Applying MSA to Eq. ( 18) is straightforward: where JHN is Eq. ( 17) calculated with r replaced by the modal size µ.One may also keep the prefactor r 2 of J H N and the r k term staying in the integral to get For the parameterization using SNAP-KT, the parameter f in Eq. ( 17) should be transformed into functions like r a or exp(bln 2 r) in order to derive a semianalytical solution for Eq. ( 18).The following is a readily available formula from Chen et al. (2008): where a 1 , a 2 , and a 3 are constants.This formula is suitable for converting the first term that contains f in Eq. ( 17) into where a 4 ≡ exp(a 1 + a 2 ln (1 − m) − a 3 ln r g ) is independent of r.However, this formula is not useful for simplifying f in the exponential term.Thus, we produced another transformation formula: where , and c 2 = Bb 2 exp(b 3 ) are all independent of r.The R 2 of fitting for Eqs.( 21) and ( 23) both reached 0.9998 for θ in the range of 1 to 110 • and q from 10 to 400; it could be more accurate  if the ranges were divided into a few sectors, each with its own fitting coefficients.With Eqs. ( 22) and ( 24), the overall nucleation rates for a spectrum of ice nuclei can be derived as For SNAP-IT, we first perform numerical integration on Eq. ( 18) and then compare the results with the modal approximation Eq. ( 20) to obtain a fitting on g 1 .The selection of the fitting parameters is not a trivial task.Hints of the proper parameters may emerge while examining the fundamental physics and its mathematical formulation.For example, one may recognize that q in Eq. ( 19) is the most pertinent parameter for heterogeneous ice nucleation.On the other hand, Eqs.
(3) and ( 10) indicate that the variance σ 2 is a key to the representation of size spectrum.Thus, we selected q≡µ/r g and σ 2 as statistical fitting parameters.This indeed results in one of the better fitting formulas: Similarly, the optimal size correction factor g 2 for SNAP-OS can be derived as Figure 3 shows that these two formulas provide reasonably good fittings.It also reveals that large corrections are necessary when q is small and, at least in the case of g 2 , σ is large.Note that there are numerous fitting formulas for our selection, and we often select those that are easier to use and can reflect physical meanings but are not of the highest accuracy.For example, in addition to maintaining the "exp(σ 2 )" dependence, Eq. ( 26) was selected to warrant a unit value toward the extreme conditions of σ → 0 and q→ ∞.
Next, we compare the four parameterization approaches (MSA and SNAPs) against the detailed numerical solution.The results for immersion freezing are shown in Fig. 4, for which the ranges of values tested are the following: 6 modal sizes (µ) between 0.02 and 4.0 µm, 10 modal widths (σ ) between 0.26 and 0.95, 8 temperatures between −5 and −40 • C, and 4 water activities between 0.82 and 1.0.The mean errors in I 0 are 317 % for MSA, 22% for SNAP-KT, 63 % for SNAP-IT, and 16 % for SNAP-OS.These errors tend to increase toward higher moments.They are 2800, 25, 60, and 12 %, respectively, for I 2 , and are 15 100, 34, 73, and 22 %, respectively, for I 3 .One can see that the SNAP-OS method performed significantly better than the other methods do.These errors seem to be large even for SNAP-OS.Fortunately, large deviations usually occur when the absolute values are close to negligible.The CPU time required for SNAP-KT, SNAP-IT, and SNAP-OS are 73, 26, and 18 % more than for the MSA method, respectively.Note that there exist feather-like features in Fig. 4, and each filament represents a set of values with different σ values.In the left panel of Fig. 4 we highlighted the MSA points with the largest σ values with filled circles.One can see that the largest error is associated with the highest σ , and the error approaches zero for a monodisperse distribution (i.e., very small σ ).Using the above example, we demonstrated the details of all SNAP methods.We will omit similar details when discussing the parameterization for other processes.

Gravitational sedimentation
The gravitational sedimentation velocity takes the form of    where V Stokes is the Stokes' law fall speed, is the Cunningham slip-flow correction, g is the normal gravitational acceleration, ρ p is the particle density, η is the dynamic viscosity of air, K N ≡ λ/r is the Knudsen number, and λ is the mean free path of air molecules.Note that C C may take a form somewhat different from Eq. (28) (cf.Seinfeld and Pandis, 2006, p. 407), but our parameterization procedure works the same with both forms.Sedimentation flux for the whole size distribution (also termed the group sedimentation flux) is expressed as As the analytical solution for this equation cannot be readily obtained, BS95 simply ignores the exponential term in Eq. ( 28) to reach the following solution: Under standard atmospheric conditions, omitting the exponential term in C C would cause an underestimation in sedimentation speed by 4 and 26 % for particles of 0.1 and 0.01 µm radii, respectively (Fig. 5a).Such underestimations actually contribute to a small absolute error in the group sedimentation flux; the percentage error is significant only for small particles, whose gravitational fall speed is low.However, an accurate description of C C may still be important for other calculations.For example, C C is an important parameter in the Brownian coagulation kernel (see Sect. 3.4).
If one wants to consider the exponential term for better accuracy, we can apply SNAP-KT by calculating C C for a realistic range of K N and then apply statistical fitting of the results using commercially available software.For example, after calculating C C for a range of K N values, their relation can be curve-fitted into the following: where a 1 = 1.43089 and a 2 = 1.0295 are the fitting coefficients.From Fig. 5a, one can see that the above fitting is quite accurate, with less than 5% error (R 2 of fitting = 0.9999) for all relevant values of K N .Adding more terms to Eq. ( 31) may give even higher accuracy, but is not necessary for practical purposes.This transformation allows Eq. ( 29) to be evaluated analytically as Note that Whitby et al. (1991) applied a similar transformation but used different a 2 values for different K N regimes to gain better accuracy.Figure 5b shows the comparison of gravitational sedimentation parameterizations.One can see that Eq. ( 32) gives good results comparing to the exact solution, whereas Eq. (30) (i.e., BS95) produces large error at small values.As SNAP-KT can already produce very good results, we will omit applying SNAP-IT and SNAP-OS to the gravitational sedimentation.

Condensation
Under the assumption of a steady-state diffusion process, the kernel of condensation growth following the two-stream Maxwellian kinetic theory with a steady-state assumption is commonly expressed as (cf.Pruppacher and Klett, 1997, p. 506) where D is the diffusion coefficient; f g is the modification due to the gas kinetic effect (Fuchs, 1959(Fuchs, , 1964)); f v is the ventilation coefficient, which can be ignored for small aerosol particles; ρ v,∞ is the ambient vapor density; and ρ v,p is the surface vapor density.The parameters D, ρ v,∞ , f g , and ρ v,p are species dependent, whereas f g and ρ v,p are also size dependent.Furthermore, ρ v,p is influenced by latent heating/cooling during condensation/evaporation.A quasianalytical solution can be obtained to account for this effect (cf.Pruppacher and Klett, 1997, p. 511), but the details will not be elaborated here.Equation ( 33) can be generalized for the simultaneous condensation of multiple species.Let the volume change due to condensation be dv = dm/ρ L , where v = 4π r 3 /3 and ρ L is the density of the condensate.From this, the bulk growth rate of the kth moments can be expressed as Note that in this formula the rate change of the total number (k = 0) for the condensation process necessarily equals zero.
If we assume that Df g ρ v,∞ − ρ v,p is size-independent, then an analytical solution can be easily derived as However, in reality the size dependence of f g and ρ v,p cannot be ignored, particularly for small aerosol particles.The Kelvin effect on ρ v,p will be further discussed in later sections.Here, we focus on the parameter for the surface gaskinetic effect, which is generally expressed as where is the vapor jump distance and is on a scale similar to that of the mean free path λ, α is the mass accommodation coefficient, and ν is the mean thermal velocity of the gas molecules (cf.Fuchs, 1959;Pruppacher and Klett, 1997).Considering the dependence of on λ, Fuchs and Sutugin (1970) provided an empirical formula for f g as a function of K N and α: It is difficult to arrive at analytical solutions to Eq. ( 34) with the formulas for f g given in Eqs. ( 36) and ( 37).An approach to resolving such a problem, as suggested by Pratsinis (1988) and adopted by the BS95 method, is to consider the harmonic mean of growth in the free-molecular regime and continuum regime: where I M,k is calculated with the free-molecular regime growth kernel K M = πr 2 αν ρ v,∞ − ρ v,p , and I C,k with the continuum regime kernel K C = 4πrD ρ v,∞ − ρ v,p .Since I M,k and I C,k can be solved analytically as a function of M k and M k−1 , respectively, Eq. ( 34) can be evaluated analytically.Although Pratsinis (1988) indicated that the harmonic mean can approximate the results well using f g from the equation developed by Fuchs and Sutugin (1970) shown in Eq. ( 37), it inevitably contains some inaccuracy, which we will evaluate below.Fukuta and Walter (1970) suggested a slightly different form of f g , which, in effect, excludes the term in Eq. ( 36), and is, for practical purposes, a harmonic mean of K M and K C : Below we omit the application of the SNAP-KT method because the fitting formula becomes too cumbersome for practical purposes.Additionally, we omit the SNAP-OS method because the SNAP-IT method is sufficient.The SNAP-IT fitting formula that we derived is as follows: where K N ≡ λ/µ represents a mean Knudsen number.Figure 6 shows the comparison between various parameterization methods for the condensation growth process.Note that the number concentration does not change during the condensation process (i.e., I 0 = 0), so only I 2 and I 3 are presented.MSA gives good results only when σ is small, but it is biased toward lower values for increasingly larger σ values (i.e., the true value increases with σ , but MSA does not).The overall error for MSA is 17 % in I 2 and 92 % in I 3 .SNAP-IT performed rather well, with 0.74 and 1.3 % error in I 2 and I 3 , respectively.The BS95 method produced significantly larger discrepancies, with 10.7 and 57.1 % error in M 2 and M 3 , respectively.However, the BS95 computation time is 21 % less than that of the SNAP-IT method.
In Fig. 6 we also plotted the numerical solutions using f g from Fukuta and Walter (1970).The strong positive biases (around 83 %) indicate a significant error associated with the harmonic mean approximation.

Brownian coagulation
Calculation of the rate change of moments caused by collision-coagulation processes involves double integrals over the size spectra of the two aerosol modes involved.For coagulation between two particles of sizes r A and r B , the coagulated particle has a size r C = (r 3 A + r 3 B ) 1/3 .It follows that the changes in their kth moments are −r k A and −r k B , respectively, for each original particle, and +r k C for the coagulated particle.With these parameters defined, the fundamental equation for coagulation between particles in the collector mode A and the contributor mode B can be expressed as where the coagulation kernel K is usually a nonlinear function of the two particle sizes and environment properties denoted by the parameter C air .Note that the coagulated particle is placed back into mode A as indicated in Eq. ( 40).
In these generalized equations one can easily verify that the number concentration (M 0 ) in the collector mode remains unchanged (i.e., I 0,A = 0) and that the total volume is conserved (i.e., I 3,A = −I 3,B ).Hence, a total of four conversion rates are needed, i.e., I 0,B , I 2,A , I 2,B , and I 3,A (or −I 3,B ).For the intramodal coagulation (i.e., A = B), the number of rates reduces to two, and all coagulation rates in Eqs. ( 41) and ( 42) should be divided by 2 to correct for double counting.
Processes contributing to aerosol coagulation include Brownian diffusion, convective Brownian diffusion enhancement, gravitational collection, turbulent inertial motion, and turbulent shear flow (Jacobson, 1997).Brownian diffusion is the dominant coagulation process for fine aerosol particles with radii typically in the range 0.01-1 µm.Here, we take this most complicated kernel as an example for parameterization, starting with the intramodal coagulation, which involves only its own moments.By analogy of gas diffusion formulation, Fuchs (1959) expressed the Brownian coagulation kernel K Br between particles A and B as where r = (r A +r B )/2, D p is the mean particle diffusion coefficient, and β represents the modification due to concentration discontinuity near the surface of the receiving particle.The mean particle diffusion coefficient is defined as D p = (D p,A +D p,B )/2, where D p,i = k B T C C,i 6π ηr i ; C C is the Cunningham slip flow correction factor, as shown in Eq. ( 28); k B is the Boltzmann constant; and η is the dynamic viscosity of air.The conventional form of β is where α p is the sticking probability (usually assumed to be unity) when two particles collide, − 2r i represents a mean coagulation distance, λ p,i = 2D p,i π ν p,i is the mean free path of the particle, and i is either A or B. The factor β has a similar form to  Eq. ( 36).However, the variables that it contains -namely δ, D p , and ν p -are all complex functions of the particles' sizes, and this makes the SNAP-KT method unfeasible to use.For this coagulation process, Pratsinis (1988) applied the harmonic-mean approximation.This approximation was also applied in the BS95 method: where I Br,M and I Br,C are the results with kernels K Br,M = 2πr 2 α p ν p and K Br,C = 8πrD p , respectively.However, the complex forms of C C , ν p , and δ still prevent the derivation of analytical solutions for I Br,M and I Br,C .Thus, following Whitby et al. (1991), BS95 made a few algebraic manipulations and combined them with lookup tables to solve the harmonic mean.For a similar reason, our parameterization for Brownian coagulation focuses on MSA and SNAP-IT but ignores SNAP-KT and SNAP-OS.There is a complication in using MSA here, because the two modal sizes used for the calculation are the same for intramodal coagulation.We found it helpful to offset the modal radius and assign r A = µ • σ 2 and r B = µ/σ 2 in Eqs. ( 40) and ( 41) for calculating Ĩk in Eq. ( 15).With this treatment, the correction factor for I 0 is obtained as www.atmos-chem-phys.net/13/10483/2013/J.-P.Chen et al.: A statistical-numerical aerosol parameterization scheme which is used further to get the correction factor for I 2 g 1,2 = g 1,0 • (a 1 + a 2 ln µ + a 3 σ 3 ). (47) Figure 7 shows the results using MSA and SNAP-IT for intramodal coagulations.Also compared is the harmonicmean approximation of BS95, as well as the numerical solutions calculated with the fifth-order Gauss-Hermite quadrature (GHQ), which is an accurate but computationally expensive option in the CMAQ model.Note that the amount of data for Brownian coagulation is much larger than that for the previous processes, so only a selected amount of data (e.g., 1 out of every 5 or 10 consecutive points) is shown to avoid clutter.One can see that BS95, GHQ, and SNAP-IT all perform reasonably well.SNAP-IT produces 3.7 and 5.9 % errors in I 0 and I 2 , respectively, which are similar to those in GHQ (4.5 and 4.0 %).The error in BS95 is about the same in I 0 (4.5 %) but somewhat larger in I 2 (22 %).The computation time used for SNAP-IT and BS95 are 12 and 10 % of that for GHQ, respectively.The intermodal Brownian coagulation involves two size distributions, so one would imagine its parameterization must be more complicated than the intramodal coagulation.However, using the SNAP-IT method, we found a rather simple but accurate formula for all intermodal rates: It turns out that g 1 mainly depends on the two spectral widths (i.e., σ A and σ B ), whereas the effects of other parameters, such as λ, have been largely reflected in the modal mean, Ĩk , and thus play little role in g 1 .Also, this formula agrees with the exp σ 2 dependence shown in Eq. ( 10).The two coefficients vary with the moments (i.e., the k value), but a 1 is consistently much smaller than a 2 (see Table A2), indicating that intermodal coagulation is more sensitive to the spectral width of the contributor mode (σ B ) than that of the collector mode (σ A ).
Figure 8 shows the accuracy of various evaluation methods for these rates.The MSA method again deviates from the numerical solution more pronouncedly at larger σ , and the mean error ranges from 18.1 to 74.1 % for various moments.The SNAP-IT method is rather accurate, having errors ranging from 2.6 to 4.5 % for the four conversion rates, which are a little better than the errors of 4.8 to 5.4 % produced by GHQ, and 4.8 to 7.4 % produced by BS95.The computation time required for SNAP-IT and BS95 are 7.8 and 7.0 %, respectively, of that for GHQ.

Other processes and diagnostic parameters
A rate process that has not been discussed earlier is aerosol scavenging by cloud drops or raindrops, which is also a type of intermodal coagulation.The mechanisms that control aerosol scavenging include Brownian diffusion, collection by phoretic forces, and gravitational collection.For the two former mechanisms, Wang et al. (1978) provided a mathematical solution that combines the two kernels, which is adopted for our parameterization.For the gravitational collection, we used the kernel from Slinn (1977).Parameterization procedures for these processes are quite similar to that for the Brownian coagulation, so only the final results are listed in Appendices A and B.
In earlier examples we showed that SNAP-IT and SNAP-OS generally outperform SNAP-KT.However, for diagnostics parameters that do not involve spectral integration obviously cannot utilize SNAP-IT or SNAP-OS.For these parameters SNAP-KT comes in handy.In fact, we have already shown parameterization on the parameter C C , which is used to derive the group sedimentation velocity, V sed , in Sect.3.2.In the below we demonstrate the application of SNAP-KT to other diagnostic parameters.
An important microphysical process that does not directly involve existing aerosol spectra is aerosol nucleation (new aerosol production).The mechanisms that control aerosol nucleation include homogeneous binary or ternary nucleation (Nair and Vohra, 1975;Coffman and Hegg, 1995) and ion-enhanced nucleation (Yu, 2006).Here we discuss the parameterization on homogeneous binary nucleation from water and sulfuric vapors as an example.The rate of binary nucleation depends mainly on the temperature and the saturation ratios of water vapor and sulfuric vapor.We will not focus on the details of the binary nucleation rates, which can be found in textbooks such as that written by Seinfeld and Pandis (2006).Instead, we will focus on a key parameter that needs to be solved by iteration: the water-sulfuric acid mixing proportion in the critical embryo.Once this parameter is obtained, the calculation of nucleation rate is straightforward.In brief, we precalculated this mixture fraction numerically for various ambient conditions and then fit the results into certain formulas, as was done earlier using the SNAP-KT methods.By applying this formula, the time required for iteration can be saved.A similar approach was applied by Kulmala et al. (1998) and Vehkamäki et al. (2002) to obtain aerosol nucleation rates.Note that, although some studies suggest that the classical binary nucleation rate may be too weak to explain observed new particle formation (e.g., Covert et al., 1992), Chen et al. (2011) indicated that earlier studies may have significantly underestimated the nucleation rates because they omitted the size dependence of surface tension.Therefore, for the binary nucleation formula given in Table A1, we adopted the method of Chen et al. (2011) for calculating the rate parameters.Another example of diagnostic parameter is the Kelvin effect, which affects the equilibrium vapor pressure of the droplet.The equilibrium radius r eq , and thus the water content of a hygroscopic particle, can be described by the Köhler theory, which is a combination of the Raoult (or solute) effect and Kelvin (or curvature) effect.With the Kelvin effect, the particles absorb less water and thus have smaller sizes (Fig. 9).The size difference due to the Kelvin effect increases with humidity, reaching about 50 % at 95 % relative humidity and near infinity as the relative humidity approaches 100 % for the case shown in Fig. 9. Apparently this effect cannot be ignored, especially for high-humidity conditions.Yet many modal aerosol models considered only water uptake due to the Raoult effect (e.g., Jacobson, 1997;Mann et al., 2010).A few that did take the Kelvin effect into account (e.g., Ghan et al., 2001) need to utilize a convenient form of the Kelvin equation that is applicable only for sufficiently dilute solutions.Normally, rigorous calculation of r eq requires numerical iteration.Here, we apply the SNAP approach to parameterize r eq as a function of the ambient humidity and tem-perature, particle dry size, and a kappa parameter that was introduced by Petters and Kreidenweis (2007) to represent particle composition.Note that for aerosol mixtures (soluble or insoluble), the overall kappa parameter can be obtained by a volume weighting of individual kappa parameters.A similar formula is obtained for calculating the wet volume of a whole aerosol mode.See Table A1 for the details of these formulas.
Another useful diagnostic parameter related to the Köhler curve is the activation cutoff size, which determines the smallest aerosol particles that can be activated into cloud drops under a certain supersaturation.Exact calculation of this cutoff size is even more tedious than obtaining r eq .Hence, it is often derived by simplifying the Köhler equation to obtain an approximate but direct relationship between the cutoff size and ambient supersaturation (cf.Pruppacher and Klett, 1997, p. 178).Our SNAP approach is well suited for parameterizing the cutoff size with high accuracy (<0.5 % error) in a way similar to that for obtaining r eq .As given in Table A1, the cutoff size is expressed as a function of the

Kohler Raoult
Köhler Raoult Fig. 9. Relationship between the ambient relative humidity and equilibrium size for an ammonium sulfate particle with 0.01 µm dry radius.The red dashed curve is the Köhler curve, which includes both the Kelvin effect and Raout effect, whereas the blue curve considers only the Raoult effect.The grey dashed line indicates 100 % relative humidity.
supersaturation, temperature, particle dry size, and the kappa parameter.
Other diagnostic parameters that we have provided in Table A1 include the modal extinction coefficient and absorption coefficient, which are important for calculating aerosol radiation effects.Another important parameter for radiation calculation is the effective radius that, under the modal assumption, has an analytical solution r e ≡ M 3 /M 2 and thus does not need parameterization.Coefficients for the parameterization formulas in Table A1 are given in Table A2.

Numerical verifications
In the previous section, we obtained fairly accurate modaltype parameterizations for aerosol microphysical processes.Additional checking of the reliability of these formulas is necessary when performing time integration, as errors may accumulate with time, which could cause numerical instability in extreme cases.

Verification with the binned parcel model
Verification of the time evolution of the size spectrum is not an easy task, especially for collision processes.A commonly accepted verification method is to use a detailed bin model that truly resolves the size distribution.The binned aerosol model used in this study is modified from the detailed cloud microphysical model of Chen and Lamb (1994), which applies a moment-conserving numerical scheme that ensures accuracy and conservation of mass and number concentra-tion.This model has been applied to various aerosol studies (cf.Chen et al., 2011;Tsai et al., 2012).
Another verification method is to obtain analytical solutions for the spectral time evolution.Such analytical solutions exist for simple collision kernels, such as the constant kernel (Bleck, 1970) and the simple mass-dependent Golovin kernel (Berry, 1967), which have been used in verifying cloud microphysics schemes (e.g., Berry, 1967;Tzivion et al., 1987;Chen and Lamb, 1994).However, there is no need to develop modal parameterization for these simple kernels because exact analytical solutions exist.Thus, timeevolving analytical solutions are typically used to verify the performance of bin models.The performance of the model we are using has been verified against these time-evolving analytical solutions for cloud microphysical processes (Chen and Lamb, 1994).We reconducted the verification for aerosol size scales and found that the bin model acquired similar high accuracy and conservation of the moments.Note that the analytical solutions mentioned above are for gamma-type size distributions.For the lognormal size distribution that we applied here, Park and Lee (2000) provided an analytical solution for constant kernel collision process.Hence, we conducted an additional verification by comparing with their analytical solution for a lognormal size distribution.The bin model produced 0.1 and 0.3 % errors in M 0 and M 2 , respectively, after a 12 h time integration.These smaller errors indicate the robustness of our bin model.
We selected Brownian coagulation (including intramodal and intermodal) for testing the time integration for its complexity.The simulations were run in parcel mode to avoid complications from other processes, such as transport and sedimentation.Results obtained using the GHQ and BS95 methods were also compared.Figure 10 shows the initial bimodal aerosol size distribution (nucleation mode and accumulation mode) and the evolved size distributions.The size distributions of the modal approaches (i.e., BS95, GHQ, and SNAP-IT) are retrieved from the three moments by assuming lognormal distribution for each mode.All modal calculations give results similar to those of the binned calculation, showing that the nucleation mode decreased significantly after one hour and essentially disappeared after six hours, whereas the accumulation mode evolved rather slowly.When looking into the details, one can find visible differences between the modal distributions and the binned calculation.For example, the BS95 and GHQ distributions deviate more obviously at the small end at one hour as well as at the large end at six hours, whereas the SNAP-IT distribution deviates more at the small end at six hours.All modal methods show fewer particles at the larger end of the accumulation mode, especially for the BS95 and GHQ methods and for the higher moments.However, such differences are not totally due to the inaccuracy of the parameterization formulas.The modal approaches retrieve the size distribution by assuming a fixed lognormal shape, which is symmetrical about the mode.However, the  binned solution indicates that the true shape is not perfectly symmetrical.
A more appropriate comparison is done by examining the evolution of the overall moments M 0 and M 2 (while M 3 is conserved).As shown in Fig. 11, M 0 of all modal calculations closely follows the binned results, with errors of 1.8, 2.1, and 2.1 % in SNAP-IT, BS95, and GHQ, respectively, after 12 h of integration.The superiority of the SNAP-IT method is more obvious in the evolution of M 2 , with a final error of 0.8 %, compared with the 2.0 % error in either BS95 or GHQ.Note that the total errors are relatively small because the accumulation mode varies rather slowly.Another simulation with nucleation mode only (i.e., intramodal coagulation) shows that the errors in GHQ and BS95 become three times larger than those in SNAP-IT (figures omitted).For SNAP-IT, GHQ, and BS95, the errors in M 0 are 0.028, 0.092 and 0.090 %, respectively, whereas for M 2 the errors are −0.03,0.103 and 0.103 %.

Verification with regional models
More laborious verifications of the SNAP method are performed here using regional models.We first incorporate the SNAP scheme into a regional atmospheric dust model of Chen et al. (2004), which originally applied 12 size bins for mineral dust.The modified dust model applies two modes of mineral dust particles.The physical processes relevant to dust are emission, transport, gravitational sedimentation and surface depositions, and for the latter two we applied the SNAP scheme.We demonstrate the performance of the SNAP scheme by simulating an East Asian dust storm event that occurred on 19 May 2005, and comparing the simulation with the binned approach.Figure 12 shows the near-surface concentration of number, surface area, and mass of the dust particles.The differences between the binned and SNAP calculations are rather small, with domain average error of 0.65, 1.74, and 8.40 % in M 0 , M 2 , and M 3 , respectively.For this         regional model simulation, the SNAP scheme requires significantly shorter computation time, about 1/3 less including all other overheads, to produce a very similar result to the binned calculation.Most of the time saving is due to the reduced computation time required for particle advection because the SNAP scheme uses 6 variables (3 moments for each mode) to describe the size distributions, as compared with the 12 variables (bins) used for the binned scheme.Figure 12 also shows an additional simulation using the sectional method but with only six bins.The computation cost for this simulation is similar to that for the modal approach because they track the same number of variables.But the six-bin calculations produced significantly larger errors, with domain average error of 31.9, 22.9, and 9.01 % in M 0 , M 2 , and M 3 , respectively.We further examine the size distributions at a location near the dust source (110 • E, 40 • N) and a downstream location between Korea and Japan (130 • E, 35 • N).As shown in Fig. 13, the SNAP size distributions are generally in good agreement with the 12-bin results.The 6-bin distribution looks similar, but its deviation from the 12-bin results is more significant, especially at the downstream location.
A second test was conducted using the CMAQ model, in which we incorporated the SNAP scheme only for the Brownian coagulation process.Three levels of nesting with 81, 27, and 9 km resolutions are applied to simulate particulate pollution over the Taiwan area during early December 2007.The simulation was conducted for eight days including spinup time, and only the last five days' results of the innermost domain were analyzed.However, verification is difficult, as there is no high-resolution binned scheme in CMAQ for verification.Nevertheless, from the analyses shown earlier in Sects.3 and 4, we know that the GHQ method is fairly accurate, so it was used as a benchmark for this comparison.Note that the modal aerosol module in CMAQ does not consider the Kelvin effect, so we also ignored it in the following simulations.Figure 14 shows the 5-day average aerosol dry mass loading simulated with SNAP, and the percent difference comparison against the GHQ method.The two schemes produced similar results.The differences are mostly less than 1 %, and reached 3 % in limited areas.This suggests that the SNAP scheme's performance is close to that of the GHQ scheme in CMAQ.The large relative error occurred mainly over areas where it is raining and the aerosol concentration is low.This also means that the absolute errors at these locations are actually rather small.
An additional test was conducted for the same case to demonstrate the Kelvin effect on aerosol processes.As discussed in Sect.3.4, the Kelvin effect reduces the water content and thus the wet size of hygroscopic aerosol particles, and this effect influences essentially all aerosol processes.Therefore, this simulation included the diagnostic formula for the equilibrium wet size, with the Kelvin effect taken into consideration (see Table A1).Figure 15 shows that when the Kelvin effect is included, aerosol number concentration varies by less than 2 %.However, changes in the higher moments are significant, with a reduction of over 30 % in the cross-section area (M 2 ) and total volume (M 3 ).Most of the changes in M 2 and M 3 were simply due to differences in water content, but the dry aerosol mass loading also changed significantly, with up to a 10 % increase or decrease at various locations.Mechanisms that may contribute to the decrease in dry aerosol volume include less solute uptake as a result of less water content, and enhanced Brownian diffusional deposition due to reduced particle size.A mechanism that may increase dry aerosol volume is reduced gravitational sedimentation, especially for large particles at high humidity.There are certainly many details worthy of discussion that are beyond the scope of this study.The purpose of the simulations here is simply to demonstrate the importance of including the Kelvin effect in the parameterization of aerosol wet size.

Conclusion
An innovative three-moment modal parameterization scheme was developed for accurate simulation of aerosol microphysical processes.Numerical calculations for the growth of a population of aerosol particles, represented by lognormal size distributions, were first performed, and then the results were analyzed by statistical fitting to generate parameterization formulas.Three different approaches were devised for this statistical-numerical aerosol parameterization, namely the kernel transformation (SNAP-KT), integral transformation (SNAP-IT), and optimal-size approximation (SNAP-OS).Another simpler method, the mean-size approximation (MSA), was taken as a no-skill reference.Each SNAP approach might be optimal for a certain process; however, we found that the integral transformation approach is suitable for most of the processes, whereas the optimal-size approximation can occasionally be applied to provide somewhat better parameterizations than SNAP-IT.Although SNAP-KT is outperformed by the other two methods, it is still very useful in obtaining parameterization for diagnostic formulas.These approaches provide parameterization formulas with-out simplifying the growth kernels, and only a minor inaccuracy resulted from the statistical fitting.Rate processes being parameterized include aerosol condensation, Brownian coagulation, sulfuric acid water binary nucleation, and dry deposition.Special attention was given to processes related to aerosol-cloud interactions, and we provided formulas for heterogeneous ice nucleation and wet scavenging, as well as a diagnostic formula for aerosol activation into cloud drops.Other diagnostic formulas provided in this work include considerations for aerosol equilibrium wet size and the Kelvin effect, as well as considerations for the group extinction and absorption coefficients.
The SNAP schemes were verified in various ways, including comparison against numerical solutions, analytical solutions, and results from a binned aerosol parcel model.All comparisons show that SNAP scheme is more accurate than the modal scheme used in CMAQ and WRF-Chem models, including the option that solves the growth integrals with a fifth-order Gauss-Hermite numerical quadrature technique.The computational efficiency of the SNAP scheme is slightly lower (10 to 20 %) than that of the fast scheme in CMAQ, which utilizes lookup tables to speed up calculation; however, it is about 15 times faster than CMAQ's numerical quadrature option.
The SNAP scheme has been implemented in an atmospheric dust regional model, and the results (including the total moments and the dust size distribution) are very close to those simulated using a binned scheme.With such modal parameterization, much computation time is saved, mainly because of the reduced number of variables that need to be considered in advection calculation.We also utilized the CMAQ model to test the integrity of the SNAP scheme, with focus on the Brownian coagulation process.The results indicate that our scheme is as reliable as the fifth-order Gauss-Hermite numerical quadrature scheme.In this model, we further applied a SNAP diagnostic formula for the commonly ignored or simplified Kelvin effect, and showed that this effect cannot be ignored in aerosol modeling.
The parameterization scheme we developed is based on lognormal size distribution.However, detailed bin model simulations indicate that the size distribution may deviate from the lognormal form.It might be worthwhile to revise the scheme based on the gamma-type function, which is suitable for describing skewed size distributions.Because it has no restriction to the number of moments used, the SNAP method can even be applied to the modified gamma distribution, which requires four moments to solve.The SNAP method also has the potential to be used for the modal parameterization of cloud microphysical processes and even other types of physical or chemical processes.Note: g 1,j and g 2,j are SNAP-IT and SNAP-OS adjustment factors for the j th moment (see Sect. 2.2); all µ and r are in m.When combined with MSA to get the full prognostic equations (i.e., Ĩ ), their R 2 values are usually higher than those shown in the last column.1: Ĩ should be calculated with r A = µ • σ 2 and r B = µ/σ 2 (see Sect. 3.3).2: X is sulfuric acid mass fraction of the critical embryo, S W is relative humidity, and S A is relative acidity.
3: Applicable at S W <100 %: r d is dry radius; µ d is dry modal value; κ = j κ j V j j V j , where κ j ≡ i j ρ j Mw ρw M j ; j is species index; V is volume fraction; i is van't Hoff factor; ρ is bulk density; M is molecular weight; and κ j = 0 for insoluble species.

Figure 1 :
Figure 1: Ratios of different mean-size approximation to the true moment as a function of the size distribution spread () for various i values in Eq. (10).

Fig. 1 .
Fig. 1.Ratios of different mean-size approximation Ĩi to the true moment I i as a function of the size distribution width (σ ) for various i values in Eq. (10).

Fig. 2 .
Fig. 2. Dependence of computation precision on bin resolution (number of bins per decade change in size) for intramodal coagulation rates.The errors are calculated by comparing with the 100bins-per-decade results.

Figure 3 :
Figure 3: Fitting surface for the correction factors for the immersion freezing process.Left:

Fig. 3 .
Fig. 3. Fitting surface for the correction factors for the immersion freezing process.Left: g 1 for SNAP-IT; right: g 2 for SNAP-OS.The dots are the original values, and the vertical bars indicate their deviation from the fitting surface.The degree of deviation is also indicated by the color of the dots: blue, green, and yellow represent less than 1, 2, and 3 standard errors, respectively, whereas red denotes greater than 3 standard error.The standard errors of fitting for the left and right panels are 0.83 and 0.04, respectively.

Figure 4 :
Figure 4: Comparing parameterized immersion freezing rates (ordinate) against the numerical solutions (abscissa).Panels from left to right are rates of the zeroth, second, and third moments, respectively.Results from MSA, SNAP-KT, SNAP-IT and SNAP-OS are represented by the blue circles, purple crosses, red dots, and green triangles, respectively.At the lower right corner of each panel is a zoom up of the central section.In the left panel, MSA points with the largest  are highlighted with cyan dots.

Fig. 4 .
Fig. 4. Comparing parameterized immersion freezing rates (ordinate) against the numerical solutions (abscissa).Panels from left to right are rates of the zeroth, second, and third moments, respectively.Results from MSA, SNAP-KT, SNAP-IT, and SNAP-OS are represented by the blue circles, purple crosses, red dots and green triangles, respectively.At the lower right corner of each panel is a close-up of the central section.In the left panel, MSA points with the largest σ are highlighted with cyan dots.

Figure 5 :
Figure 5: (a) Cunningham slip-flow correction as a function of the Knudsen number K N (left ordinate).The exact solution, BS95 and SNAP-KT results are given as the black solid line, blue

Fig. 5 .
Fig. 5. (a) Cunningham slip-flow correction (left ordinate) as a function of the Knudsen number K N .The exact solution, BS95 and SNAP-KT results are given as the black solid line, blue squares, and red dots, respectively.Also shown on the right ordinate are the ratios of BS95 (blue curve) and SNAP-KT (red curve) results to the exact solution.(b) Comparison of parameterized group sedimentation velocity (ordinate) against the exact numerical solution (abscissa).

6Figure 6 :Fig. 6 .
Figure 6: Comparison of parameterized diffusion growth rates (ordinate) with numerical solutions (abscissa): Left: second moment growth rate I 2 (unit: m 2 /particle/s); Right: third moment growth rate I 3 (unit: m 3 /particle/s).MSA is shown in blue open circle, SNAP-IT is in red dot, and BS95 is in green triangle.Also shown are the numerical solutions using from Fukuta and Walter [1970](grey square; labeled as FW).At the lower right corner of each panel is a zoom up of the central section.All rates have been normalized by total number concentration.

Figure 7 :Fig. 7 .
Figure 7: Comparison between various intra-modal Brownian coagulation rates from MSA (blue open circles), SNAP-IT (red dots), Gauss-Hermit quadrature (GHQ; green triangles), and BS95 (purple crosses).Left: rates for (unit: 1/s); Right: rates for (unit: m 2 /particle/s).At the lower right corner of each panel is a zoom up of the central section.All rates have been normalized by total number concentration.
ambient relative humidity and equilibrium size for an ammonium e with 0.01 m dry radius.The red dashed curve is the Köhler curve which includes in effect and Raout effect, whereas the blue curve considers only the Raoult effect.ed line indicates 100% relative humidity.

10Figure 10 :Fig. 10 .
Figure 10: Simulation of the evolution of size distribution due to Brownian coagulation using SNAP-IT (red dashed curve), BS95 (blue dotted curve), GHQ (green dash-dotted curve) and binned model (thick black curve).Thin solid curves indicate the initial size distribution.The left and right pannels are 1 hr and 6 hr results, respectively.Panels from top down are the number, surface area and volume density distributions.

Figure 12 :
Figure 12: Simulated near-surface mineral dust concentrations using the 12-bin sectional scheme

Fig. 12 .
Fig.12.Simulated near-surface mineral dust concentrations using the 12-bin sectional scheme (left column), and the difference (in %) from it using the SNAP scheme (middle) and the 6-bin sectional scheme (right).Domain mean errors are given in the parentheses on the lower right corner of each panel.From top down are the number concentration (M 0 ), surface area concentration (M 2 ), and mass concentration (M 3 ).Only the areas with significant dust concentrations (mass >20 µg m −3 ) are analyzed.

Figure 13 :
Figure 13: Dust particle size distribution calculated with SNAP (red curve), 12-bin scheme (filled triangles) and 6-bin scheme (empty triangles) at two selected locations (geographical coordinates given in the lower left corner) in Fig. 12.

Fig. 13 .
Fig. 13.Dust particle size distribution calculated with the SNAP (red curve) 12-(filled triangles) and 6-bin schemes (empty triangles) at two selected locations (geographical coordinates given in the lower left corner) in Fig. 12.

Figure 14 :
Figure 14: Left: simulated aerosol mass loading over the Taiwan area using SNAP for Brownian

Figure 15 :
Figure 15: Changes in aerosol moments due to the inclusion of Kelvin effect.Panels from left to

Fig. 14 .
Fig. 14.Left: simulated aerosol mass loading over the Taiwan area using SNAP for Brownian coagulation in the CMAQ model.Other panels from left to right: percentage difference between SNAP and GHQ in M 0 , M 2 and M 3 .

Figure 15 :
Figure 15: Changes in aerosol moments due to the inclusion of Kelvin effect.Panels from left to

Fig. 15 .
Fig. 15.Changes in aerosol moments due to the inclusion of the Kelvin effect.Panels from left to right are percent change in number ( M 0 ), surface area of wet particles ( M 2 wet), volume of wet particles ( M 3 wet), and volume of dry particles ( M 3 dry).

Table A2 .
Coefficients for SNAP formulas.