Formulation and test of an ice aggregation scheme for two-moment bulk microphysics schemes

A simple formulation of aggregation for twomoment bulk microphysical models is derived. The solution involves the evaluation of a double integral of the collection kernel weighted with the crystal size (or mass) distribution. This quantity is to be inserted into the differential equation for the crystal number concentration which has classical Smoluchowski form. The double integrals are evaluated numerically for log-normal size distributions over a large range of geometric mean masses. A polynomial fit of the results is given that yields good accuracy. Various tests of the new parameterisation are described: aggregation as stand-alone process, in a box-model, and in 2-D simulations of a cirrostratus cloud. These tests suggest that aggregation can become important for warm cirrus, leading even to higher and longer-lasting in-cloud supersaturation. Cold cirrus clouds are hardly affected by aggregation. The collection efficiency is taken from a parameterisation that assumes a dependence on temperature, a situation that might be improved when reliable measurements from cloud chambers suggests the necessary constraints for the choice of this parameter.


Introduction
Cirrus clouds, in particular at temperatures higher than −40 • C, often contain very large ice crystals with maximum dimensions exceeding 1 mm (Heymsfield and McFarquhar, 2002, Fig. 4.6).These large crystals generally have complex shapes (Field and Heymsfield, 2003, Fig. 3), and many of them seem to be aggregates of simpler crystals, although one has to be careful in identifying irregular crystals with ag-gregates (Bailey and Hallett, 2009).But also in cold cirrus clouds (T < −40 • C) aggregated ice crystals can be found (e.g., Kajikawa and Heymsfield, 1989;Connolly et al., 2005;Bailey and Hallett, 2009), indicating that aggregation might also play a role for the cold temperature regime.
The process of ice aggregation was already investigated in the 19th century.From September 1842 on, Faraday made a series of experiments in order to investigate the ability of ice to stick onto other ice particles (Faraday, 1859), which was called "regelation" by Tyndall (1857).During this time, the accepted explanation was the so-called pressure melting proposed by Thomson (1859Thomson ( , 1860)); the main idea is that sufficient compressive forces exist at the contact region, causing melting if the ice particles are brought together.However, results by Nakaya and Matsumoto (1954) show that the required pressures are far to high to be realistic under atmospheric conditions.Faraday (1859) proposed the existence of a so-called "liquid-like" layer at the ice surface, which solidifies in case of contact with another piece of ice.This approach was supported about 100 yr later by Weyl (1951) and Fletcher (1962).Additionally, measurements by Nakaya and Matsumoto (1954) and Hosler and Hallgren (1960) indicate a temperature dependence of the sticking ability, which could be explained by the liquid layer on top of the ice crystals.Kingery (1960) proposed a different way to explain ice aggregation, namely ice sintering.Two (spherical) ice particles attach at a single point, which is not a thermodynamically stable state; in order to minimize surface free energy, a neck between the spheres is formed, thus the two particles stick together.This process of ice sintering, of course, would be supported by the "liquid-like" layer, as proposed.The details of E. Kienast-Sjögren et al.: Ice aggregation in two-moment schemes this process are discussed by Hobbs (1965).From measurements by Kumai (1964) (image reprinted in Hobbs, 1965) and Hobbs and Mason (1964) the existence of small quasispherical ice crystals (diameter ∼ 10 µm) sticking together is evident, even at cold temperatures down to −37 • C or even below.For large ice particles mechanical interlocking (Jiusto and Weickmann, 1973) may play a role, especially for large dendritic snow flakes.In summary, it is likely that ice sintering in combination with the "liquid-like" layer is the main process for ice aggregation at small sizes; the modern term for "liquid-like" layer is quasi-liquid-layer (QLL).At high temperatures near melting point or even down to −15 • C, the existence of QLLs is quite evident (Kahan et al., 2007;Sazaki et al., 2012); however, it is not clear if, for the cold temperature regime T < −40 • C, the QLL still exists; also from a theoretical point of view, this question is still undecided (Ryzhkin and Petrenko, 2009).
Although details are still unclear, it is evident that aggregation can only occur when ice crystals collide.Collisions can be caused by a variety of processes, for instance turbulent motions and gravitational settling of crystals.For crystals larger than a few µm, gravitational settling is the most efficient aggregation process (Jacobson, 2005, Fig. 15.7).Turbulent fluctuations in clouds can however enhance the gravitational aggregation process as a result of synergy between dynamics and microphysics (Sölch and Kärcher, 2011): cirrus clouds developing in a steady uplift situation have a thin nucleation zone at their top.New crystals form there as soon as the ongoing cooling drives the relative humidity over the nucleation threshold.The number of new ice crystals is a strong function of dS i /dt, the rate at which the supersaturation increases at the threshold.Turbulent motions lead to variations in dS i /dt, thus, in consequence to variations in the number concentration of new crystals.If dS i /dt is by chance particularly low, only few crystals form and they grow subsequently in highly supersaturated air with only weak competition for the excess water vapour.Thus they first grow large by deposition, obtain large fall speeds, and can then collect many ice crystals on their way from cloud top to base.
The dominance of gravitational collection has some consequences: (i) the importance of aggregation decreases with altitude (thus with decreasing temperature) because the absolute humidity decreases (roughly exponentially) and therefore mean crystal dimensions decrease with altitude; (ii) aggregation is more important in deep than in shallow (ice) clouds; (iii) aggregation is more important in well-developed than in young (ice) clouds.
Since for atmospheric investigations the evolution of an ensemble of ice particles must be evaluated, the stochastic collection equation must be investigated.This type of equation was already investigated by Smoluchowski (1916Smoluchowski ( , 1917) ) in order to describe Brownian coagulation analytically.For investigating the evolution of a size distribution, the stochastic collection equation must be solved or treated in a numerical way.In former studies, different treatments of modelling ice aggregation has been obtained.Some authors have simulated aggregation via spectral models (Khain and Sednev, 1996;Cardwell et al., 2003) or even by single particle tracking (Sölch and Kärcher, 2011).Field and Heymsfield (2003) believe that size distributions of ice crystals are dominated by depositional growth for small particles (e.g. up to 100 µm) and dominated by aggregation for larger particles which they underpin by demonstrating that the crystal size distributions for large crystals display a scaling behaviour."Scaling" is a modern expression for the attainment of a self-preserving size distribution (SPD) (see Pruppacher and Klett, 1997, ch. 11.7.2, see also Sect.4.4): the SPD theory suggests that the process of coagulation makes a crystal population loose memory of its initial size distribution and attaining asymptotically a size distribution of a relatively simple form.The further evolution of the latter with time can be described simply by scaling transformations, that is, when the x (size) and y (number) axes are transformed with two simple functions of time (x (t) = x f x (t), y (t) = y f y (t)), the size distribution is represented by a constant curve in this changing coordinate system.Such scaling behaviour in ice clouds has been demonstrated by several researchers and traced back to a dominance of aggregation processes (e.g.Westbrook et al., 2007).
This behaviour is also a kind of justification for the modelling aggregation processes via bulk models, i.e. using (fixed) size-distributions and predicting general moments such as number and mass concentration, respectively.This was done in some former studies (see, e.g., Passarelli, 1978;Lin et al., 1983;Mitchell, 1988;Levkov et al., 1992;Ferrier, 1994;Lawson et al., 1998;Field and Heymsfield, 2003;Morrison et al., 2005).However, the main focus in all these model studies as well as in most observational studies (see, e.g., Field et al., 2006;Connolly et al., 2012) is on the "high" temperature regime, i.e.T > −30 • C, where precipitation is mainly formed via the ice phase.There are only few observations in the cold temperature regime, mostly in convective outflow cirrus clouds (Connolly et al., 2005), indicating that aggregation of ice crystals happens.However, even in cold cirrus clouds formed in situ, aggregates are sometimes found (Kajikawa and Heymsfield, 1989), thus aggregation might occur at these temperature and might have an impact on microphysical properties.
The bulk ice microphysics scheme by Spichtinger and Gierens (2009a) so far did not represent ice aggregation; as only cold cirrus clouds have been simulated, this was tolerable.However, as seen above, aggregation can be an important process and a complete cirrus microphysics scheme must have a representation of it.This new treatment is described in this study.It has to be noted that although the scheme can treat multiple classes of ice (e.g.ice formed by homogeneous nucleation and ice formed by heterogeneous nucleation), it is not yet possible to compute aggregation between these different classes.Therefore we describe here aggregation between crystals of a single class of ice.

Mathematical formulation of aggregation for two-moment bulk microphysics schemes
Bulk microphysics schemes do not explicitly resolve the size spectrum of the modelled hydrometeors as spectral models (also known as bin models) or even models following single particles do.Instead, bulk models use only some low order moments of the a priori assumed size distribution and predict their temporal evolution subject to microphysical processes as nucleation, depositional or condensational growth (and evaporation or sublimation), sedimentation, and aggregation.All these processes can be formulated by first considering the process for a single particle (or two particles in the case of aggregation) and then computing integrals over the assumed size distribution.Two-moment schemes consider the evolution of the zeroth and another low-order moment, which are proportional to the particle number concentration (N) and mass concentration (q c ).We use a crystal mass distribution f (m) instead of a size pdf, thus we use the zeroth and first moment of f (m).In the following we present the theory for an assumed mass distribution, which we use in two versions, namely n(m)dm is the number concentration of particles having masses between m and m+dm, and f (m)dm is the normalised version of this, namely f (m) = n(m)/N where N = n(m)dm (the total number concentration irrespective of particle mass).This implies f (m)dm = 1.Evidently, aggregation does not change the mass concentration, thus the prognostic equation for q c is simply As this paper deals with aggregation alone, we will drop in the following the lower index referring to the process.In order to formulate the differential equation for N we start by writing down the following master-equation (e.g.Pruppacher and Klett, 1997) Here, K(m, m ) is the so-called coagulation kernel (i.e. the rate at which the crystal concentration changes due to aggregation per unit concentration of crystals of mass m and per unit concentration of crystals of mass m ).The first rhs integral describes the formation of particles of mass m from aggregation of two smaller particles, and the second rhs integral describes the aggregation of particles of mass m with other crystals of arbitrary mass, which leads to a loss of particles of mass m.Note that we can extend the 1st integral upper limit to infinity without changing its value because n(m−m ) is zero for negative arguments.This fact will be used below.
As stated above, in order to find the prognostic equation for N we have to integrate, i.e.

∂N (t
It is easy to see that every combination of m and m in the first integral occurs as well in the second one.Thus, both integrals are equal (I 1 = I 2 ), apart from the factor 1/2, so that the result is The derivation of this result is as follows: first we interchange the order of integration in the first integral (I 1 ), resulting in Now we substitute x for m − m in the inner integral, with dx = dm.The integral limits can still be set to zero and infinity, because n vanishes for negative arguments.Therefore Since it does not matter whether we write the integration variable as x or as m and because the integrand is symmetric in its two variables, we see that I 1 equals the second integral from above, and this completes our proof of Eq. ( 3).Now we go on using the normalised mass distribution.In this form, Eq. (3) reads In a mathematical sense, the double integral over the kernel function is nothing else than its expectation value for the given distribution f (m, t).This is usually notated as K (t) where we have retained the time dependence for clarity.The prognostic equation for N (t) is therefore A similar formulation was already given by Murakami (1990), however for the conversion from ice to snow and without the statistical interpretation of K .Assuming that K (t) is a constant during a single time step t in the bulk scheme, there is a formal analytical solution of the form where t * is any appropriate time within the time step.The form of this solution has already been obtained by Smoluchowski (1916Smoluchowski ( , 1917) ) for Brownian coagulation (see also Pruppacher and Klett, 1997, Sect. 11.5).It is seen that for long timesteps such that t 1/[N(t) K (t * )] the final N(t + t) becomes independent of the initial N(t).With the new value of N and the (here unchanged) value of q c we can compute the updated mean mass of f (m, t + t), i.e. we can compute the updated mass distribution.
The prognostic equation for N is the desired result, and the solution can in principle be computed for arbitrary forms of the mass distribution and for arbitrary coagulation kernels.Nevertheless, the necessary integrations are tedious and it may be justified to construct a look-up table where the required values of K (t) can be read off.The computation of the integrals for special choices of f (m) and the kernel is demonstrated next.

Choice of a mass distribution
In principle we could use any probability density function on R + for f (m).Following Spichtinger and Gierens (2009a), we use here a log-normal distribution, i.e.
Here, log denotes natural logarithm.The normalised f (m) has two parameters.The first parameter is the modal mass (or geometric mean) m m , which is updated after every time step by the prognostic values of number and mass concentration, respectively; the second parameter is the geometric standard deviation σ m , which is usually fixed or formulated as a function of the mean mass.The mass distribution is then given by normalisation with number concentration, i.e. n(m) = N •f (m).The general kth moment of the mass distribution is denoted by µ k [m] and for log-normal distributions we obtain Number and mass concentration are prognostic variables in our scheme, represented by the general moments , respectively, with a mean mass of m = q c /N = m m exp 1 2 (log σ m ) 2 .Aggregation would tend to lead to a deviation from the log-normal distribution towards an exponential one.This effect cannot be taken into account in our model and would require some development, for instance introduction of an ice class "aggregates" with exponential distribution and development of an aggregation scheme between small ice crystals (log-normally distributed) and aggregates.All this is future work.

Choice of an aggregation kernel
We assume that ice crystals aggregate in particular when large crystals fall through an ensemble of small crystals, when they collide and stick together.This particular mechanism is called gravitational collection and can be described by the following form of a collection kernel (Pruppacher and Klett, 1997, p. 569) Here, R and r are the "radii" (see below) of the larger and smaller colliding ice crystals, respectively, such that the first factor on the rhs is the geometric cross section for the collision.The second factor is the absolute difference of the fall speeds of the two crystals, that is, the speed of their relative motion.Because of hydrodynamic (and potentially other) effects it is not just the geometric cross section that determines whether two crystals collide, and even if they collide they do not need to stick together.Therefore a correction factor E(R, r) is applied which accounts for these effects.E is usually called collision or collection efficiency.Choices of E will be presented below.
The next problem to solve is to formulate the collection kernel for non-spherical ice crystals instead of spheres.Here we do this for hexagonal cylinders, a common shape for ice crystals as used for instance in Spichtinger and Gierens (2009a).For convenience, the size is replaced by the particle mass.This can be done since mass (m) and size (L) are related (see e.g.Heymsfield and Iaquinta, 2000), usually expressed via power laws (e.g.m = αL β ).Using these relations, we obtain the following expression for the surface of a hexagonal ice crystal of mass m where ρ b = bulk density of ice = 0.81 × 10 3 kg m −3 .Assuming randomly oriented columns (analogous to the usual approximation for radiation parameterisation in Ebert and Curry, 1992) we obtain r 2 = A(m) 4π by replacing the ice crystals surface by the surface of a sphere, such that We replace the radii (R, r) with the corresponding masses (M, m) and obtain The terminal velocity of each of the falling particles can be described using the following power law The parameters α, β, γ and δ needed to translate K(R, r) into K(M, m) are given in Spichtinger and Gierens (2009a).Note that the parameters usually are constants for values of m in a certain mass interval.Thus, this leads to a generic splitting of the integrals, as can be seen in the next section.A correction factor for the terminal velocity is added in order to consider density changes, leading to different aerodynamic drag.The correction factor is represented as follows (Heymsfield and Iaquinta, 2000) with constants p 0 = 300 hPa, T 0 = 233 K, a = −0.178,b = −0.394.Since the correction factor depends only on temperature and pressure, respectively, it can be treated as a constant for the calculations of the integrals.

Computation of the integrals
The integral values ( K in m 3 s −1 ) were calculated using an adaptive Simpson quadrature (e.g.Lyness, 1969) for different values of the geometrical standard deviation σ m .The calculated integral results for σ m = 2.85 are indicated as squares in Fig. 1.The calculation was done for T = 233 K and 300 hPa, i.e. c(T , p) = 1 in Eq. ( 14).For other choices of atmospheric parameters the values change moderately.For reasons given below we compute the integrals with a collision efficiency of unity.The integration was conducted up to a modal mass of 10 6 ng = 10 −6 kg (equivalent to a length ∼ 8 mm) as this is the upper limit for an aggregated particle, which will be used in the later parameterisation.The calculated integral values were divided into four ranges because of mass dependent coefficients α, β, γ and δ.For each range a polynomial was fitted through the calculated values.The combination of these polynomials is displayed as solid line in Fig. 1.The four polynomial ranges correspond to changes in the growth and sedimentation behaviour of ice crystals (e.g.changes of the parameters α, . . ., δ) as indicated in Spichtinger and Gierens (2009a): 1. 1 × 10 −4 ng < m m ≤ 2.5 × 10 −3 ng: mainly hexagonal ice crystals with aspect ratio 1.
E.Kienast-Sjögren et al.: Ice aggregation in 2-moment schemes 5 We replace the radii (R,r) with the corresponding masses (M,m) and obtain The terminal velocity of each of the falling particles can be described using the following power law The parameters α, β, γ and δ needed to translate K(R,r) into K(M,m) are given in Spichtinger and Gierens (2009a).Note that the parameters usually are constants for values of m in a certain mass interval.Thus, this leads to a generic splitting of the integrals, as can be seen in the next section.
A correction factor for the terminal velocity is added in order to consider density changes, leading to different aerodynamic drag.The correction factor is represented as follows (Heymsfield and Iaquinta, 2000) with constants p 0 = 300hPa, T 0 = 233K, a = −0.178,b = −0.394.Since the correction factor depends only on temperature and pressure, respectively, it can be treated as a constant for the calculations of the integrals.

Computation of the integrals
The integral values ( K in m 3 s −1 ) were calculated using an adaptive Simpson quadrature (e.g.Lyness, 1969) for different values of the geometrical standard deviation σ m .The calculated integral results for σ m = 2.85 are indicated as squares in Figure 1.The calculation was done for T = 233K and 300hPa, i.e. c(T,p) = 1 in eq. ( 14).For other choices of atmospheric parameters the values change moderately.For reasons given below we compute the integrals with a collision efficiency of unity.The integration was conducted up to a modal mass of 10 6 ng = 10 −6 kg (equivalent to a length ∼ 8mm) as this is the upper limit for an aggregated particle, which will be used in the later parameterization.The calculated integral values were divided into four ranges because of mass dependent coefficients α, β, γ and δ.For each range a polynomial was fitted through the calculated values.The combination of these polynomials is displayed as solid line in Figure 1.The four polynomial ranges correspond to changes in the growth and sedimentation behaviour of ice crystals (e.g.changes of the parameters α,...,δ) as indicated in Spichtinger and Gierens (2009a) 1. 1•10 −4 ng < m m ≤ 2.5•10 −3 ng: Mainly hexagonal ice crystals with aspect ratio 1.As can be seen in Figure 1, the polynomial fits to the exact integral values are very good; the maximum deviation is less than 10% and the mean deviation even less than 2%.Thus the polynomials can be used as an accurate solution while saving computing time.These fits were calculated for a whole range of lognormal size distributions with geometric standard deviations in the range σ m = 1.90 − 5.29.The coefficients for these fits are given in the appendix, table 2. In figure 2 the kernels for these vales are displayed against modal mass m m of the distribution (upper panel) as well as against mean mass of the distribution m = m m exp(0.5 • (logσ m ) 2 ) (lower panel).As expected by theory, the kernels are larger for a wider distribution; this behaviour can be seen clearly in the upper panel of fig. 2. In fact, the kernels for a fixed modal mass vary over more than one order of magnitude for different geometric standard deviations.This indicates that the width of the size distribution is very important for the strength of the aggregation process.

Comparison with other model parameterisations
For testing our new parameterisation, we first compare the derived kernels with already existing parameterisations.There were some former attempts in deriving aggregation parameterisations for hydrometeors, especially snow and ice crystals.Passarelli (1978) derived an analytical kernel for 3. 4×10 2 ng < m m ≤ 1×10 4 ng: the sedimentation velocity gets larger.
As can be seen in Fig. 1, the polynomial fits to the exact integral values are very good; the maximum deviation is less than 10 % and the mean deviation even less than 2 %.Thus the polynomials can be used as an accurate solution while saving computing time.These fits were calculated for a whole range of log-normal size distributions with geometric standard deviations in the range σ m = 1.90-5.29.The coefficients for these fits are given in the appendix, Table A1.In Fig. 2 the kernels for these vales are displayed against modal mass m m of the distribution (upper panel) as well as against mean mass of the distribution m = m m exp(0.5 • (log σ m ) 2 ) (lower panel).As expected by theory, the kernels are larger for a wider distribution; this behaviour can be seen clearly in the upper panel of Fig. 2. In fact, the kernels for a fixed modal mass vary over more than one order of magnitude for different geometric standard deviations.This indicates that the width of the size distribution is very important for the strength of the aggregation process.

Comparison with other model parameterisations
For testing our new parameterisation, we first compare the derived kernels with already existing parameterisations.There were some former attempts in deriving aggregation parameterisations for hydrometeors, especially snow and ice crystals.Passarelli (1978)  Since the mean mass depends on the geometric width of the distribution, the variation of the kernels due to different width is smaller.Additionally, the aggregation kernel as derived by Schumann (2012) is shown for comparison.The corresponding sizes range between a modal length of 0.6 and 8000 µm.
aggregating ice crystals, leading to an unhandy expression of hypergeometric functions.However, he assumed an exponential size distribution, because he was interested in aggregation of snow, i.e. aggregation at warm temperatures.
The assumption of exponential distributions simplified the expressions.This procedure is not viable in our case because we use lognormal distributions for ice crystals (as justified in Spichtinger and Gierens, 2009a).Mitchell (1988) derived an aggregation kernel in a similar way as Passarelli (1978), using exponential distributions and, similar to our treatment, power laws for terminal velocities.Again, the aggregation parameterisation was derived for the warm temperature range, with the special assumption of an exponential distribution.Ferrier (1994) used an approach similarly to ours, using gamma size distribution.He evaluated the double integrals numerically in order to create a look-up table.However, again this parameterisation was made for the warm temperature regime.Many other models rely on these de-scribed parameterisations (e.g., Morrison et al., 2005, using Passarelli's parameterisation), maybe with modifications.
For the low temperature regime (T < −40C • ) only one recent parameterisation by Schumann ( 2012) is available.Schumann (2012) estimated the aggregation kernel (his equation 52, translated into our formulation) to be using the volume mean radius r = 3 m 4πρ bs with a bulk ice density of 917 kg m −3 .For comparison we used our terminal velocity formulation, the volume mean radius was derived from the mean mass of an ice crystal; the efficiency is set to be E ≡ 1.In fig. 2 (lower panel) some kernels for our formulation (σ m = 1.9/2.85/5.29,representing narrow/medium/broad size distributions) are shown in comparison with the kernel as derived by Schumann (2012).The kernels are plotted against the mean mass m.At least in the mass range 10 −14 ≤ m ≤ 10 −9 kg, the qualitative agreement is quite good, although the kernel by Schumann ( 2012) is about 5 times higher than our parameterisations.In the mass ranges m ≤ 10 −14 kg and m ≥ 10 −9 kg there is a larger overestimation compared to our parameterisations.These overestimations are not crucial, since (1) the larger masses do not apply for the parameterisation by Schumann (2012) which is used for contrails, where small to moderate crystal masses prevail, and (2) for the very small masses the aggregation time scales are much larger than cloud and contrail lifetimes (see below).

Time scale analysis
In order to estimate the possible impact of aggregation on ice crystal number concentrations, we estimate the time scales of aggregation In fig. 3 the timescales of aggregation are displayed for a typical size distribution with geometric standard deviation value of σ m = 2.85 and for typical ice crystal number concentrations in cirrus clouds (see, e.g., Krämer et al., 2009) in the range between N = 10 4 m −3 = 10L −1 and N = 10 7 m −3 = 10 cm −3 , respectively.Since aggregation is a pure sink for ice crystal number concentrations, the time scale is negative; however, in fig. 3 we show absolute values of τ for a better representation.
It is evident, that only in a very narrow range set by parameters number concentration and ice crystal size, respectively, aggregation might play a role.Since the life time of cirrus clouds might be in the order of one day (e.g.Spichtinger et al., 2005), this time interval might serve as an upper limit for the impact of aggregation on cirrus clouds.Since in clouds with only a few ice crystals (e.g.N = 10 4 m −3 ) the ice crystals can growth to large sizes -at least at high temperatures exponential size distribution, because he was interested in aggregation of snow, i.e. aggregation at warm temperatures.
The assumption of exponential distributions simplified the expressions.This procedure is not viable in our case because we use log-normal distributions for ice crystals (as justified in Spichtinger and Gierens, 2009a).Mitchell (1988) derived an aggregation kernel in a similar way as Passarelli (1978), using exponential distributions and, similar to our treatment, power laws for terminal velocities.Again, the aggregation parameterisation was derived for the warm temperature range, with the special assumption of an exponential distribution.Ferrier (1994) used an approach similarly to ours, using gamma size distribution.He evaluated the double integrals numerically in order to create a look-up table.However, again this parameterisation was made for the warm temperature regime.Many other models rely on these de-scribed parameterisations (e.g., Morrison et al., 2005, using Passarelli's parameterisation), maybe with modifications.
For the low temperature regime (T < −40 • C) only one recent parameterisation by Schumann ( 2012) is available.Schumann (2012) estimated the aggregation kernel (his Eq. 52, translated into our formulation) to be using the volume mean radius r = 3 m 4πρ bs 1 3 with a bulk ice density of 917 kg m −3 .For comparison we used our terminal velocity formulation, the volume mean radius was derived from the mean mass of an ice crystal; the efficiency is set to be E ≡ 1.In Fig. 2 (lower panel) some kernels for our formulation (σ m = 1.9/2.85/5.29,representing narrow/medium/broad size distributions) are shown in comparison with the kernel as derived by Schumann (2012).The kernels are plotted against the mean mass m.At least in the mass range 10 −14 ≤ m ≤ 10 −9 kg, the qualitative agreement is quite good, although the kernel by Schumann (2012) is about 5 times higher than our parameterisations.In the mass ranges m ≤ 10 −14 kg and m ≥ 10 −9 kg there is a larger overestimation compared to our parameterisations.These overestimations are not crucial, since (1) the larger masses do not apply for the parameterisation by Schumann (2012) which is used for contrails, where small to moderate crystal masses prevail, and (2) for the very small masses the aggregation timescales are much larger than cloud and contrail lifetimes (see below).

Timescale analysis
In order to estimate the possible impact of aggregation on ice crystal number concentrations, we estimate the timescales of aggregation In Fig. 3 the timescales of aggregation are displayed for a typical size distribution with geometric standard deviation value of σ m = 2.85 and for typical ice crystal number concentrations in cirrus clouds (see, e.g., Krämer et al., 2009) in the range between N = 10 4 m −3 = 10 L −1 and N = 10 7 m −3 = 10 cm −3 , respectively.Since aggregation is a pure sink for ice crystal number concentrations, the timescale is negative; however, in Fig. 3 we show absolute values of τ for a better representation.
It is evident that only in a very narrow range set by parameter number concentration and ice crystal size, aggregation might play a role.Since the life time of cirrus clouds might be in the order of one day (e.g.Spichtinger et al., 2005), this time interval might serve as an upper limit for the impact of aggregation on cirrus clouds.Since in clouds with only a few ice crystals (e.g.N = 10 4 m −3 ) the ice crystals can grow Fig. 3. Time scales of aggregation for a medium width of the ice crystal size distribution (σm = 2.85) and for typical ice crystal number concentrations as found in cirrus clouds at cold temperatures (see, e.g., Krämer et al., 2009).The corresponding sizes range between a modal length of 0.6 and 8000 µm.
-this regime can be effectively influenced by aggregation.Cirrus clouds containing many ice crystals are usually dominated by small crystals.Thus, although ice crystal number concentration is able to reduce the aggregation time scales by orders of magnitudes, an upper limit in crystal size due to thermodynamic constraints leads to less effective aggregation.We will see later, that this very simple estimation from time scale analysis will be coroborrated by detailed tests.

Various tests
The change in particle number density per timestep is described in (5) as The solution of this can either be achieved through an exact solution involving separation of variables with the following result or through the following Euler approximation Tests have shown that both methods give practically identical results.The following tests have been performed with the Euler approximation.Further tests have shown that the polynomial approximation of the K integrals was a sufficient approximation (see above), so the following tests have been performed using the polynomial approximation.

Maximum aggregation
The new formulation of the aggregation process was tested in MATLAB with different start values for the particle number density ranging from 10 4 to 5 • 10 4 particles per cubic metre.We let aggregation occur as the only process.We set E ≡ 1 for these tests, that is, the following results show a maximum effect of aggregation.The aggregation was run for 1000 s, i.e. approximately 17 minutes.If aggregation occurs, the particle number density N will decrease and the mean mass ( m = q c /N ) increase.Using the expressions for the log-normal distribution we compute then a new modal mass and the corresponding new f (m) is used in the next timestep.For the calculations, we set an upper boundary for the mass of the aggregated particles at 10 −6 kg, which corresponds to a particle size of 8 mm.Larger ice crystals will not occur in the model.
Figure 4 shows the results of these tests.For small initial modal masses (e.g. 10 −11 kg, green line) and starting with, for example, 10 5 particles m −3 , after simulating 1000 s we still have 10 5 particles m −3 .Thus almost no aggregation occurs.If the initial modal mass is instead increased to 10 −9 kg (blue line) while still starting with 10 5 m −3 particles, we only have about 10 4 m −3 particles left after 1000 s.Thus about 90% of the particles have aggregated.
As expected, when small crystals (i.e.small modal masses) predominate, nothing happens.The probability for collision is negligible and the relative fall speeds are low.The larger the particles get, the more they aggregate.For the largest initial modal masses (i.e. 10 −9 kg) the particles stick together very fast so that the iteration has to be stopped before reaching 1000 s.
Timesteps from 0.1 s (timescale for microphysics) to 10 s (timescale for dynamics) gave very similar results, that is, the parameterization is consistent and convergent with different timesteps.We chose to use a timestep of 1 s.

Introducing temperature dependency
Experience from field measurements suggests that aggregation occurs more efficiently in warmer than in colder air (e.g.Kajikawa and Heymsfield, 1989).This behaviour is also consistent with the possible existence of a QLL on top of ice crystals, even at low temperatures.The temperature dependence is expressed by the following parameterization for the collision efficiency of ice crystals (Lin et al., 1983;Levkov et al., 1992) which is independent of the crystal masses and dependent only on temperature (the original papers do not mention whether the parameterisation is based on measurements) With this parameterization, E can be taken out from the integral calculation and treated as a prefactor.Note here, that Fig. 3. Timescales of aggregation for a medium width of the ice crystal size distribution (σ m = 2.85) and for typical ice crystal number concentrations as found in cirrus clouds at cold temperatures (see, e.g., Krämer et al., 2009).The corresponding sizes range between a modal length of 0.6 and 8000 µm.
to large sizes -at least at high temperatures -this regime can be effectively influenced by aggregation.Cirrus clouds containing many ice crystals are usually dominated by small crystals.Thus, although ice crystal number concentration is able to reduce the aggregation timescales by orders of magnitudes, an upper limit in crystal size due to thermodynamic constraints leads to less effective aggregation.We will see later that this very simple estimation from timescale analysis will be corroborated by detailed tests.

Various tests
The change in particle number density per time step is described in Eq. ( 5) as The solution of this can either be achieved through an exact solution involving separation of variables with the following result or through the following Euler approximation Tests have shown that both methods give practically identical results.The following tests have been performed with the Euler approximation.Further tests have shown that the polynomial approximation of the K integrals was a sufficient approximation (see above), so the following tests have been performed using the polynomial approximation.

Maximum aggregation
The new formulation of the aggregation process was tested in MATLAB with different start values for the particle number density ranging from 10 4 to 5 × 10 4 particles per cubic metre.We let aggregation occur as the only process.
We set E ≡ 1 for these tests, that is, the following results show a maximum effect of aggregation.The aggregation was run for 1000 s, i.e. approximately 17 min.If aggregation occurs, the particle number density N will decrease and the mean mass ( m = q c /N) increase.Using the expressions for the log-normal distribution we compute then a new modal mass and the corresponding new f (m) is used in the next time step.For the calculations, we set an upper boundary for the mass of the aggregated particles at 10 −6 kg, which corresponds to a particle size of 8 mm.Larger ice crystals will not occur in the model.
Figure 4 shows the results of these tests.For small initial modal masses (e.g. 10 −11 kg, green line) and starting with, for example, 10 5 particles m −3 , after simulating 1000 s we still have 10 5 particles m −3 .Thus almost no aggregation occurs.If the initial modal mass is instead increased to 10 −9 kg (blue line) while still starting with 10 5 m −3 particles, we only have about 10 4 m −3 particles left after 1000 s.Thus about 90 % of the particles have aggregated.
As expected, when small crystals (i.e.small modal masses) predominate, nothing happens.The probability for collision is negligible and the relative fall speeds are low.The larger the particles get, the more they aggregate.For the largest initial modal masses (i.e. 10 −9 kg) the particles stick together very fast so that the iteration has to be stopped before reaching 1000 s.Timesteps from 0.1 s (timescale for microphysics) to 10 s (timescale for dynamics) gave very similar results, that is, the parameterisation is consistent and convergent with different timesteps.We chose to use a time step of 1 s.

Introducing temperature dependency
Experience from field measurements suggests that aggregation occurs more efficiently in warmer than in colder air (e.g.Kajikawa and Heymsfield, 1989).This behaviour is also consistent with the possible existence of a QLL on top of ice crystals, even at low temperatures.The temperature dependence is expressed by the following parameterisation for the collision efficiency of ice crystals (Lin et al., 1983;Levkov et al., 1992) which is independent of the crystal masses and dependent only on temperature (the original papers do not mention whether the parameterisation is based on measurements)  Axis) and after 1000 s of iteration (N , y-Axis).For small modal masses, the particle number density hardly changes during the iteration and the plot is almost a straight line.As reference for no aggregation, a black line is plotted.As expected, larger modal masses result in increased aggregation.Thus the particle number density decreases during integration, which is shown in the graph.The corresponding size to the modal masses plotted ranges between 6 and 345 µm.Note that the tests have been performed with E ≡ 1, that is, the maximum effect of aggregation is seen here.The red line is almost equal to the black line, thus can hardly be seen in the plot.
experimental evidence of the exact form of the temperature dependence is not given.Nevertheless, this is the only temperature dependence we found from literature, which also seems to be reasonable.The temperature dependent collision efficiency is 0.21 ≤ E(T ) ≤ 0.44 for typical temperatures in ice clouds (210 ≤ T ≤ 240K).Thus we expect reduced aggregation effects on N compared to the previous tests when we introduce the new factor.Figure 5 shows the final crystal concentrations as before for aggregation without temperature dependency (E = 1) and for different temperatures.As expected, with decreasing temperature aggregation becomes less important.Even at the highest considered temperature (240 K) the reduction of the aggregation effect is considerable, in particular for high initial number densities.This finding is a good argument for ignoring aggregation in simulations of cold cirrus clouds where not only E is small but also the median crystal masses are smaller than in warm cirrus.

Test of aggregation within a box model
For investigating the impact of aggregation idealized box model simulations were carried out.Here, we use a maximum efficiency E ≡ 1 as well as a temperature dependent efficiency E(T ) in order to see a realistic impact of aggregation on cirrus clouds.

Model description and setup
In this section we test the effect of aggregation on the ice crystal number concentration in the framework of a box model (Spichtinger and Gierens, 2009a) which we consider a more realistic test in that various microphysical processes can act simultaneously.The box is thought to represent an initially cloud free air parcel which is lifted with a constant vertical velocity.During the cooling procedure, homogeneous freezing of aqueous solution droplets (short "homogeneous nucleation", parameterized after Koop et al., 2000) will occur, i.e. ice crystals are formed.In the supersaturated environment the ice crystals grow by diffusional growth (based on approximations by Koenig, 1971) to larger sizes.The parameterisations for both processes are described in detail in Spichtinger and Gierens (2009a).For determining the impact of aggregation for different temperature and velocity regimes, many idealized box simulations were carried out.Each simulation starts at p = 300hPa.The initial temperature is given by T = 210/220/230/240K, the vertical velocity range is given by w = 0.01/0.02/0.05/0.1/0.2/0.5/1/2/5m s −1 .The total simulation time is calculated by t sim = 720m/w; this procedure ensures that at the end of the simulation the same environmental conditions (T,p) are reached for each inital Fig. 4. Particle number density N 0 at the start of the iteration (x axis) and after 1000 s of iteration (N, y axis).For small modal masses, the particle number density hardly changes during the iteration and the plot is almost a straight line.As reference for no aggregation, a black line is plotted.As expected, larger modal masses result in increased aggregation.Thus the particle number density decreases during integration, which is shown in the graph.The corresponding size to the modal masses plotted ranges between 6 and 345 µm.Note that the tests have been performed with E ≡ 1, that is, the maximum effect of aggregation is seen here.The red line is almost equal to the black line, thus can hardly be seen in the plot.
With this parameterisation, E can be taken out from the integral calculation and treated as a prefactor.Note here that experimental evidence of the exact form of the temperature dependence is not given.Nevertheless, this is the only temperature dependence we found from literature, which also seems to be reasonable.The temperature dependent collision efficiency is 0.21 ≤ E(T ) ≤ 0.44 for typical temperatures in ice clouds (210 ≤ T ≤ 240 K).Thus we expect reduced aggregation effects on N compared to the previous tests when we introduce the new factor.Figure 5 shows the final crystal concentrations as before for aggregation without temperature dependency (E = 1) and for different temperatures.As expected, with decreasing temperature aggregation becomes less important.Even at the highest considered temperature (240 K) the reduction of the aggregation effect is considerable, in particular for high initial number densities.This finding is a good argument for ignoring aggregation in simulations of cold cirrus clouds where not only E is small but also the median crystal masses are smaller than in warm cirrus.

Test of aggregation within a box model
For investigating the impact of aggregation idealised box model simulations were carried out.Here, we use a maxi-  Axis) and after 1000 s of iteration (N , y-Axis).For small modal masses, the particle number density hardly changes during the iteration and the plot is almost a straight line.As reference for no aggregation, a black line is plotted.As expected, larger modal masses result in increased aggregation.Thus the particle number density decreases during integration, which is shown in the graph.The corresponding size to the modal masses plotted ranges between 6 and 345 µm.Note that the tests have been performed with E ≡ 1, that is, the maximum effect of aggregation is seen here.The red line is almost equal to the black line, thus can hardly be seen in the plot.
experimental evidence of the exact form of the temperature dependence is not given.Nevertheless, this is the only temperature dependence we found from literature, which also seems to be reasonable.The temperature dependent collision efficiency is 0.21 ≤ E(T ) ≤ 0.44 for typical temperatures in ice clouds (210 ≤ T ≤ 240K).Thus we expect reduced aggregation effects on N compared to the previous tests when we introduce the new factor.Figure 5 shows the final crystal concentrations as before for aggregation without temperature dependency (E = 1) and for different temperatures.As expected, with decreasing temperature aggregation becomes less important.Even at the highest considered temperature (240 K) the reduction of the aggregation effect is considerable, in particular for high initial number densities.This finding is a good argument for ignoring aggregation in simulations of cold cirrus clouds where not only E is small but also the median crystal masses are smaller than in warm cirrus.

Test of aggregation within a box model
For investigating the impact of aggregation idealized box model simulations were carried out.Here, we use a maximum efficiency E ≡ 1 as well as a temperature dependent Particle number density N 0 at the beginning of the iteration (x-Axis) and after 1000 s of iteration (N, y-Axis) for a modal mass of 10 −9 kg, which corresponds to a particle with modal length 345 µm.As a reference for no aggregation, a black line is plotted.Adding a temperature dependency shows a weakened aggregation effect with decreasing temperatures.The lower the tempertures the more particles are present at the end of the simulation.
efficiency E(T ) in order to see a realistic impact of aggregation on cirrus clouds.

Model description and setup
In this section we test the effect of aggregation on the ice crystal number concentration in the framework of a box model (Spichtinger and Gierens, 2009a) which we consider a more realistic test in that various microphysical processes can act simultaneously.The box is thought to represent an initially cloud free air parcel which is lifted with a constant vertical velocity.During the cooling procedure, homogeneous freezing of aqueous solution droplets (short "homogeneous nucleation", parameterized after Koop et al., 2000) will occur, i.e. ice crystals are formed.In the supersaturated environment the ice crystals grow by diffusional growth (based on approximations by Koenig, 1971) to larger sizes.The parameterisations for both processes are described in detail in Spichtinger and Gierens (2009a).For determining the impact of aggregation for different temperature and velocity regimes, many idealized box simulations were carried out.Each simulation starts at p = 300hPa.The initial temperature is given by T = 210/220/230/240K, the vertical velocity range is given by w = 0.01/0.02/0.05/0.1/0.2/0.5/1/2/5m s −1 .The total simulation time is calculated by t sim = 720m/w; this procedure ensures that at the end of the simulation the same environmental conditions (T,p) are reached for each inital Fig. 5. Particle number density N 0 at the beginning of the iteration (x axis) and after 1000 s of iteration (N , y axis) for a modal mass of 10 −9 kg, which corresponds to a particle with modal length 345 µm.As a reference for no aggregation, a black line is plotted.Adding a temperature dependency shows a weakened aggregation effect with decreasing temperatures.The lower the temperatures the more particles are present at the end of the simulation.
mum efficiency E ≡ 1 as well as a temperature dependent efficiency E(T ) in order to see a realistic impact of aggregation on cirrus clouds.

Model description and set-up
In this section we test the effect of aggregation on the ice crystal number concentration in the framework of a box model (Spichtinger and Gierens, 2009a) which we consider a more realistic test in that various microphysical processes can act simultaneously.The box is thought to represent an initially cloud free air parcel which is lifted with a constant vertical velocity.During the cooling procedure, homogeneous freezing of aqueous solution droplets (short "homogeneous nucleation", parameterised after Koop et al., 2000) will occur, i.e. ice crystals are formed.In the supersaturated environment the ice crystals grow by diffusional growth (based on approximations by Koenig, 1971) to larger sizes.The parameterisations for both processes are described in detail in Spichtinger and Gierens (2009a).For determining the impact of aggregation for different temperature and velocity regimes, many idealised box simulations were carried out.Each simulation starts at p = 300 hPa.The initial temperature is given by T = 210/220/230/240 K, the vertical velocity range is given by w = 0.01/0.02/0.05/0.1/0.2/0.5/1/2/5m s −1 .The total simulation time is calculated by t sim = 720m/w; this procedure ensures that at the end of the simulation the same environmental conditions (T , p) are reached for each initial temperature run (or equivalently, an altitude difference of z sim = 720 m is reached).For each set-up three scenarios were calculated: without aggregation, with temperaturedependent aggregation and with maximum aggregation.For deriving the maximum impact of aggregation, the box is a closed system, i.e. ice crystals stay in the box.However, for more realistic treatment, we have to consider sedimentation.The process of sedimentation is treated in the box model as described in Spichtinger and Cziczo (2010).Here, we repeat some key features of this treatment.Time and space is discretized by time steps t and space increments z, respectively; t n := n• t.Generally, the sedimentation process for a quantity ψ (e.g.ice mass mixing ratio q c or ice crystal number concentration N) can be described in the following way with the terminal velocity v ψ for the quantity ψ; α ψ = v ψ t z denotes the respective Courant number and F top ψ is the flux through the top of the layer.For our purpose, we use fluxes for quantities ice mass mixing ratio q c and ice crystal number density N; the mass and number weighted terminal velocities are then given by the following expressions whereas n(m) denotes the mass distribution.Since F top ψ is usually unknown for a single box, we can assume, that the flux through the top is given by a fraction of the flux through the bottom, i.e., F top ψ = f sed F bottom ψ with the sedimentation parameter f sed .Using this ansatz, it is possible to categorise different sedimentation scenarios: f sed = 1.0 corresponds to no sedimentation as net effect.The flux exiting the lower part of the box has the same magnitude as the flux entering the top of the box.
f sed = 0.9 corresponds to regions in the middle and lower part of the cloud.The flux leaving the region is almost but not completely balanced by the flux entering that region from above.
f sed = 0.5 corresponds to the top region of a cloud.The flux leaving that region is only half balanced by a flux from above.
These three scenarios are used in the simulations in order to investigate the impact of aggregation under different, but realistic conditions.
In the following we discuss a typical scenario of a steady updraught of w = 0.05 m s −1 at high temperatures (i.e. initial temperature T = 240 K).

Results for an initial temperature of 240 K
As a baseline experiment we consider first a case without the effect of sedimentation, that is, with f sed = 1.The temporal variation of the crystal number density in the box is shown in Fig. 6.As the box is cooled down from 240 K, the threshold supersaturation for homogeneous nucleation is reached about 120 min later, and the nucleation burst leads to a high crystal concentration (≈ 2 × 10 4 m −3 ).The further development depends strongly on whether we include aggregation or not.Without aggregation, the crystal concentration stays essentially constant (a weak reduction is caused by the expansion of the lifting box).With aggregation, however, the crystal concentration decreases strongly (by about 95 %) within two hours.In the scenario with a temperature-dependent aggregation the reduction of the number concentration is less drastic, but still of the order 85 %.This happens in this academic case because the large crystals effectively stay within the box (the crystals leaving the box are replaced by identical crystals entering from above) and aggregate over the whole simulation time.
Now we turn sedimentation on, allowing the large crystals to leave the box without complete replacement from above.For the top region of the clouds (f sed = 0.5), the results are displayed in Fig. 6.Again, nucleation occurs after 120 min and a large number of ice crystals appear.These grow by vapour deposition and RH i decreases.Sedimentation is a much more important process now than aggregation, which can be seen from the timescale of the decrease of N which occurs much faster than in the previous case where sedimentation was switched off.The two curves representing the cases with and without aggregation are almost identical, also after further nucleation bursts occur.This is possible because RH i starts to increase again because of the ongoing cooling when N is sufficiently diminished.
In the middle of the cloud, where a good deal of falling ice crystals are replaced by crystals falling from above (f sed = 0.9), aggregation is a bit more important than at the top of the cloud.This is shown in Fig. 6.We see that the reduction of N after the nucleation burst is slower than at the top of the cloud because more falling crystals are replaced.Aggregation accelerates the reduction of N such that a certain level of N is now reached some ten minutes earlier in the case with aggregation than without.The scenario with temperature-dependent aggregation lies between the two others.Secondary nucleation bursts do not occur in either case in the middle of the cloud until the end of the simulation.peratures, for all choices of f sed and both aggregation scenarios.This is indeed what we found from our box model simulations: the curves for the cases with and without aggregation differ not drastically, thus we refrain from showing the simulations in detail.

Maximum impact of aggregation as derived from box model simulations
For determining the maximum impact of aggregation for different temperature and velocity regimes, we investigate the box model simulations with sedimentation factor f sed = 1, i.e. no sedimentation is allowed.In order to determine the maximum impact of aggregation we define an aggregation factor f agg : The factor is given by the ratio of ice crystal number concentrations of the aggregation run and the reference run (without aggregation): with t =end of simulation, and for temperature-dependent simulations the analogon is formed.In fig.7 the aggregation factor is shown for maximum aggregation (top panel) and temperature-dependent aggregation (bottom panel).
Figure 7 can be interpreted as follows: -For low vertical velocities at warm temperatures, homogeneous nucleation does produce just low number concentrations (see, e.g., fig.12 in Spichtinger and Gierens, 2009a, for a relationship between w and N ).Since there is only few competition on the available water vapour, these few ice crystals can grow to large sizes; thus aggregation is quite effective as we know from time scale analysis in sec.3.5.This can be seen in an aggregation factor approaching values at f aggr ∼ 0.01 − 0.1 -As soon as we proceed to colder temperatures and/or higher vertical velocities the picture changes.Under such conditions, much more ice crystals are produced in homogeneous nucleation events.Thus, although the ice crystals still live in a highly supersaturated environment, they have to compete for the available water vapour.Thus, they grow, but only to smaller sizes.Since for the efficiency of aggregation the size is much more important than the number concentration, in this regimes aggregation is less efficient.This leads to larger aggregation factors.For very cold temperatures, when ice crystals stay really at small sizes, aggregation is not important anymore, i.e. f aggr ∼ 0.9 − 1.
-The temperature dependence of the collision efficiency reduces the values but not the qualitative result.Thus, the aggregation factors are closer to 1 than in the case of E ≡ 1, however, the structure of the curves does not change

Lower initial temperatures
The box model simulations at 240 K have shown that aggregation generally has a weak effect on the evolution of crystal number concentrations.Since the collision efficiency decreases exponentially with decreasing temperature, we therefore expect a negligible effect of aggregation at lower temperatures, for all choices of f sed and both aggregation scenarios.This is indeed what we found from our box model simulations: the curves for the cases with and without aggregation do not differ drastically, thus we refrain from showing the simulations in detail.

Maximum impact of aggregation as derived from box model simulations
For determining the maximum impact of aggregation for different temperature and velocity regimes, we investigate the box model simulations with sedimentation factor f sed = 1, i.e. no sedimentation is allowed.In order to determine the maximum impact of aggregation we define an aggregation factor f agg : the factor is given by the ratio of ice crystal number concentrations of the aggregation run and the reference run (without aggregation): with t = end of simulation, and for temperature-dependent simulations the analogon is formed.In Fig. 7 the aggregation factor is shown for maximum aggregation (top panel) and temperature-dependent aggregation (bottom panel).
Figure 7 can be interpreted as follows: -For low vertical velocities at warm temperatures, homogeneous nucleation does produce just low number concentrations (see, e.g., Fig. 12 in Spichtinger and Gierens, 2009a, for a relationship between w and N).
Since there is only few competition on the available water vapour, these few ice crystals can grow to large sizes; thus aggregation is quite effective as we know from timescale analysis in Sect.3.5.This can be seen in an aggregation factor approaching values at f aggr ∼ 0.01-0.1 -As soon as we proceed to colder temperatures and/or higher vertical velocities the picture changes.Under such conditions, many more ice crystals are produced in homogeneous nucleation events.Thus, although the ice crystals still live in a highly supersaturated environment, they have to compete for the available water vapour.Thus, they grow, but only to smaller sizes.Since for the efficiency of aggregation the size is much more important than the number concentration, in this regimes aggregation is less efficient.This leads to larger aggregation factors.Fig. 7. Maximum impact of aggregation for box simulations using a constant updraught (i.e.cooling rate) at different initial temperatures.The impact of aggregation is determined by the aggregation factor defined in eq. ( 24).All simulations are without sedimentation in order to give the maximum impact.Top panel: Aggregation factor for different realistic vertical velocities (0.01 ≤ w ≤ 5 m s −1 ) at different inital temperatures.The collision efficiency is set to E ≡ 1. Bottom panel: same as in top panel, but for temperature dependent collision efficiency E = E(T ) as given by eq. ( 20).
Note, that the shape of the curves for different temperature regimes is nearly the same; actually, the curve of aggregation factor for temperature T init = 240K could be shifted to the left in order to represent the curves for other temperatures even quantitatively.

Test of aggregation within a 2D model: Simulations of synoptically driven cirrostratus
In order to investigate the impact of aggregation in a more realistic situation we implemented the new aggregation parameterization into the EULAG model including the already mentioned bulk microphysics scheme (Spichtinger and Gierens, 2009a).We investigate typical formation conditions for stratiform cirrus clouds, i.e. a synoptic scale updraught.
In the next subsection we present the setup of the simulations.Then we will present and discuss the results.Table 1.Initial vertical positions and temperature ranges for icesupersaturated layers (low/middle/high)

Setup
For simulating stratiform cirrus clouds as typical for midlatitudes we specify vertical profiles of temperature and pressure as shown in fig.8, respectively; additionally, we prescribe an ice-supersaturated layer with vertical extension of ∆z = 1.5km at different altitudes (top of layer at z top = 9/10/11km, i.e. low/middle/high).The vertical extension and the corresponding temperature ranges are presented in table 1.We use a 2D domain (x-z-plane) in the troposphere with a horizontal extension L x = 12.7km (∆x = 100m) and a vertical extension 4 ≤ z ≤ 14km (∆x = 50m).At initialisation the potential temperature field is superimposed by Gaussian noise with standard deviation σ θ = 0.025K.We choose a moderate wind shear for horizontal wind, i.e. du/dz = 10 −3 s −1 with u(z = 0) = 0 m s −1 , leading to a maximum wind of u max ≈ 10 m s −1 at z = 14km.The whole 2D domain is lifted with a constant vertical velocity.In order to investigate different synoptic conditions we choose two values w = 5cm s −1 and w = 8cm s −1 .In order to obtain similar conditions at the end of the simulations (i.e. the same vertical distance of lifting ∆z lif t = 720m or equivalently a cooling of ∆T ≈ 7.04K), the simulation time is adjusted; in case of w = 5cm s −1 the total simulation time is ∆t = 240min, whereas for w = 8 cm s −1 the total simulation Fig. 7. Maximum impact of aggregation for box simulations using a constant updraught (i.e.cooling rate) at different initial temperatures.The impact of aggregation is determined by the aggregation factor defined in Eq. ( 24).simulations are without sedimentation in order to give the maximum impact.Top panel: aggregation factor for different realistic vertical velocities (0.01 ≤ w ≤ 5 m s −1 ) at different initial temperatures.The collision efficiency is set to E ≡ 1. Bottom panel: same as in top panel, but for temperature dependent collision efficiency E = E(T ) as given by Eq. ( 20).
ice crystals stay really at small sizes, aggregation is not important anymore, i.e. f aggr ∼ 0.9-1.
-The temperature dependence of the collision efficiency reduces the values but not the qualitative result.Thus, the aggregation factors are closer to 1 than in the case of E ≡ 1, however, the structure of the curves does not change.
Note that the shape of the curves for different temperature regimes is nearly the same; actually, the curve of aggregation factor for temperature T init = 240 K could be shifted to the left in order to represent the curves for other temperatures even quantitatively.Fig. 7. Maximum impact of aggregation for box simulations using a constant updraught (i.e.cooling rate) at different initial temperatures.The impact of aggregation is determined by the aggregation factor defined in eq. ( 24).All simulations are without sedimentation in order to give the maximum impact.Top panel: Aggregation factor for different realistic vertical velocities (0.01 ≤ w ≤ 5 m s −1 ) at different inital temperatures.The collision efficiency is set to E ≡ 1. Bottom panel: same as in top panel, but for temperature dependent collision efficiency E = E(T ) as given by eq. ( 20).
Note, that the shape of the curves for different temperature regimes is nearly the same; actually, the curve of aggregation factor for temperature T init = 240K could be shifted to the left in order to represent the curves for other temperatures even quantitatively.

Test of aggregation within a 2D model: Simulations of synoptically driven cirrostratus
In order to investigate the impact of aggregation in a more realistic situation we implemented the new aggregation parameterization into the EULAG model including the already mentioned bulk microphysics scheme (Spichtinger and Gierens, 2009a).We investigate typical formation conditions for stratiform cirrus clouds, i.e. a synoptic scale updraught.
In the next subsection we present the setup of the simulations.Then we will present and discuss the results.Table 1.Initial vertical positions and temperature ranges for icesupersaturated layers (low/middle/high)

Setup
For simulating stratiform cirrus clouds as typical for midlatitudes we specify vertical profiles of temperature and pressure as shown in fig.8, respectively; additionally, we prescribe an ice-supersaturated layer with vertical extension of ∆z = 1.5km at different altitudes (top of layer at z top = 9/10/11km, i.e. low/middle/high).The vertical extension and the corresponding temperature ranges are presented in table 1.We use a 2D domain (x-z-plane) in the troposphere with a horizontal extension L x = 12.7km (∆x = 100m) and a vertical extension 4 ≤ z ≤ 14km (∆x = 50m).At initialisation the potential temperature field is superimposed by Gaussian noise with standard deviation σ θ = 0.025K.We choose a moderate wind shear for horizontal wind, i.e. du/dz = 10 −3 s −1 with u(z = 0) = 0 m s −1 , leading to a maximum wind of u max ≈ 10 m s −1 at z = 14km.The whole 2D domain is lifted with a constant vertical velocity.In order to investigate different synoptic conditions we choose two values w = 5cm s −1 and w = 8cm s −1 .In order to obtain similar conditions at the end of the simulations (i.e. the same vertical distance of lifting ∆z lif t = 720m or equivalently a cooling of ∆T ≈ 7.04K), the simulation time is adjusted; in case of w = 5cm s −1 the total simulation time is ∆t = 240min, whereas for w = 8 cm s −1 the total simulation Fig. 8. Initial vertical profiles for the simulations; left: temperature, middle: pressure, right: relative humidity with respect to ice for different set-ups (low, middle and high altitude range, corresponding to high, medium and low temperature range, see Table 1).

Test of aggregation within a 2-D model: simulations of synoptically driven cirrostratus
In order to investigate the impact of aggregation in a more realistic situation, we implemented the new aggregation parameterisation into the EULAG model including the already mentioned bulk microphysics scheme (Spichtinger and Gierens, 2009a).We investigate typical formation conditions for stratiform cirrus clouds, i.e. a synoptic scale updraught.In the next subsection we present the set-up of the simulations.Then we will present and discuss the results.

Set-up
For simulating stratiform cirrus clouds as typical for midlatitudes, we specify vertical profiles of temperature and pressure as shown in Fig. 8; additionally, we prescribe an icesupersaturated layer with vertical extension of z = 1.5 km at different altitudes (top of layer at z top = 9/10/11 km, i.e. low/middle/high).The vertical extension and the corresponding temperature ranges are presented in Table 1.We use a 2-D domain (x-z-plane) in the troposphere with a horizontal extension L x = 12.7 km ( x = 100 m) and a vertical extension 4 ≤ z ≤ 14 km ( x = 50 m).At initialisation the potential temperature field is superimposed by Gaussian noise with standard deviation σ θ = 0.025 K.We choose a moderate wind shear for horizontal wind, i.e. du/dz = 10 −3 s −1 with u(z = 0) = 0 m s −1 , leading to a maximum wind of u max ≈ 10 m s −1 at z = 14 km.a constant vertical velocity.In order to investigate different synoptic conditions we choose two values w = 5 cm s −1 and w = 8 cm s −1 .In order to obtain similar conditions at the end of the simulations (i.e. the same vertical distance of lifting z lift = 720 m or equivalently a cooling of T ≈ 7.04 K), the simulation time is adjusted; in case of w = 5 cm s −1 the total simulation time is t = 240 min, whereas for w = 8 cm s −1 the total simulation time is t = 150 min.As the results turned out to be similar for the chosen vertical velocities in terms of impact of aggregation, we will concentrate on a detailed investigation of the case w = 5 cm s −1 .For investigating the impact of aggregation we use three different setups: in the reference case, we switch off aggregation; in scenario "temperature-dependent" we use the full aggregation parameterisation including the temperature dependency, as described in Sect.4.1.2.In order to see the maximum effect of aggregation, we use the scenario "maximum impact", i.e. here the aggregation has efficiency E ≡ 1.We assume that ice forms by homogeneous nucleation only, parameterised after Koop et al. (2000).The background aerosol (sulfuric acid) is prescribed with a log-normal distribution with (dry) modal radius r m = 25 nm and geometrical standard deviation of σ r = 1.5.

Results and discussion
In general the simulations behave similar to those carried out by Spichtinger and Gierens (2009b) for the case of pure homogeneous nucleation: as the domain is lifted it cools by adiabatic expansion and the relative humidity increases until crystals are formed at the threshold for homogeneous nucleation.Figure 9 shows part of the temporal evolution of the reference simulation in time steps of 30 min, starting at t = 60 min, i.e. the state for t = 60/90/120/150 min is displayed.The formation of a quite homogeneous cirrostratus can be seen to occur after about 2 h simulation time by increasing ice water content IWC = q c • ρ air (black isolines).Some structure is formed by the horizontal wind driving small circulations inside the layer.In the further evolution, the vertical depth of the cloud layer is extended due to sedimenting ice crystals, reaching about z ∼ 3 km at the end of the simulation at t = 240 min.Figure 10 shows the results at the end of the simulations.Mean values of ice water content, ice crystal number concentration and relative humidity with respect to ice are shown, averaged over the domain.As expected, the impact of aggregation increases with temperature, even in the cases with E ≡ 1 where the aggregation efficiency itself has no temperature dependence.Obviously the remaining factors in the collision kernel contribute significantly to the temperature dependence and this is due to the fact that crystals can grow larger in the higher absolute humidity environment at higher temperatures.Thus, the geometrical factor evidently grows with temperature.The factor depending on the difference of terminal fall speeds increases on average with mean crystal size if the width of the size 12 E.Kienast-Sjögren et al.: Ice aggregation in 2-moment schemes time is ∆t = 150min.As the results turned out to be similar for the chosen vertical velocities in terms of impact of aggregation, we will concentrate on a detailed investigation of the case w = 5 cm s −1 .For investigating the impact of aggregation we use three different setups: In the reference case, we switch off aggregation; in scenario "temperature-dependent" we use the full aggregation parameterization including the temperature dependency, as described in sec.4.1.2.In order to see the maximum effect of aggregation, we use the scenario "maximum impact", i.e. here the aggregation has efficiency E ≡ 1.We assume that ice forms by homogeneous nucleation only, parameterized after Koop et al. (2000).The background aerosol (sulphuric acid) is prescribed with a lognormal distribution with (dry) modal radius r m = 25nm and geometrical standard deviation of σ r = 1.5, respectively.

Results and discussion
In general the simulations behave similar to those carried out by Spichtinger and Gierens (2009b) for the case of pure homogeneous nucleation: As the domain is lifted it cools by adiabatic expansion and the relative humidity increases until crystals are formed at the threshold for homogeneous nucleation.Fig. 9 shows part of the temporal evolution of the reference simulation in time steps of 30 min, starting at t = 60min, i.e. the state for t = 60/90/120/150min is displayed.The formation of a quite homogeneous cirrostratus can be seen to occur after about 2 h simulation time by increasing ice water content IWC = q c • ρ air (black isolines).Some structure is formed by the horizontal wind driving small circulations inside the layer.In the further evolution, the vertical depth of the cloud layer is extended due to sedimenting ice crystals, reaching about ∆z ∼ 3km at the end of the simulation at t = 240min.Fig. 10 shows the results at the end of the simulations.Mean values of ice water content, ice crystal number concentration and relative humidity wrt ice are shown, averaged over the domain.As expected, the impact of aggregation increases with temperature, even in the cases with E ≡ 1 where the aggregation efficiency itself has no temperature dependence.Obviously the remaining factors in the collision kernel contribute significantly to the temperature dependence and this is due to the fact that crystals can grow larger in the higher absolute humidity environment at higher temperatures.Thus, the geometrical factor evidently grows with temperature.The factor depending on the difference of terminal fall speeds increases on average with mean crystal size if the width of the size distribution does so.In the formulation of Spichtinger and Gierens (2009a) the width of the mass distribution (i.e. the square root of the second central moment) is proportional to the mean mass.Therefore higher temperatures lead to more aggregation also via the terminal velocities in this model.Aggregation increases the average mass of ice crystals and so leads to stronger sedimentation which has an effect on ice mass and number concentration in the simulated clouds.Mean crystal number concentrations and ice water contents are strongly (up to and partly exceeding a factor 2) reduced by aggregation in the simulation using the highest temperature.The effects are of similar quality in the colder cases, but smaller.Effects on the mean profiles of relative humidity are distribution does so.In the formulation of Spichtinger and Gierens (2009a) the width of the mass distribution (i.e. the square root of the second central moment) is proportional to the mean mass.Therefore higher temperatures lead to more aggregation also via the terminal velocities in this model.
Aggregation increases the average mass of ice crystals and therefore leads to stronger sedimentation which has an E. Kienast-Sjögren et al.: Ice aggregation in 2-moment schemes 13  effect on ice mass and number concentration in the simulated clouds.Mean crystal number concentrations and ice water contents are strongly (up to and partly exceeding a factor 2) reduced by aggregation in the simulation using the highest temperature.The effects are of similar quality in the colder cases, but smaller.Effects on the mean profiles of relative humidity are present, but here it is more instructive to look at the statistics (see below, Fig. 11).The pdfs of number concentration of ice crystals display in our simulations broad maxima at around N i ∼ 100 L −1 .Whereas there is hardly any effect on the statistics of number concentration in the low temperature case, aggregation shifts the peaks to lower values and broadens them.As expected, this effect becomes more prominent at higher temperatures.This is clearly a signature of the aggregation-enhanced sedimentation (see also discussion in Spichtinger and Gierens, 2009a).The high number concentration tails of the distributions are merely little affected by aggregation; this is quite plausible because high number concentrations are usually coupled with small ice crystals, thus aggregation is weakly effective in this range.
Also the total surface area of the ice crystals decreases by aggregation.This and the above mentioned increase of sedimentation fluxes diminish the sink for supersaturation and higher relative humidities are maintained over longer periods of time compared to the reference cases without aggregation.The statistics of relative humidities typically peak at values slightly above 100 %.This is best seen in the cold case where aggregation has hardly any effect on the pdf of RH i .At the higher temperatures, where aggregation becomes more efficient we see the peaks shifted to higher values (i.e. more significant "quasi equilibrium" supersaturation), clearly an effect of the less effective sink for water vapour in a cloud affected by aggregation (see also the discussion in Spichtinger and Cziczo, 2010).This effect is most pronounced in those regions of a cloud that otherwise approach ice saturation most quickly, typically the middle part of the cloud.Thus aggregation contributes to ice supersaturation within relatively warm cirrus clouds.Cold cirrus is hardly affected by aggregation according to our simulations (and under the condition that the gravitational collection kernel is the only relevant one for cirrus clouds).
In the simulations with stronger updraught we can see differences in details, but qualitatively they behave similar to the simulations shown.Therefore, we do not deem it necessary to present them here.

Scaling size distribution
The log-normal crystal mass distribution is used in the scheme of Spichtinger and Gierens (2009a) and the question arises whether this is an appropriate choice in situations where aggregation dominates ice growth.In such cases the evolving crystal size distribution apparently can be mapped onto a universal shape: , where m 0 is a mass unit (to make the prefactor dimensionless) and ψ(x) is a function that depends on time only via the ratio m/m m (t) (Field and Heymsfield, 2003;Westbrook et al., 2004Westbrook et al., , 2007;;Sölch and Kärcher, 2011).The temporal evolution of the size distribution can thus be captured by scaling both axes with a function of the time-varying modal mass, the m axis with its inverse and the f axis with its θ-th power.Although the log-normal distribution can be treated in this way (with θ = 1 for constant geometric width), it is not a true solution for gravitational aggregation.Researchers, guided by numerical simulations, tend to use for ψ(x) functions with an exponential upper tail (e.g.Field and Heymsfield, 2003, use a gamma distribution), however, it is not mathematically proven that the gravitational collection kernel has a scaling solution at all (Aldous, 1999).One condition for this is the homogeneity of the kernel function K(m, M).
From the mass-length and length-fall speed relations we can see that the kernel can only be a homogeneous function if the exponents in these relations are constant over all sizes.This is only so for spheres.The exponents change, however, with size for ice crystals and the shape of the crystals, for instance expressed by the aspect ratio, changes with size (see, for instance, Heymsfield and Iaquinta, 2000).In our case with hexagonal columns it is the function A(m), the surface of the crystal, where m appears with two different exponents (for the basal and prism faces, respectively); A(m) is nonhomogeneous.Hence, for ice the collection kernel is not a homogeneous function and therefore there is no true scaling.
If it is nevertheless possible to match observed size distributions to a universal scaling function, this suggests that there is an approximate scaling for large ( 100 µm) crystals, which requires that the mentioned exponents are mainly constant in this range.The upper tail of a log-normal distribution is slightly bent upward on a semilog plot, that is, it does not have a perfect exponential; but it might be possible to fit the observed size distributions of large crystals with log-normals as well as with gamma-distributions (in particular, given the noise in the data).Unless there is a mathematical proof that ice aggregation leads to a specific size distribution other than the log-normal there is no need to abandon it.

Conclusions
We have derived from the master-equation for coagulation a simple formulation of aggregation for two-moment bulk microphysical models.So far we developed the formulation for aggregation of crystals belonging to the same class only (the microphysics scheme of Spichtinger and Gierens (2009a) allows more than one class of ice).A more general formulation with aggregation of crystals from different classes is rather a numerical than a mathematical problem, but a difficult one; thus, it is beyond the scope of our study.The core of the present formulation is a double integral of the collection kernel weighted with the crystal size (or mass) distribution, which is the expectation value of the kernel.This quantity is to be inserted into the differential equation for the crystal number concentration which is of a form that was already derived by Smoluchowski (1916Smoluchowski ( , 1917)).The double integrals are evaluated numerically for log-normal size distributions over a large range of geometric mean masses.The direct evaluation of the integrals within a cloud simulation run takes a lot of computing time and is not recommended.Instead the pre-calculated results can either be read from a look-up table or -even better -a polynomial fit of the results can be used that yields good accuracy.
We have tested the new parameterisation in various environments: stand-alone (to see how the solution of the differential equation behaves and to test the polynomial fits), in a box-model (where aggregation occurs simultaneously with other microphysical processes and where a first check can be made, whether and when aggregation is important), and in a 2-D simulation of a cirrostratus cloud (where additionally, cloud dynamics can enhance or dampen the effects of aggregation).Overall these tests suggest that aggregation can become important at (relatively) warmer cirrus temperatures, affecting not only ice number and mass concentrations, but leading also to higher and longer-lasting in-cloud supersaturation.Sedimentation fluxes are increased when aggregation is switched on.Cold cirrus clouds are hardly affected by aggregation.The temperature dependence originates not only from an assumed temperature dependence of the collection efficiency but also from the other factors in the collision kernel: Higher temperatures imply larger ice crystals and larger spread in terminal velocities (if the assumed type of size distribution is such that the width of it increases with increasing mean size).From timescale analysis the importance of aggregation can be derived depending on number concentration and size of the ice crystals.For cold clouds it is often justified to ignore aggregation when the research focus is on ice mass and number densities.However, when the focus of research is crystal habits and their effect on radiation, aggregation should not be ignored since cirrus clouds usually contain complex, irregular and imperfect ice crystals as reported by Bailey and Hallett (2009).The authors have shown that even cold cirrus contains complex ice crystals that may be the result of aggregation.They occur predominantly at high supersaturation, but supersaturation does not directly appear in the formulation of the kernel function.One might test whether an extension of the formulation of the collection efficiency (i.e.E = E(T , S)) yields better results; however, this kind of further development as well as sensitivity studies is far beyond the scope of this study and is left for future work.Rather it is desirable to measure collection efficiencies in big cloud chambers.By doing so, one should simultaneously test whether aggregation is the only process that leads to complex forms of ice crystals.This is not probable since rosette shaped crystals are too regular to be formed by random collisions.There is obviously a gap in our understanding and much occasion remains for further experimental research before we should develop our numerical formulations into unjustified detail.

Fig. 1 .
Fig. 1.K integral values (in m 3 s −1 ) as a function of the modal mass (in kg) of the crystal mass distribution and for σ m = 2.85.The corresponding sizes range between a modal length of 0.6 and 8000 µm.Exact values from a numerical integration are shown as squares.The solid lines represent polynomial fits in 4 different mass ranges, indicated by the black lines.

Fig. 1 .
Fig. 1.K integral values (in m 3 s −1 ) as a function of the modal mass (in kg) of the crystal mass distribution and for σ m = 2.85.The corresponding sizes range between a modal length of 0.6 and 8000 µm.Exact values from a numerical integration are shown as squares.The solid lines represent polynomial fits in 4 different mass ranges, indicated by the black lines.

Fig. 2 .
Fig. 2. Fits on calculated aggregation kernels for different geometric standard deviations.Upper panel: Aggregation kernels vs. modal mass.The change due to the width of the distribution can be seen clearly.Lower panel: Aggregation kernels vs. mean mass.Since the mean mass depends on the geometric width of the distribution, the variation of the kernels due to different width is smaller.Additionally, the aggregation kernel as derived bySchumann (2012) is shown for comparison.The corresponding sizes range between a modal length of 0.6 and 8000 µm.

Fig. 2 .
Fig.2.Fits on calculated aggregation kernels for different geometric standard deviations.Upper panel: aggregation kernels vs. modal mass.The change due to the width of the distribution can be seen clearly.Lower panel: aggregation kernels vs. mean mass.Since the mean mass depends on the geometric width of the distribution, the variation of the kernels due to different width is smaller.Additionally, the aggregation kernel as derived bySchumann (2012) is shown for comparison.The corresponding sizes range between a modal length of 0.6 and 8000 µm.

Fig. 4 .
Fig.4.Particle number density N0 at the start of the iteration (x-Axis) and after 1000 s of iteration (N , y-Axis).For small modal masses, the particle number density hardly changes during the iteration and the plot is almost a straight line.As reference for no aggregation, a black line is plotted.As expected, larger modal masses result in increased aggregation.Thus the particle number density decreases during integration, which is shown in the graph.The corresponding size to the modal masses plotted ranges between 6 and 345 µm.Note that the tests have been performed with E ≡ 1, that is, the maximum effect of aggregation is seen here.The red line is almost equal to the black line, thus can hardly be seen in the plot.

Fig. 5 .
Fig.5.Particle number density N0 at the beginning of the iteration (x-Axis) and after 1000 s of iteration (N, y-Axis) for a modal mass of 10 −9 kg, which corresponds to a particle with modal length 345 µm.As a reference for no aggregation, a black line is plotted.Adding a temperature dependency shows a weakened aggregation effect with decreasing temperatures.The lower the tempertures the more particles are present at the end of the simulation.

Fig. 4 .
Fig.4.Particle number density N 0 at the start of the iteration (x-Axis) and after 1000 s of iteration (N , y-Axis).For small modal masses, the particle number density hardly changes during the iteration and the plot is almost a straight line.As reference for no aggregation, a black line is plotted.As expected, larger modal masses result in increased aggregation.Thus the particle number density decreases during integration, which is shown in the graph.The corresponding size to the modal masses plotted ranges between 6 and 345 µm.Note that the tests have been performed with E ≡ 1, that is, the maximum effect of aggregation is seen here.The red line is almost equal to the black line, thus can hardly be seen in the plot.
Fig.5.Particle number density N 0 at the beginning of the iteration (x-Axis) and after 1000 s of iteration (N, y-Axis) for a modal mass of 10 −9 kg, which corresponds to a particle with modal length 345 µm.As a reference for no aggregation, a black line is plotted.Adding a temperature dependency shows a weakened aggregation effect with decreasing temperatures.The lower the tempertures the more particles are present at the end of the simulation.

Fig. 6 .
Fig. 6.Ice crystal number concentration as a function of time for different sedimentation factors.Top panel: f sed = 1.0, i.e. no sedimentation, leading to the strongest impact of aggregation; middle panel: f sed = 0.5, representing the upper part of a cloud; bottom panel:f sed = 0.9, representing the lower part of a cloud.Simulations without aggregation are indicated by red lines, simulations including aggregation with maximum strength are represented by blue lines and simulations with temperature-dependent aggregation are indicated by gree lines.All simulations starts at T = 240K, p = 300hPa and are driven by a constant vertical updraught of w = 0.05 m s −1 .

Fig. 6 .
Fig. 6.Ice crystal number concentration as a function of time for different sedimentation factors.Top panel: f sed = 1.0, i.e. no sedimentation, leading to the strongest impact of aggregation; middle panel: f sed = 0.5, representing the upper part of a cloud; bottom panel: f sed = 0.9, representing the lower part of a cloud.Simulations without aggregation are indicated by red lines, simulations including aggregation with maximum strength are represented by blue lines and simulations with temperature-dependent aggregation are indicated by green lines.All simulations starts at T = 240 K, p = 300 hPa and are driven by a constant vertical updraught of w = 0.05 m s −1 .
For very cold temperatures, when E.Kienast-Sjögren et al.: Ice aggregation in 2

Fig. 8 .
Fig. 8. Initial vertical profiles for the simulations; left: temperature, middle: pressure, right: relative humidity wrt ice for different setups (low, middle and high altitude range, corresponding to high, medium and low temperature range, see tab. 1).

Fig. 8 .
Fig. 8. Initial vertical profiles for the simulations; left: temperature, middle: pressure, right: relative humidity wrt ice for different setups (low, middle and high altitude range, corresponding to high, medium and low temperature range, see tab. 1).

Fig. 10 .Fig. 11 .
Fig. 10.Vertical profiles of mean ice crystal number concentration (left), ice water content (middle) and relative humidity wrt ice (right) at the end of the simulation (t = 240min, w = 5 cm s −1 ) for different temperature regimes.Top row: low temperature conditions, middle row: medium temperature conditions, bottom row: high temperature conditions

Fig. 10 .
Fig. 10.Vertical profiles of mean ice crystal number concentration (left), ice water content (middle) and relative humidity with respect to ice (right) at the end of the simulation (t = 240 min, w = 5 cm s −1 ) for different temperature regimes.Top row: low temperature conditions, middle row: medium temperature conditions, bottom row: high temperature conditions.

Table 1 .
Initial vertical positions and temperature ranges for icesupersaturated layers (low/middle/high).

www.atmos-chem-phys.net/13/9021/2013/ Atmos. Chem. Phys., 13, 9021-9037, 2013 9034 E. Kienast-Sjögren et al.: Ice aggregation in two-moment schemes
Fig. 10.Vertical profiles of mean ice crystal number concentration (left), ice water content (middle) and relative humidity wrt ice (right) at the end of the simulation (t = 240min, w = 5 cm s −1 ) for different temperature regimes.Top row: low temperature conditions, middle row: medium temperature conditions, bottom row: high temperature conditions Statistics of ice crystal number concentrations (top row) and relative humidity wrt ice (bottom row) for different temperature regimes (left: low temperature conditions, middle: medium temperature conditions, right: high temperature conditions).Statistics of ice crystal number (top row) and relative humidity with respect to ice (bottom row) for different temperature regimes (left: low temperature conditions, middle: medium temperature conditions, right: high temperature conditions).

Table A1 .
Coefficients a k , k = 1, . . ., 4 for the fitting polynomials as given by Eq. (A1) for different geometric standard deviations σ m of the underlying ice crystal mass distribution of log-normal type.