Revisiting Twomey’s approximation for peak supersaturation

Twomey’s seminal 1959 paper provided lower and upper bound approximations to the estimation of peak supersaturation within an updraft and thus provides the first closed expression for the number of nucleated cloud droplets. The form of this approximation is simple, but provides a surprisingly good estimate and has subsequently been employed in more sophisticated treatments of nucleation parametrization. In the current paper, we revisit the lower bound approximation of Twomey and make a small adjustment that can be used to obtain a more accurate calculation of peak supersaturation under all potential aerosol loadings and thermodynamic conditions. In order to make full use of this improved approximation, the underlying integro-differential equation for supersaturation evolution and the condition for calculating peak supersaturation are examined. A simple rearrangement of the algebra allows for an expression to be written down that can then be solved with a single lookup table with only one independent variable for an underlying lognormal aerosol population. While multimodal aerosol with N different dispersion characteristics requires 2N+1 inputs to calculate the activation fraction, only N of these one-dimensional lookup tables are needed. No additional information is required in the lookup table to deal with additional chemical, physical or thermodynamic properties. The resulting implementation provides a relatively simple, yet computationally cheap, physically based parametrization of droplet nucleation for use in climate and Numerical Weather Prediction models.


Background
Atmospheric aerosols are fundamental to the formation of clouds.They provide the substrate onto which liquid droplets can form, overcoming the energy barriers associated with clustering of water molecules (see Mason, 1971;Pruppacher and Klett, 1997, hereafter PK97).Not only do they facilitate the formation of cloud, but they also determine the initial sizes and numbers of droplets as the cloud forms.Variations in the underlying aerosol population; through numbers, sizes, dispersion characteristics and chemical composition; combined with variations in the dynamic and thermodynamic properties of the ambient air; lead to wide ranging differences in the number concentrations of droplets within clouds.Differences in number further impact the physical evolution of the clouds, their lifetime and their interaction with radiation.These impacts are most famously seen through the Twomey (1977) and Albrecht (1989) effects, but also through other mechanisms once feedbacks with the dynamics (e.g.Ackerman et al., 2004;Jiang et al., 2006;Wang and Feingold, 2009) or the ice phase (Lohmann and Feichter, 2005) are taken into account.Collectively, the interactions between aerosol and clouds represent a significant source of uncertainty in our ability to model the climate system and, moreover, to estimate the impacts of anthropogenic climate change (IPCC, Summary for Policymakers, 2007).
The approach generally used to ascertain how many of the total population of the underlying interstitial aerosol will grow into droplets is based upon solving a system for the time variation of supersaturation.Supersaturation is key, since a given aerosol particle will activate to form a droplet at a particular critical supersaturation, s cr , and so the evolution of supersaturation determines which and how many aerosol can activate.For a volume of air undergoing a cooling, the cooling results in an increase in the saturation of the volume.Conversely any condensation that takes place as a result of the increased saturation will itself lead to a decrease in the saturation.An equation that describes the contention Published by Copernicus Publications on behalf of the European Geosciences Union.dσ. (1) Here, the cooling term is represented by α and for an updraft undergoing constant vertical velocity, w, can be written as α = ψ 1 (T , p)w.The condensation term is derived from the equations for diffusional growth kinetics of the droplets, given the differential activity spectrum, ϕ(σ ), which expresses the number of cloud nuclei activated in the interval [σ, σ + δσ ].The terms ψ 1 and ψ 2 along with other notation used in this paper are defined in Table A1.We note that in order to obtain Eq. ( 1) it is assumed that (i) the initial particle size is negligible compared to the droplet size, (ii) particles grow in equilibrium with the environmental relative humidity, and (iii) kinetic and inertial limitations to droplet growth are ignored.These and further approximations are discussed in Ghan et al. (2011), while more detail surrounding the derivation and approximations used to obtain Eq. ( 1) is provided in Appendix A. We also note that although kinetic effects are not explicitly treated in Eq. ( 1), the parametrization of Fountoukis and Nenes (2005) is employed, which uses a size averaged diffusivity coefficient to account for changes in diffusivity with droplet size (cf.Shipway and Abel, 2010).
Equation ( 1) is the bedrock of Twomey's analysis which first enabled the estimation of the fraction of aerosol that would grow into cloud droplets.Once the peak supersaturation has been achieved, i.e. condensation terms start to dominate over the adiabatic source term, and assuming droplets remain in equilibrium with their environment, no more aerosol will be activated.By finding the peak supersaturation, s max , i.e. when ds dt = 0 and the number of nucleated droplets, N c , can simply be determined by the differential activity spectrum as Equations ( 1) and (2), or equations that closely resemble them, have subsequently been used in parametrization schemes which aim to either improve or extend Twomey's approximation (e.g. Feingold and Heymsfield, 1992;Cohard et al., 1998;Nenes and Seinfeld, 2003;Shipway and Abel, 2010;Morales Betancourt and Nenes, 2014).These parametrizations use various methods to tackle the problem of evaluating the integrals contained within the condensation term, but many (e.g.Cohard et al., 1998;Shipway and Abel, 2010;Morales Betancourt and Nenes, 2014) exploit the same lower bound approximation that Twomey developed for evaluating the integral of supersaturation, i.e.
Twomey obtained this bound using a simple geometric argument that considers the areas of triangles bounded above by a line whose slope is given by ds dt = α.Since ds dt is a monotonically decreasing function, the area represented is guaranteed to provide a lower bound approximation to the integral under the supersaturation curve (see the red hatching in Fig. 1).It is somewhat surprising though that this lower bound has also proven to be a relatively good approximation, and what is more it makes the solution of Eq. ( 1) far more tractable.
The following sections revisit this lower bound approximation, develop a slight improvement and subsequently employ this revised approximation in a new parametrization scheme for inclusion in cloud resolving models and General Circulation Models.Section 2 develops an equation set, based on Eq. ( 1), which enables a rapid solution to find the peak supersaturation.The method of solution further allows for a more general approximation to that described by Eq. ( 4) and so Sect. 3 investigates the possibility of improving on Twomey's lower bound.Section 4 then provides some sample results from testing the newly developed parametrization.

A reduced lookup table approach
The peak supersaturation, s max , satisfies the condition ds dt s max = 0 (5) and so Eq. ( 1) with the left hand side set to zero provides an equation for determining peak supersaturation, which can in turn be used to evaluate cloud droplet activation.Retaining Twomey's approximation in Eq. ( 4) and rearranging Eq. ( 1) we find Twomey used a differential activity spectrum of the form ϕ(s) = kCs k−1 .There are limitations to this expression as discussed by Ghan et al. (1993), but moreover, it is desirable to link the activity spectrum directly to the underlying dry aerosol characteristics without the need for a priori information to determine the parameters C and k.To this end (cf.Shipway and Abel, 2010), we use the expression of Khvorostyanov and Curry (2006)  function: This represents a superposition of N modes relating to the modes of an underlying aerosol population.It is common and convenient to represent aerosol in a modal form in aerosol-chemistry models (e.g.Mann et al., 2010), and although the following methodology does not preclude other representations, we consider only lognormal modal representations here.Khvorostyanov and Curry (2006) relate the dry aerosol properties to the critical saturation properties via the relations and Here A k is the Kelvin curvature parameter.The variables r d0 and σ d represent the mean geometric radius and dispersion of the dry aerosol size distribution and are represented by a subscript i for each lognormal mode in Eq. ( 7).Likewise, N a,i are the number concentrations for each mode of the dry aerosol.The relations in Eqs. ( 8) and ( 9) link the underlying dry aerosol physicochemical properties, via b and β, to the activity spectrum (see Khvorostyanov and Curry, 2006;Shipway and Abel, 2010, for more comprehensive details).In essence, b represents the solubility of the soluble fraction, while the representation of the distribution of this soluble fraction within the particle volume is determined by β.In the application demonstrated here, we assume the soluble fraction of the aerosol is uniformly distributed within the particle volume and so β = 0.5.Under this assumption, the parameter b is equivalent to the more commonly used hygroscopicity parameter κ (Petters and Kreidenweis, 2007).
The resulting expression when using this information in Eq. ( 1) becomes The first thing to note about this expression is that all the thermodynamic and dynamic information, i.e. temperature, pressure and vertical velocity, is held on the left hand side of this equation.Thus for fixed aerosol characteristics, the right hand side need only be calculated once to be used in a variety of thermodynamic states.However, if aerosol characteristics also vary in space and time (as is usually the case in a GCM) we cannot exploit this feature of the equation without precalculating all possible aerosol states.
Nevertheless, the computationally challenging part of evaluating the right hand side comes from the integral terms.In this form, in addition to s max , there is a dependence within the integral on aerosol characteristics through the critical supersaturations, s 0,i , and the standard deviations, σ s,i , but not number concentration.This leaves us with three parameters for each aerosol mode; however, a further reduction in the dimension of the problem can be obtained by making a substitution of the form σ = s 0,i σ for each integral within each summand, which then reduces the problem to a form where and where σ , x and y in Eq. ( 13) are dummy variables.Now if we wish to evaluate Eq. ( 11) for a range of thermodynamic, dynamic and aerosol conditions, we need only calculate A and the a i and make use of a precalculated lookup table for I (x, y).For single mode aerosol, the right hand side of Eq. ( 11) reduces to a single term and since A, a i , s 0,i and ln(σ s,i ) are all known, the peak supersaturation can be easily recovered from using an inverse lookup procedure.For a multimodal distribution, Eq. ( 11) must be solved using an iterative approach.Although such an iterative approach makes the method more expensive, the cost of evaluating the right hand side of Eq. ( 11) is now very low and, since I (x, y) is monotonic and well behaved, convergence is rapid and the method remains computationally inexpensive.In the application used throughout the remainder of this paper, a second order Householder's method (i.e.Halley's method) is applied and produces rapid convergence to a high degree of accuracy.

Revisiting Twomey's approximation
The method used in Sect. 2 retained the original approximation of Twomey as given by Eq. ( 4).Unlike approaches such as those adopted by Shipway and Abel (2010) or Morales Betancourt and Nenes (2014), this method no longer seeks to evaluate the inner integral with a closed analytical form, but rather precalculates the integral numerically.This gives us the freedom to choose an expression for the approximation to the supersaturation integral which might make analytical progress more difficult, but provides a more accurate approximation.
Figure 1 illustrates the approximation (4), which can be interpreted as a right trapezoid whose uprights have length s and σ with the top most edge having a gradient α (red hatched region).The width of the trapezoid is thus determined by this gradient.If the local gradient at s deviates significantly from α, i.e. as s approaches s max , the width of the trapezoid will be significantly less than (t − τ ) and the approximation will significantly underestimate the integral.
However, if we maintain this trapezoidal representation, then it can be clearly seen that taking the gradient of the upper edge to be equal to the mean gradient of the supersaturation curve, i.e. (s−σ )  (t−τ ) , will lead to a much more accurate approximation (black hatched region).The problem now becomes one of determining an expression for the gradient of the supersaturation curve that we can make use of in our evaluation of Eq. (1).

A parametrization for ds dt
A feature of the evolution of the supersaturation gradient is that we know that for small values of time, ds dt ≈ α, while when s = s max , ds dt = 0. Thus we might expect a scaled form for ds dt to be given by 1 α It is noted that this functional form, i.e. with argument s s max , can be very easily used in the analysis of Sect.2, since the required scaling by s 0,i within each integral will return a function in terms of x = s max s 0,i .To investigate the behaviour of ds dt , we perform a numerical integration of the full expression for Eq. ( 1) given the activity spectrum in Eq. ( 7) using a range of aerosol loadings and updrafts.A composite of the data from these numerical integrations is shown by the grey shaded area in Fig. 2, where the data are normalized by the anticipated scalings.It is striking from these data how little spread there is, suggesting that a single fitted curve may suffice to represent all regimes.
We choose to fit a curve through these data of the form where a least squares optimization on the envelope of data leads to the choice that µ ≈ 3 and λ ≈ 0.6.These values are used in the final parametrization.The functional form chosen in 16 is motivated solely by inspection of the data in Fig. 2; however, further discussion of this approximation, and the evaluation of µ and λ is given in Appendix B.

A revised approximation for the inner integral of s
Having facilitated a suitable expression for ds dt , we can now exploit this to obtain a revised form for Twomey's approximation in Eq. (4).For this, we simply take the gradient of the topmost edge of the trapezoid in Fig. 1 (which is assumed to be α in Twomey's approximation) to be the mean gradient given by ds dt or when s = s max , simply ds dt Our revised approximation for the inner integral of s, becomes In order to understand the impact of this revised approximation, the relation ( 19) is applied to the full form of Eq. ( 1) and a numerical integration is performed.This has been carried out for a range of scenarios (see Sect. 4), but an example evolution of supersaturation is provided in Fig. 3.This figure shows three evolution curves when using the Whitby marine aerosol loadings (Whitby, 1978) under conditions with T = 279 K, p = 1000 mb and an updraft velocity of w = 0.5 m s −1 .The solid line represents the evolution when the full equation set is evaluated numerically with no approximations to the inner integral.The dashed line shows the effect of using Twomey's original approximation, Eq. ( 4), and demonstrates that using a lower bound approximation to the inner integral necessarily underestimates the contribution from the condensation term and thus overestimates the peak supersaturation (NB when we employ the approximation terms we only do so up until peak supersaturation is achieved and thus the downward branch of the supersaturation curve is not simulated).The dotted curve shows how the supersaturation evolves if we replace Twomey's approximation with that from Eq. ( 19) where the mean gradient of ds dt is calculated using the parametrization (Eq.16).A similar integration where the gradient is calculated explicitly as the curve evolves results in a curve that almost exactly overlies the dotted curve and is thus not shown in Fig. 3.These curves demonstrate the benefits of using the revised parametrization in Eq. ( 19), but further suggest that using the parametrization in Eq. ( 16) will suffice for an approximation to the gradient when applied to the full parametrization.For this particular example, the revised approximation leads to a peak supersaturation of 0.518 % compared with the value from the integration with no approximation of 0.508 %, i.e. a 2 % overes-  1).Broken lines show numerical integration using various approximations to the inner integral.

A revised parametrization to determine s max
Section 2 provides a new parametrization for the determination of peak supersaturation and thus the number of activated droplets.By itself, it represents a fast and accurate method, which has an accuracy comparable with many existing parametrizations (see Sect. 4).By exploiting the revised approximation in Eq. ( 20) we can potentially improve the accuracy further with minimal modification.In fact, the only change that is necessary, is to update the expression for the integral (Eq.13) which is used in the lookup tables and is independent from any properties of the aerosol or thermodynamics.The final expression for I (x, y) if we make use of Eq. ( 20) and rescale according to σ = s 0,i σ , simply becomes where µ and λ are taken to be 3 and 0.6, respectively.Precalculation of I (x, y) allows a rapid solution to Eq. ( 11), where y may be taken as a continuous variable or, if the modal representation of the aerosol has fixed dispersion characteristics, can simply use appropriate discrete values. .Evaluation of peak supersaturation for the parametrization with the original (Eq.4) and revised (Eq.20) approximations to the integral under the supersaturation curve.These are compared to the full numerical solution of Eq. ( 1).Results from other well-used parametrization schemes are also shown, together with results from the parcel model used in Ghan et al. (2011).Each panel provides results from different aerosol scenarios as described by Whitby (1978) and range from relatively clean marine conditions to highly polluted conditions.

Sample results
The recent paper of Ghan et al. (2011, hereafter G11) provides an overview of droplet nucleation parametrizations and their relative performance across a range of aerosol loadings and dynamical situations.We refrain from providing an exhaustive evaluation in the current expository text, but note that in the majority of scenarios tested within the Ghan paper, the parametrization presented here performs as well, if not better than the Shipway and Abel ( 2010) parametrization (labelled "Shipway" in the plots of G11).The Shipway and Abel (2010) parametrization uses a similar approach to that used in the current work when considering the representation of the differential activity spectrum, while it also uses the Twomey approximation to determine peak supersaturation.
However, as a sample demonstration, plots are presented in Figs. 4 and 5, which use the Whitby aerosol loadings (Whitby, 1978) over a range of updraft velocities.These loadings use three modes of aerosol (Aitken, accumulation and coarse) and represent samples from relatively clean marine conditions to very polluted urban conditions.These figures show a comparison of various parametrizations and nu-merical solutions for the evaluation of s max and the fraction of activated aerosol, respectively.The red and green curves respectively represent the solution obtained with the look-up table procedure described here using Twomey's approximation and the newly developed approximation in Eq. ( 20).The dashed curves represent the peak supersaturation obtained when a full integration of Eq. ( 1) is performed without any further approximation.The remaining curves are taken from the paper of G11 for comparison and represent parametrizations of Abdul-Razzak and Ghan (2000, "ARG") and Fountoukis and Nenes (2005, "Nenes"); two parametrizations that are widely used in climate models.The diamonds show the solution obtained from the numerical parcel model used in G11, which represents the complete evolution of droplet growth without the approximations implicit in the derivation of Eq. (1) (see Appendix A).
The first thing to note from these plots is the discrepancy between the full numerical implementation and the results obtained from the numerical integration of Eq. ( 1).For marine aerosol loadings, these two curves are very closely aligned across all updraft values, suggesting that in this regime with lower aerosol number concentrations, the errors associated with the approximations described in Appendix A are unimportant.However as aerosol numbers are increased, there may be an increased competition between the neglect of the influence of the initial droplet radius and the y term (again see Appendix A): for the Urban aerosol, with very high number concentrations, the coarse mode aerosol dominate the nucleation process, compared to the Clean Continental and Background scenarios in which more accumulation mode aerosol are activated (see .This suggests that the critical radius for activation is larger in the Urban scenario and that the influence of the initial droplet radius is greater.However, the larger number of activated aerosol and the low peak supersaturations in the Continental and Background scenarios compared to the Marine scenario, might suggest that the y term becomes more influential (see Appendix A).This speculation doesn't take into account the way kinetic effects are treated.A more extensive investigation into the attribution of the differences between results from Eq. (1) and a full numerical treatment would shed light onto areas for future improvement of parametrization schemes, but this is beyond the scope of the current investigation.Despite these differences, it is apparent from the data in Fig. 4, that Twomey's Eq. ( 1) provides a very good estimate of peak supersaturation.If we now turn our attention to the red and green curves representing the parametrizations developed in this work, we immediately see the benefit of the re-vised approximation to Twomey's original.Across all scenarios, the revised approximation (green curve) sits very close to the dashed line of Eq. (1), while Twomey's lower bound approximation (red curve), as one would expect, consistently overestimates the peak supersaturation.That said, this latter overestimation in certain circumstances can be seen to compensate for some of the differences between the parcel model and Eq.(1) discussed above.
A more comprehensive analysis of the performance of this and other parametrizations has been conducted by Partridge (2015) and will be reported elsewhere in the literature.

Summary
A new method for parametrizing peak supersaturation, s max , and thus droplet nucleation, has been developed.This method involves reducing the number of degrees of freedom within the integral expressions involved in the equation for s max to 2 in the case of N aerosol modes with variable dispersion characteristics (as opposed to 2N + 1 before the transformation) and only 1 for fixed dispersion.This reduction makes it feasible to precalculate the integral expressions using numerical methods and then reuse these data in the form of lookup tables when solving for s max .
This method in itself represents a much more computationally efficient method for determining aerosol activation and links through to the underlying aerosol physicochemistry.However, another benefit is that it is not constrained to use the fixed form of the lower bound approximation of Twomey (1959), which is frequently used to make analytic integration feasible, and so a more accurate approximation is derived.This newly derived approximation perhaps lacks the elegance of Twomey's original estimate, but is shown to faithfully reproduce the evolution of supersaturation and the calculation of the peak supersaturation across a range of scenarios.
The underlying equation set on which the parametrization is built makes a number of physical assumptions, which are common to many well used activation parametrizations (Ghan et al., 2011).These assumptions and the resulting consequences are discussed, but the accuracy and computational efficiency of the resulting parametrization make it suitable for inclusion in both GCMs and high resolution models that use a modal representation of aerosol.

Appendix A: Derivation of Twomey's supersaturation equation
A complete derivation of Eq. ( 1) is provided in chapter 13 of PK97, but we briefly summarize some of the approximations that are necessary to obtain this form.The starting point is with a more general expression for the condensation term (cf.Eqs.13-29 of PK97), namely where w L is the total condensed water mixing ratio.Given the differential activity spectrum, ϕ(σ ), the rate of change in the total condensate mass can be expressed as where dm dt represents the rate of growth of condensate mass on each droplet of a given size as it becomes active.For spherical drops, m = 4 3 π ρ w a 3 , where a is the drop radius.The rate of growth of a droplet by diffusional growth can be approximated by The expression for y is given in Table A1 and represents the Kelvin effect, which describes the way saturation vapour pressure changes due to the curved surface of a droplet , and Raoult's law, which accounts for the change in vapour pressure due to the solute.We note that for large values of the droplet radius, both terms in y will become small.Thus the first approximation we will make is that y can be neglected in Eq. (A3).We also note that the diffusivity and thermal conductivity terms within the expression for G(T , p) are actually dependent on the size of the droplet.This dependence is also neglected in the following derivation, although in the application of the parametrization the approach of Fountoukis and Nenes (2005) is used to determine a representative mean diffusivity (cf.Shipway and Abel, 2010).
These approximations remove the dependence of the right hand side of Eq. (A3) on a and so we can integrate to find To complete the derivation of Eq. ( 1) we must make further approximations, but it is worthwhile considering the impact that these may have before we do.The assumption that y s, will serve to overestimate the rate of condensation (similarly when the approximation is applied to derive Eq.A3).This assumption becomes decreasingly valid for small droplet sizes and low saturation conditions.I.e. for a very large number of small droplets that depress the peak supersaturation, this approximation may lead to an overestimate of the condensation rate, and thus an underestimate of the peak supersaturation.Twomey employs the further assumption that the initial droplet size is small and can be neglected.In this case, the validity of the approximation is weakened in conditions of low supersaturations, but with larger droplets dominating the condensation process.Large numbers of large aerosol will likely lead to approximation of this term underestimating the condensation rate and overpredicting the peak supersaturation.Having made these two further approximations, the form for Eq.(A3) becomes By substituting Eq. ( A6) into (A1), and making the final approximation that s 1 within the bracketed term of (A1), we recover the Eq. ( 1).where 2 F 1 (a, b; c; z) is the Gauss hypergeometric function (Abramowitz and Stegun, 1964, chapter 15).This by itself provides us with an approximate evolution of s, but since it is in terms of s max does not provide any insight into the problem at hand of determining s max .However, Pinsky et al. (2013) use a scale analysis to demonstrate that the time at which peak supersaturation is achieved, t max is well approximated by where C 1 = 1.058 and C 2 = 1.904 are universal coefficients that are independent of any physical variables.Taking s = s max and using Eq.(B2) in Eq. (B1) leads us to an implicit relation for µ and λ: 2 F 1 λ, which, given the argument is unity in the hypergeometric function, further reduces to (see Abramowitz and Stegun, 1964, Eq. 15.1.20).Figure B1 shows the solution space for µ and λ along with least squares error estimates of µ given λ using the data in Fig. 2. The point of intersection suggests values of µ = 2.64 and λ = 0.515 would be consistent with both minimizing to the numerical data and the scale analysis of Pinsky et al. (2013).Comparison of the solid and dashed lines in Fig. 2, provides reassurance that the final parametrization employed in Eq. ( 21) is consistent.The result of Pinsky et al. (2013) further provides a basis for testing the validity of the parametrization of the gradient of ds dt between the activation time, τ , and the time of peak supersaturation approximated by the mean gradient in Eq. ( 18).There it is assumed that the monotonic behaviour of the ds dt curve allows for the mean ds dt value to be used as an approximation to (s−σ )  (t−τ ) .If we consider the extreme ends of the curve, i.e. (s max −0) (t max −0) , then the relation (B2) suggests that Thus the mean ds dt value in Eq. ( 18) represents only a slight underestimate of the gradient as depicted in Fig. 1.

Figure 2 .
Figure 2. Scaled dsdt as a function of s s max evaluated from numerical solution of Eq. (1) using a range of updraft velocities and the various Whitby aerosol loadings.The solid line represents the parametrization (Eq.16) with µ = 3 and λ = 0.6.The dashed line uses the coefficients determined in Appendix B.

Figure 3 .
Figure 3. Evolution of supersaturation using Whitby marine aerosol loading.Solid line shows numerical integration of the full equation set in Eq. (1).Broken lines show numerical integration using various approximations to the inner integral.

Figure 5 .
Figure 5.As Fig. 4 but showing the resulting fraction of aerosol which are activated to become droplets.
A4)    where τ is the activation time of a given droplet.Equations (A4) and (A3) together allow us to write the mass change

Figure B1 .
Figure B1.Solution space for the parameters µ and λ as determined from condition (B4) (solid curve) and least squares fit of data in Fig.2to the functional form Eq. (16) (dashed curve).

Table A1 .
List of symbols.