On the analytic approximation of bulk collision rates of non-spherical hydrometeors

Analytic approximations of the binary collision rates of hydrometeors are derived for use in bulk microphysical parameterizations. Special attention is given to nonspherical hydrometeors like raindrops and snowflakes. The terminal fall velocity of these particles cannot be sufficiently well approximated by power-law relations which are used in most microphysical parameterizations, and therefore an improved formulation is needed. The analytic approximations of the bulk collision rates given in this paper are an alternative to look-up tables and can replace the Wisner approximation, which is used in many atmospheric models.


Introduction
The approximation of bulk collision rates is a classic problem in the formulation of cloud microphysical parameterizations for atmospheric models.The most common formulation is the continuous growth equation which applies to particles of very different sizes and fall speeds; i.e., the size and the fall speed of the smaller and thus slower falling particle is neglected (Rogers and Yau, 1996;Pruppacher and Klett, 1997;Straka, 2009).In practice, collisions of different particles of similar size and fall speed do occur, and one needs a formulation which is more general than the simple continuous growth equation.The standard approach for the collision rate of two ensembles of precipitation-sized particles goes back to Wisner et al. (1972), who used the ansatz that the velocity difference can be approximated by the difference of weighted means which then simplifies the solution of the collision integrals.If the weighting functions are properly chosen, this Wisner approximation can recover the continuous growth equations as asymptotic limits.The main disadvantage of this approach is that the bulk collision rate of the Wisner approximation becomes zero, if the difference of the weighted mean sedimentation velocities is zero, while for broad size distributions as they occur in nature the bulk collision rate is always non-zero.Verlinde et al. (1990) studied analytic solution of the bulk collision integrals for generalized gamma distributions and terminal fall velocities which can be approximated by powerlaw relations.They provide exact solutions for this problem, but those include the general hypergeometric function and are therefore difficult and expensive to evaluate.Gaudet andSchmidt (2005, 2007) derived a generalized form of the Wisner approximation, which overcomes some of the deficiencies of the classic approximation.Straka and Gilmore (2006) discussed the effect of non-spherical raindrops on the bulk collision rate and argued that the usual approximation of spherical geometry in combination with a power-law fall speed is sufficiently accurate, but that this is the result of the cancellation of two relatively large errors.Ćurić and Janc (2010) showed that the use of size distributions without an upper cutoff diameter can lead to significant errors, especially for raindrops and hail.
An alternative to the Wisner approximation was suggested by Murakami (1990) and Mizuno (1990) and later by Seifert andBeheng (2006, SB2006 hereafter).The latter used the variance (or standard deviation) to approximate the bulk difference of the sedimentation velocities.Like most earlier Published by Copernicus Publications on behalf of the European Geosciences Union.

A. Seifert et al.: Approximation of bulk collision rates
parameterizations SB2006 limited their formulas to powerlaw relations for the terminal fall velocities.
In this paper we revisit the variance approximation introduced by SB2006.Especially for raindrops and snowflakes the power-law relations for the terminal fall velocity are not valid for large particles.Therefore new equations are derived which yield a better approximation for the collision rates that involve these particles.Additionally, the SB2006 variance ansatz is improved by introducing a more general weighting factor.For raindrops we include the non-spherical geometry of large raindrops for consistency with the more accurate terminal fall velocity approximation.
An alternative to analytic approximations of the integrals is look-up tables as they are used or recommended by many previous studies in their attempt to improve the parameterization of the collision integrals (e.g., Walko et al., 1995;Straka and Gilmore, 2006;Saleeby and Cotton, 2008;Thompson et al., 2008).In our opinion, look-up tables have some disadvantages and their use should be limited to cases where analytic solutions or approximations are not available.The reasons which made us discard the use of look-up tables are as follows: 1. Look-up tables increase the complexity of the modeling system because an automatic pre-processing becomes necessary which guarantees that the look-up tables are consistent with the microphysical parameters chosen for the specific simulation.This can become especially cumbersome in operational numerical weather prediction where reproducibility is of the essence.
2. Large multi-dimensional look-up tables can become inefficient because of the additional memory access, especially when the tables are larger than the cache size.On today's supercomputers atmospheric models are often memory bandwidth limited.Doing more floating point operations without additional memory access can thus be more efficient than using look-up tables.
3. Look-up tables are usually not included in the publications, nor are the corresponding source codes.This makes it difficult to reproduce the results from such models and hinders scientific progress.
4. While analytic approximations allow further theoretical studies, e.g., to explore the sensitivity to certain parameters or assumptions, this becomes limited to numerical studies once look-up tables have been introduced.
Therefore analytic approximations are in our opinion an important part of microphysical parameterizations.Another alternative to look-up tables are fits with, e.g., rational functions as used, for example, by Frick et al. (2013).In the end the question of which implementation strategy is most efficient depends very much on the application and the hardware architecture of the supercomputer and the processor.The paper is organized as follows: in Sect. 2 we review and discuss the geometry and terminal fall velocity of raindrops and snowflakes and their approximation in bulk microphysical parameterizations.Sections 3 and 4 shortly introduce the Wisner and variance approximation of the bulk collision rates.In Sects.5 and 6 we derive explicit parameterizations of the bulk collision rates of the collision between graupel and rain, and between snow and rain.Section 7 deals with the self-collection of snow.In Sect.8 the previously introduced as well as other binary collision interactions are discussed by use of quantitative error measures.Some conclusions are given in Sect.9.The Supplement contains auxiliary figures and further technical details on the collision interactions which are not in discussed in the main text.

Geometry and fall speeds of raindrops and snowflakes
Many hydrometeors in the atmosphere can be approximated as spheres, and this approximation is often made in bulk microphysical models.As is well known, the assumption of sphericity does not hold for raindrops larger than 1 mm diameter because the aerodynamic pressure forces lead to an oblate shape and oscillations (e.g., Beard and Chuang, 1987;Szakall et al., 2010).The axis ratio of raindrops where D r,min and D r,max are the minimum and maximum dimensions of the raindrop, is shown in Fig. 1a based on the theoretical model of Beard and Chuang (1987), the polynomial fit of Chuang and Beard (1990), and the parameterization following Khvorostyanov and Curry (2002) given by The length scale D η is a constant which we chose to be 5.2 mm.Compared to Khvorostyanov and Curry (2002) this is a slightly higher value, but it fits the data of Beard and Chuang (1987) better than the 4.7 mm used by Khvorostyanov and Curry (2002).Assuming a perfect oblate spheroid, the mass of the raindrop is given by x r = (π/6)ρ w ηD 3 r,max , i.e., given the equivalent diameter of the raindrop Eq. ( 2) constitutes an implicit equation for the maximum dimension.
For the parameterization of the bulk collision rate we need a simpler explicit equation for D r,max as a function of the equivalent diameter D r .In the following we use  2) with Dη = 5.2 mm, the results of the Beard and Chuang (1987) model and the polynomial fit of Chuang and Beard (1990).For the maximum diameter as a function of the equivalent diameter (right) the explicit relation Dr,max = Dr exp(ωrDr) with ωr = 33 m −1 is used as an approximation to the implicit relation given by Eq. ( 2).
asymptotics, i.e., a constant value at large diameters.Most of these relations are of the form or some variant thereof.Here D r is again the equivalent diameter of raindrops, and α r , β r and γ r are constant coefficients.In the following we will call Eq. ( 4) an Atlas-type fall speed relation.In contrast to that, bulk microphysical parameterizations traditionally use the less accurate power law approximation for the terminal fall velocity (e.g.Kessler, 1969;Liu and Orville, 1969;Lin et al., 1983, etc.), which increases without bounds for large diameters.
Over the last decades an aerodynamic theory has been developed, which predicts the terminal fall velocity of arbitrarily shaped particles over the whole size range, i.e, a large range of Reynolds numbers.The basic idea to apply boundary layer theory to the problem goes back to Abraham (1970).Later Böhm (1989Böhm ( , 1992) ) showed that the boundary layer theory not only provides a powerful approach to treat the large range of Reynolds numbers, but can also overcome the difficulties of arbitrarily shaped complicated particles.Böhm's approach was then further developed by Mitchell (1996) as well as Khvorostyanov andCurry (2002, 2005).
Figure 2 shows the terminal fall velocity of raindrops using the empirical relation of Beard (1976).The theoretical model of Khvorostyanov and Curry (2005, KC05 hereafter) agrees reasonably well with the empirical data, but underestimates the fall velocity of raindrops between 2 and 5 mm diameter.
The simple power law formula is as quite crude approximation and can only give a rough estimate over a limited size range.Especially for large raindrops the power law overestimates dramatically.In contrast, the Atlas-type fall speed relation provides a good approximation over a large size range.
The only disadvantage is that it underestimates the terminal fall velocity of small drops, i.e., it does not approach the Stokes law and becomes even negative for very small drops.Nevertheless, for water drops larger than 0.1 mm diameter, i.e. for raindrops, it provides a very good parameterization.
Snowflakes show a similar behavior to raindrops in the sense that for large snowflakes the terminal fall velocity becomes approximately independent of size.The reason for this is that the particle bulk density decreases with increasing size of the snowflake, i.e., snowflakes grow faster than constant density spheres.The geometry of snowflakes can be described by a quadratic mass-size-relation (Locatelli and Hobbs, 1974;Brandes et al., 2007) x s = a s D 2 s (5) with the maximum dimension D s , the mass of the snowflake x s and a constant parameter a s .The quadratic mass-sizerelation is typical for growth regimes dominated by aggregation as shown by Westbrook et al. (2004a,b).Such relations Fig. 1.Parameterization relations for non-spherical geometry of raindrops.The axis ratio as a function of equivalent diameter (left) following from Eq. ( 2) with D η = 5.2 mm, the results of the Beard and Chuang (1987) model and the polynomial fit of Chuang and Beard (1990).For the maximum diameter as a function of the equivalent diameter (right) the explicit relation D r,max = D r exp(ω r D r ), with ω r = 33 m −1 , is used as an approximation to the implicit relation given by Eq. ( 2).
with ω r = 33 m −1 , and Fig. 1b suggests that this is a sufficiently accurate approximation.
The fact that raindrops deviate from a spherical shape results in an increased cross-sectional area and a higher drag; thus it decreases the terminal fall velocity compared to a spherical particle of the same mass.Together with the high Reynolds number of a falling raindrop, i.e., the effects of turbulence, this explains the well-known empirical result that the terminal fall velocity of raindrops becomes independent of drop size for large drops, i.e., that the terminal fall velocity approaches a constant value.Many parameterizations of the terminal fall velocity of raindrops take this fact into account, e.g., the formulas by Best (1950), Atlas et al. (1973), Rogers et al. (1993) or Kogan and Shapiro (1996) give the correct asymptotics, i.e., a constant value at large diameters.Most of these relations are of the form or some variant thereof.Here D r is again the equivalent diameter of raindrops, and α r , β r and γ r are constant coefficients.In the following we will call Eq. (4) an Atlas-type fall speed relation.In contrast to that, bulk microphysical parameterizations traditionally use the less accurate power-law approximation for the terminal fall velocity (e.g., Kessler, 1969;Liu and Orville, 1969;Lin et al., 1983), which increases without bounds for large diameters.
Over the last decades an aerodynamic theory has been developed which predicts the terminal fall velocity of arbitrarily shaped particles over the whole size range, i.e, a large range of Reynolds numbers.The basic idea to apply boundary layer theory to the problem goes back to Abraham (1970).Later Böhm (1989Böhm ( , 1992) ) showed that the boundary layer theory not only provides a powerful approach to treat the large range of Reynolds numbers but can also overcome the difficulties of arbitrarily shaped complicated particles.Böhm's approach was then further developed by Mitchell (1996) as well as Khvorostyanov andCurry (2002, 2005).
Figure 2 shows the terminal fall velocity of raindrops using the empirical relation of Beard (1976).The theoretical model of Khvorostyanov and Curry (2005, KC05 hereafter) agrees reasonably well with the empirical data, but underestimates the fall velocity of raindrops between 2 and 5 mm diameter.The simple power-law formula is as quite crude approximation and can only give a rough estimate over a limited size range.Especially for large raindrops the power law overestimates dramatically.In contrast, the Atlas-type fall speed relation provides a good approximation over a large size range.
Here and in the following we have used the parameters The only disadvantage is that they underestimate the terminal fall velocity of small drops; i.e., they do not approach the Stokes law and even become negative for very small drops.Nevertheless, for water drops larger than 0.1 mm diameter, i.e., for raindrops, it provides a very good parameterization.
Snowflakes show a similar behavior to raindrops in the sense that for large snowflakes the terminal fall velocity becomes approximately independent of size.The reason for this is that the particle bulk density decreases with increasing size of the snowflake; i.e., snowflakes grow faster than constant density spheres.The geometry of snowflakes can be described by a quadratic mass-size relation (Locatelli and  2002) can explain the empirical behavior, but significant differences exist (grey dashed line).The power law approximation (blue dotted line) is in general inappropriate, but the Atlas-type relation (red dash-dotted line) gives a good approximation for raindrops larger than 0.1 mm.
Table 1.Coefficients for the mass-size relation x = âD b, the maximum dimension as function of particle mass D = ax b , the power law terminal fall velocity v(x) = αx β of particles with mass x and the shape parameters of the Gamma distributions f (x) = Ax ν exp(−Bx ξ ) and f (Deq) = N0D µ eq exp(−λDeq).Note that for the raindrops only the spherical geometry is given here, but the non-spherical correction is taken into account explicitly by Dr,max = Dr exp(ωrDr).For the area-size relation we give the pre-factor γ A in the formula A = γ A D 2 where A is the cross sectional area and D is the maximum dimension.difference.As mentioned in the introduction, Verlinde et al. (1990) derived a general analytic solution for particle distributions in form of gamma distributions combined with power law relations for the terminal fall velocities.However, this analytic solution leads to the general hypergeometric function which by itself is very difficult to evaluate.In practice, the analytic solution is computationally as expensive as the numerical calculation of the integral itself.Therefore many schemes do still apply the approximation introduced by Wisner et al. (1972) who replaced the actual terminal fall velocity by some bulk mean value v which does not depend on D.
Doing the same also for the collision efficiency yields x n j dD i dD j , and with Fig. 2. Terminal fall velocity of raindrops as a function of equivalent diameter using different approaches.The empirical relation of Beard (black solid line) is regarded as the reference.The aerodynamic theory of Khvorostyanov and Curry (2002) can explain the empirical behavior, but significant differences exist (grey dashed line).The power-law approximation (blue dotted line) is in general inappropriate, but the Atlas-type relation (red dash-dotted line) gives a good approximation for raindrops larger than 0.1 mm.Hobbs, 1974;Brandes et al., 2007), with the maximum dimension D s , the mass of the snowflake x s and a constant parameter a s .The quadratic mass-size relation is typical for growth regimes dominated by aggregation as shown by Westbrook et al. (2004a, b).Such relations are used in many bulk microphysical parameterizations; e.g., typical values for a s are between 0.038 kg m −2 (Doms and Schättler, 2002;Baldauf et al., 2011) and a s = 0.069 kg m −2 (Wilson and Ballard, 1999;Field et al., 2005).Heymsfield (2003, his Eqs. 4 and 5) suggests that the lower value might be more appropriate for the tropics while the higher value is consistent with data from mid-latitudes.
Figure 3 shows the terminal fall velocity of low-density snowflakes with a s = 0.038 kg m −2 .To apply the KC05 aerodynamic model, an additional assumption about the cross-sectional area A s has to be made.Based on the observations of Field et al. (2008), namely the data presented in their Fig. 3, we chose for simplicity a constant area ratio âs = A s /[(π/4)D 2 s ] = 0.45.This should provide a good estimate for the small and medium-sized snowflakes, but might overestimate the cross-sectional area for snowflakes larger than 3 mm diameter and lead to an underestimation of the terminal fall velocity.Using, on the other hand, the area-size relation A s = 0.2285D 1.88 s of Mitchell (1996) for aggregates would result in about 20 % lower terminal fall velocities than applying âs = 0.45.The uncertainty in the terminal fall velocity of large snowflakes is obviously large, even for a given mass-size relationship.
For the Atlas-type approximation of resulting terminal fall velocity of snowflakes, we have the choice between a formulation using equivalent diameter where D eq is the diameter of a liquid sphere, i.e., and a formulation with maximum dimension Note that we use D eq for the equivalent diameter of snow and D r for the equivalent diameter of raindrops.This is only an attempt to make some of the equations more readable, and both equivalent diameters are of course identical.
To optimize the parameters in these fits, we apply the Nelder-Mead downhill simplex method (Press et al., 1992), and with the resulting parameter values we find that the formulation using the equivalent diameter D eq provides the more accurate fit to the KC05 terminal fall velocity.Especially from the log-log plot in Fig. 3b it becomes obvious that any power-law approximation can only provide a good fit for a small size range of snowflakes, but not for the whole relevant sizes range from, at least, 0.1 mm to 10 mm maximum dimension.
For all hydrometeors used in this study, the geometries, terminal fall velocities and some further assumptions are summarized in Table 1.Hail is assumed as spherical particles with a (constant) particle density of 920 kg m −3 .The cloud ice is hexagonal plates using the relations as given by Mitchell (1996), and for graupel we use a mass-size relation for lump graupel of Locatelli and Hobbs (1974).Note that these are not necessarily the particle properties assumed in the SB2006 two-moment scheme (or later publications using that scheme), but instead the particle properties here are chosen to span a wide, but typical, range of parameters.Some more details are discussed in Sects. 2 and 3 of the Supplement.

Wisner approximation
The classic gravitational collection kernel as it is found in most textbooks (e.g., Rogers and Yau, 1996;Pruppacher and Klett, 1997) is given by with particle diameters D i , D j ; terminal fall velocities v i , v j ; and the collision-coalescence or aggregation efficiency E ij .
For oblate spheroids the relevant particle diameters are the maximum dimensions, while for more complicated geometries one would have to use an area-equivalent spherical diameter perpendicular to the fall direction.Assuming, without loss of generality, that i is the collecting species and j the one which is collected, the spectral collection rate for the species j , i.e., the loss term in the budget equation of the particle size distribution f j (D j ), is given by As usual, the particle size distribution f (D) is defined as the number of particles per unit volume in the size range [D, D+ dD].Multiplication with x n j and integration over the internal coordinate D j leads to the bulk collision rates for the integral moments of the collected species x n j dD i dD j ; i.e., M j,1 is the mass density L j of the collected species and M j,0 = N j its number density.The corresponding tendencies of the collecting species are ∂N i ∂t = 0 and The integral in Eq. ( 11) is in general very hard to solve analytically because of the absolute value of the fall speed difference.As mentioned in the Introduction, Verlinde et al. (1990) derived a general analytic solution for particle distributions in the form of gamma distributions combined with power-law relations for the terminal fall velocities.However, this analytic solution leads to the general hypergeometric function, which by itself is very difficult to evaluate.In practice, the analytic solution is computationally as expensive as the numerical calculation of the integral itself.Therefore many schemes do still apply the approximation introduced by Wisner et al. (1972), who replaced the actual terminal fall velocity by some bulk mean value v which does not depend on D. Doing the same also for the collision efficiency yields and with we can write this as With the usual assumptions about the particle size distributions, the remaining integral in C n,ij can be solved quite easily.Wisner et al. (1972) specified both bulk terminal fall velocities as the mass-weighted fall speeds.A more detailed analysis of the asymptotic behavior, which should recover the continuous growth equations, shows that the bulk terminal fall velocity of the collecting particles should be weighted with D 2 while the bulk fall speed of the collected particles has to be weighted with D 2 x (Seifert, 2002).When using a two-moment scheme the equation for the number densities should apply D 2 -weighted fall speeds for both species.The bulk velocities in the Wisner approximation are therefore calculated by and n = 0 applies to number density, and n = 1 to the mass density equation.Note that for the collecting species the bulk velocity does not depend on n; i.e., the same velocity is used for all moments.With the usual assumptions about the particle size distributions, the remaining integral in C n,ij can be solved quite easily.Wisner et al. (1972) specified both bulk terminal fall velocities as the mass-weighted fall speeds.A more detailed analysis of the asymptotic behavior, which should recover the continuous growth equations, shows that the bulk terminal fall velocity of the collecting particles should be weighted with D 2 while the bulk fall speed of the collected particles has to be weighted with D 2 x (Seifert, 2002).When using a two-moment scheme the equation for the number densities should apply D 2 -weighted fall speeds for both species.The bulk velocities in the Wisner approximation are therefore calculated by and n = 0 applies to number density, n = 1 to mass density equation.Note that for the collecting species the bulk velocity does not depend on n, i.e., the same velocity is used for all moments.

Variance approximation
The Wisner approximation has one major disadvantage: The collision rate becomes zero when the two bulk velocities are equal.For polydisperse particle size distributions this is, of course, not consistent with the original collection equation.
The true bulk collection rate will have a minimum somewhere, but it never becomes zero.One could argue that the small collision rates close to that minimum can be neglected anyway, but this may be a false conclusion.For example, when graupel grows by collection of raindrops, then the small graupel particles have fall speeds which are smaller than those of raindrops, while the large graupel may have higher fall speeds than rain, i.e., during this growth the graupel has to go through the minimum in the collision rate.If the bulk collision rate becomes zero due to the Wisner approximation, it might significantly slow down the growth of graupel because just by collecting rain it is impossible to overcome this singularity.Murakami (1990) and Mizuno (1990) suggested a new formulation of the Wisner approximation in which they replaced the difference of the weighted means in Eq. ( 15) by an ad-hoc parameterization of the velocity difference.
with α M 90 = β M 90 = 1 and γ M 90 = 0.04 as given by Murakami (1990) while Mizuno (1990) uses different values.This parameterization is, for example, used by Mansell et al. (2010) as well as Milbrandt and Yau (2005b).The disadvan-Fig.3. Terminal fall velocity of snowflakes with the mass-size relation x = 0.038D 2 s and a constant area ratio of âs = 0.45.As a reference we use the aerodynamic theory of Khvorostyanov andCurry (2002, 2005) (black solid line).The approximations are a power-law (blue dashed) and an Atlas-type relation using either equivalent diameter (red dashed) or the maximum dimension (orange dashed).
Table 1.Coefficients for the mass-size relation x = âD b, the maximum dimension as a function of particle mass D = ax b , the power-law terminal fall velocity v(x) = αx β of particles with mass x and the shape parameters of the gamma distributions f (x) = Ax ν exp(−Bx ξ ) and f (D eq ) = N 0 D µ eq exp(−λD eq ).Note that for the raindrops only the spherical geometry is given here, but the non-spherical correction is taken into account explicitly by D r,max = D r exp(ω r D r ).For the area-size relation we give the pre-factor γ A in the formula A = γ A D 2 , where A is the cross-sectional area and D is the maximum dimension.

Variance approximation
The Wisner approximation has one major disadvantage: the collision rate becomes zero when the two bulk velocities are equal.For polydisperse particle size distributions this is, of course, not consistent with the original collection equation.The true bulk collection rate will have a minimum somewhere, but it never becomes zero.One could argue that the small collision rates close to that minimum can be neglected anyway, but this may be a false conclusion.For example, when graupel grows by collection of raindrops, then the small graupel particles have fall speeds which are smaller than those of raindrops, while the large graupel may have higher fall speeds than rain; i.e., during this growth the graupel has to go through the minimum in the collision rate.If the bulk collision rate becomes zero due to the Wisner approximation, it might significantly slow down the growth of graupel because just by collecting rain it is impossible to overcome this singularity.Murakami (1990) and Mizuno (1990) suggested a new formulation of the Wisner approximation in which they replaced the difference of the weighted means in Eq. ( 15) by an ad hoc parameterization of the velocity difference.
with α M90 = β M90 = 1 and γ M90 = 0.04 as given by Murakami (1990), while Mizuno (1990) uses different values.This parameterization is, for example, used by Mansell et al. (2010) as well as Milbrandt and Yau (2005b).The disadvantage of this ad hoc formulation is that the necessary coefficients cannot be derived, and it is questionable whether one set of coefficients is good for all possible particle types and combinations.
Inspired by the parameterization of Murakami and Mizuno, SB2006 introduced an approximation in which they parameterize the bulk velocity difference by the square root of the second moment of the velocity differences, a quantity which can be calculated analytically.
For the SB2006 variance formulation we write the Wisner approximation as and the bulk velocity difference is parameterized as 20) , with the normalization factor N n,ij given by Here we have introduced an additional exponent m following a suggestion of Blahak (2012).SB2006 originally proposed m = 1, but this choice is not the only possible one.For the exponent m values between one and two seem reasonable.This exponent modifies the weight of both size distributions, and, e.g, by giving more weight to higher moments it essentially shifts the minimum of the collision rate along the internal coordinate (see Sect. 5 and Fig. 16 of the Supplement for further details).

Graupel-rain collection rates
To specify the integrals of the previous section for the riming rate of graupel collecting rain, we make the following assumptions First, we solve for the integral including the cross-sectional area With the bulk quantities defined as xj = L j /N j and Dj = D j ( xj ) = a j xb j j for any species j ∈ {r, g, i, s} we find where δ * g,2 is equal to δ 0 g of Eq. ( 90) of SB2006.Here (ξ ) is the gamma function which we evaluate using the method given in Press et al. (1992).With Fortran 2008, (ξ ) becomes part of the Fortran standard and optimized codes should become readily available.Nevertheless, the evaluation of the gamma function remains time consuming in these relations.In most two-moment bulk schemes the shape parameters like µ r , ν g , or ξ g and the particle properties that determine, e.g., b g or β g are constant during a simulation, and therefore parameters like δ * g,k or ϑ * p,m,k , as given by Eq. ( 33) below, have to be calculated only once at the initialization of the microphysics scheme.However, the behavior of the particle sedimentation can be improved by using a diagnostic µ-λ or shape-slope relation, i.e., by relating the shape parameters of the size distributions to the slope or mean size of the distribution (Milbrandt and Yau, 2005a;Seifert, 2008;Milbrandt and McTaggart-Cowan, 2010).Unfortunately, this makes it then necessary to recalculate the coefficients which include the gamma function during runtime.A compromise is to use the shape-slope relation only for precipitation-sized particles outside the cloud where the gravitational sorting is dominant, but to revert to a constant shape parameter inside the clouds where the size distribution is dominated by collisioncoalescence and other growth processes.By doing so, most parameters become again constant coefficients, and only a few change with time.Note that with the new parameterizations suggested here, fewer gamma functions occur in the relations compared to the SB2006 parameterization, which is based purely on power laws.
For the bulk velocity difference the assumptions lead to which is a coefficient similar, but not equal, to Eq. ( 92) of SB2006.Using these approximations the rate equations for all moments can be parameterized.For number and mass densities of rain we find .
The resulting approximations using these equations in comparison with the Wisner approximation and the equations given by SB2006 are shown in Figs. 4 and 5, in which the normalized bulk collision rates are compared with the numerical solution.Note that the normalized bulk collision rates have units of m s −1 ; i.e., they can also be interpreted as collection velocities.These normalized rates or characteristic bulk collection velocities of number and mass give a better visual impression of the agreement with the numerical solution than the collision rates themselves, because the trivial dependencies on the number or mass density and the cross-sectional area have been removed.For the numerical solution we have applied the Berry and Reinhardt (1973) higher-order finite-difference method on a logarithmic grid with 450 bins (mass doubling every 8th bin).For the particle size distribution of graupel a generalized gamma distribution is assumed with a graupel mass density L g = 1 g m −3 and shape parameters ν g = 1 and ξ g = 1.
For rain we assume a gamma distribution in equivalent diameter with µ r = 2 and L r = 1 g m −3 .The Figs. 4a and 5a show the well-known problems of the Wisner approximation, namely that the collision rate becomes zero where the correct solution has its local minimum.This is rectified by the variance approximation, and visually the SB2006 formulation is superior to the Wisner approximation using power-law fall speeds.Using the more accurate Atlas-type approximation of the terminal fall velocity of raindrops results in a much better agreement with the numerical solution, especially for large mean raindrop diameters.This is the case for both the Murakami and Mizuno (Figs. 4c and  5c) as well as the variance approximation (Figs.4d and 5d).Note that this improvement is only achieved by taking into account both the Atlas-type fall speed and the non-spherical geometry of the raindrops, while the approximation that apply power-law fall speeds gives better results when combined with a purely spherical geometry as already pointed out by Straka and Gilmore (2006).This error compensation is probably not coincidence, but reflects the consistency of the assumptions.The agreement of the new variance formulation with the numerical solution is very good for the mass collision rate, but larger errors do occur for the number rate, as is also the case for the other approximations.For both collision rates the optimal values for the tuning parameter, m = 1.6 for the mass rate and m = 2 for the number rate, have been used for the new variance approximation (cf.Figs. 13 and 14 of the Supplement and corresponding text).The parameterization of Murakami and Mizuno gives a result for the mass rate which is almost as good as for our new scheme, but is slighly worse for the number rate.Note that we did not retune and optimize the parameters of the Murakami-Mizuno ansatz, but simply applied the values of Murakami (1990).A quantitative error analysis is discussed in Sect.8.

Snow-rain collection rates
The collection rate of snow and rain is an example of the case of both terminal fall velocities being approximated by Atlas-type relations.For the raindrops we do again take into account the non-spherical correction, but for the snowflakes this is already included in the x s ∼ D 2 s relationship.Combined with two gamma distributions in equivalent diameter the assumptions are A. Seifert et al.: Approximation of bulk collision rates 11 a) Wisner approx.using power-law fall speed b) SB2006 using power-law fall speed c) Murakami-Mizuno approx.using Atlas-type fall speed d) Variance approx.using Atlas-type fall speed Fig. 4. Normalized mass collision rate for graupel and rain using different approximations (dashed) compared to a numerical solution of the collision integral (solid) as a function of the raindrop mean volume diameter for different mean volume diameters of the graupel size distribution.
tational expense might be acceptable.The error of the approximations as measured by the symmetric mean absolute percentage error (SMAPE) is in general below 10 %.Given the numerous uncertainties and assumptions in such schemes like particle geometries, terminal fall velocity, collision and sticking efficiencies, particle size distributions etc., this error seems acceptable.
To achieve the best possible result for a specific collection rate a calibration of the ansatz using the exponent m, cf.Eq. ( 20), is necessary, but in most cases m = 1.5 for mass and m = 2 for number rates of the interaction of two different species, and m = 1 for selfcollection provides a good f s (D eq ) = N 0,s D µ s eq exp(−λ s D eq ) (39) x s = a s D 2 s (40) Note that both D eq and D r are equivalent diameters, i.e., the diameter of a liquid water sphere.Here D eq is used for the equivalent diameter of snowflakes and D r for the equivalent diameter of raindrops, while D s and D r,max are the corresponding maximum dimensions. A. Seifert et al.: Approximation of bulk collision rates a) Wisner approx.using power-law fall speed b) SB2006 using power-law fall speed c) Murakami-Mizuno approx.using Atlas-type fall speed d) Variance approx.using Atlas-type fall speed approximation.
Whether the analytic approximations resulting from the variance approach are computationally more efficient than a look-up table will depend on the specific implementation and even more on the detailed architecture, cache size and memory bandwidth of the processor, i.e., whether the perfor-mance of the model is already limited by memory bandwidth.In such cases the look-up table leads to additional memory access which can become a serious performance issue and the additional floating point operations necessary for the analytic equations of the variance approximation can then be the smaller computational burden.With these assumptions we find the bulk collision rates for the number density How the revised and improved formulations of the terminal fall velocities and the bulk collision rates affect the simulations of clouds and precipitation will be investigated in a future study.
The Fortran code to evaluate, test and optimize the variance approximation of the bulk collision integrals for a given set of particle assumptions is available from the corresponding author.It is planned to make this code publicly available as part of the UCLA-LES code at http://gitorious.org/uclales when the new relations have been implemented and tested in the SB two-moment microphysics scheme.Fig. 6.Normalized mass collision rate for snow and rain using different approximations (dashed) compared to a numerical solution of the collision integral (solid) as a function of the raindrop mean volume diameter for different mean volume diameters of the snow size distribution.
Figures 6 and 7 compare the four different parameterizations for the bulk collision rates.For both rain and snow, we assume a gamma distribution in equivalent diameter with a mass density of L r = L s = 1 g m −3 , and the shape parameters are µ r = 2 and µ s = 2.For the collision rate of rain and snow the standard Wisner approximation with power-law fall gives reasonable results except for the minimum of the collision rate which occurs for drizzle drops.This is fixed by the SB2006 approximation, which gives a good approximation of the whole size range.The Murakami-Mizuno approximation with Atlas-type fall speeds can improve the collision rates for the very large raindrops when compared to SB2006, but suffers from the underestimation of the minimum.The new variance approximation is therefore the best parameterization as it combines both improvements over the classic Wisner approximation.To achieve the good agreement, the calibration exponent m is necessary; here we find m = 1.5 for mass and m = 2 for number.
7 Self-collection of snow Self-collection rates, i.e., loss of particles due to collisions with particles of the same species, cannot be parameterized by the Wisner approximation.Therefore most doublemoment schemes revert to look-up tables and some apply the analytic solution of Passarelli (1978) or the more general one of Verlinde et al. (1990).The analytic solutions are restricted to power-law fall speed relations and include the Gaussian hypergeometric function, which makes them expensive in the case of a time-dependent shape parameter µ s .
With the assumptions f s (D eq ) = N 0,s D µ s eq exp(−λ s D eq ) is relatively simple and computationally efficient.A comparison with the numerical solution of the integral for L s = 1 g m −3 is given in Fig. 8 for three values of the shape parameter µ s .This shows that the parameterization provides a good but not perfect approximation and is able to capture the dependency on the shape parameter µ s correctly.For the calibration exponent we have chosen m = 1.Taking into account the Atlas-type fall speed relation is clearly superior to the simple power-law fall speed of the SB2006 formulation which is shown for comparison.The numerical solution using power fall speeds should be very close to the analytic solutions of Passarelli (1978) and Verlinde et al. (1990) and provide also the reference for the SB2006 approximation.This shows that the main advantage of the refined variance approximation comes from the use of the Atlas-type velocity relations.
Fig. 8. Normalized bulk number collision rate in m s −1 for the self-collection of snow as a function of the mean maximum diameter of snowflakes.Shown are the results for three different values of the gamma shape parameter, µ s = 2, µ s = 6 and µ s = 10, for the SB2006 approximation (dotted), the revised variance approximation using the Atlas-type terminal fall velocity (dashed) and the numerical solution of the integral with Atlas-type (solid) and powerlaw (dash-dotted) fall speed.

Other collection rates
The previous sections have shown three examples of hydrometeor collision rates which include non-spherical particles, and applying an Atlas-type fall speed approximation and the variance formulation of the differential fall velocity seems to be an appropriate parameterization.This basically applies to all collision rates that include either raindrops or snowflakes.Further examples are discussed in the Supplement.Here we show some quantitative error measures which help to summarize the quality of the various approximations for such collision interactions.As error measures we have chosen the root mean square error (RMSE) and the symmetric mean absolute percentage error (SMAPE) Here φ can either be N for the number rate or L for the mass rate, and i, j identify the chosen binary collision interaction.Using a simple relative error instead of SMAPE would Whether the analytic approximations resulting from the variance approach are computationally more efficient than a look-up table will depend on the specific implementation and even more on the detailed architecture, cache size and memory bandwidth of the processor, i.e., whether the performance of the model is already limited by memory bandwidth.In such cases the look-up table leads to additional memory access, which can become a serious performance issue, and the additional floating point operations necessary for the analytic equations of the variance approximation can then be the smaller computational burden.
How the revised and improved formulations of the terminal fall velocities and the bulk collision rates affect the simulations of clouds and precipitation will be investigated in a future study.
The Fortran code to evaluate, test and optimize the variance approximation of the bulk collision integrals for a given set of particle assumptions is available from the corresponding author.It is planned to make this code publicly available as part of the UCLA-LES code at http://gitorious.org/uclales when the new relations have been implemented and tested in the SB two-moment microphysics scheme.

Fig. 1 .
Fig.1.Parameterization relations for non-spherical geometry of raindrops.The axis ratio as a function of equivalent diameter (left) following from Eq. (2) with Dη = 5.2 mm, the results of theBeard and Chuang (1987) model and the polynomial fit ofChuang and Beard (1990).For the maximum diameter as a function of the equivalent diameter (right) the explicit relation Dr,max = Dr exp(ωrDr) with ωr = 33 m −1 is used as an approximation to the implicit relation given by Eq. (2).
Fig. 2. Terminal fall velocity of raindrops as function of equivalent diameter using different approaches.The empirical relation of Beard (black solid line) is regarded as the reference.The aerodynamic theory of Khvorostyanov Curry (2002) can explain the empirical behavior, but significant differences exist (grey dashed line).The power law approximation (blue dotted line) is in general inappropriate, but the Atlas-type relation (red dash-dotted line) gives a good approximation for raindrops larger than 0.1 mm.
Fig. 3. Terminal fall velocity of snowflakes with the mass-size relation x = 0.038D 2s and a constant area ratio of âs = 0.45.As reference we use the aerodynamic theory ofKhvorostyanov and Curry (2002, 2005) (black solid line).The approximations are a power law (blue dashed) and an Atlas-type relation using either equivalent diameter (red dashed) or the maximum dimension (orange dashed).

Fig. 4 .
Fig. 4.Normalized collision rate for graupel and rain using different approximations (dashed) compared to a numerical solution of the collision integral (solid) as a function of the raindrop mean volume diameter for different mean volume diameters of the graupel size distribution.

Fig. 5 .
Fig. 5.As Fig. 4, but for the normalized number collision rate of graupel and rain.

Fig. 5 .
Fig. 5.As Fig. 4 but for the normalized number collision rate of graupel and rain.

Fig. 6 .
Fig.6.Normalized mass collision rate for snow and rain using different approximations (dashed) compared to a numerical solution of the collision integral (solid) as a function of the raindrop mean volume diameter for different mean volume diameters of the snow size distribution.

Fig. 7 .
Fig. 7.As Fig. 6 but for the normalized number collision rate of snow and rain.

Fig. 9 .
Fig. 9. Root mean square error (RMSE, left) and symmetric mean absolute percentage error (SMAPE, right) of the normalized bulk number and mass collection rates for seven different binary collision interactions.