Analytical treatment of ice sublimation and test of sublimation parameterisations in two–moment ice microphysics models

We derive an analytic solution to the spectral growth/sublimation equation for ice crystals and apply it to idealised cases. The results are used to test parameterisations of the ice sublimation process in two–moment bulk microphysics models. Although it turns out that the relation between number loss fraction and mass loss fraction is not a function since it is not unique, it seems that a functional parameterisation is the best that one can do in a bulk model. Testing a more realistic case with humidity oscillations shows that artificial crystal loss can occur in simulations of mature cirrus clouds with relative humidity fluctuating about ice saturation.


Introduction
Cloud microphysical models can essentially be grouped into bulk and bin models.Bulk models describe a cloud in terms of gross quantities like total water or ice mass and total number concentration of hydrometeors.Although an assumption of the underlying size or mass spectrum of the hydrometeors is usually made, its evolution is not completely described since only the low-order moments are predicted in the model equations.Bin models instead describe a cloud in much more detail by explicitly accounting for the size or mass spectrum.
Here the mass spectrum is divided into a large number of bins and the microphysical processes directly change the number of particles in the various bins (some recent examples can be found in Lin et al., 2002;Monier et al., 2006;Leroy et al., 2007Leroy et al., , 2009)).While bulk microphysics models are designed rather for computational speed than for an exact description of the evolution of a cloud's mass spectrum, it is vice versa for the bin models which pay for more exact results with Correspondence to: K. Gierens (klaus.gierens@dlr.de)higher computational costs, that is, longer computing time and higher memory requirements.For many applications (for instance in an operational weather forecast model) it is too expensive to use a bin model, hence bulk models are needed.Nonetheless, their behaviour should be tested under a large variety of situations such that their strenghts and weaknesses become understood as much as possible.
Cloud microphysical models that use a two-moment bulk scheme predict both the temporal and spatial variations of the total number and total mass of the droplets and ice crystals.The temporal change of number and mass concentrations are determined by process rates which are sometimes difficult to formulate.One example is deposition.When ice crystals grow by vapour deposition, the ice mass grows accordingly while the number of crystals is constant.This is the easy case.The difficult case arises in subsaturated conditions.When ice crystals sublimate, the ice mass decreases accordingly, while the decrease rate of crystal number cannot be formulated in a straightforward way.The prognostic equations for ice number and mass concentration (for recent examples cf.Morrison and Grabowski, 2008;Spichtinger and Gierens, 2009) contain deposition terms, DEP (for mass concentration) and N DEP (for number concentration) in the nomenclature of Spichtinger and Gierens (2009).
The problem is now, that there is an equation for DEP , but not for N DEP .Harrington et al. (1995) have performed a large series of simulations to find an appropriate parameterisation for N DEP .Spichtinger and Gierens (2009) roughly followed their results and used the following simple approximation to calculate N DEP : where f m is the mass fraction sublimated in the current time step, f n is the desired number fraction, and α is a constant parameter.Harrington et al. (1995) suggest a range of 1≤α≤1.5 for this parameter, which is set to α = 1.1 by Published by Copernicus Publications on behalf of the European Geosciences Union.

7482
K. Gierens and S. Bretl: Analytical treatment of sublimation Spichtinger and Gierens (2009).Morrison and Grabowski (2008) use α = 1, i.e. f n = f m .The constancy of α is certainly an oversimplification, since this would mean that ice crystals in fall streaks sublimate in a similar way than ice crystals deep inside a cloud where small humidity fluctuations around ice saturation might lead to crystal loss.Hence, we deemed it worth to investigate this process in more detail and to test potential alternative parameterisations.
For this purpose we will use an analytic solution to the problem, proceeding from the spectral form of the growth/sublimation equation, to be derived next in Sect. 2. From the analytical solution we compute number and mass loss fractions (Sect.3) and derive timescales (Sect.4).Alternative parameterisations will be tested in Sect. 5. A discussion and conclusions are presented in the final Sects.6 and 7.

Time dependent crystal mass distribution
As ice crystals in a cloud generally have a large variety of sizes, bulk microphysics models take account of this variety by implicitly or explicitly assuming an underlying crystal size distribution, which can be expressed as a crystal mass distribution or mass spectrum, f (m).Crystal growth by vapour deposition or crystal sublimation changes the mass spectrum, so that it becomes also a function of time, i.e. f (m, t).The temporal change of the mass spectrum due to deposition/sublimation can be described by an equation that has the form of a continuity equation in mass space (see, e.g., Wacker and Herbert, 1983): f (m, t) dm can equivalently be interpreted as the time dependent probability that an ice crystal selected randomly in a cloud has a mass between m and m + dm.The probabilistic interpretation will turn out useful for the analytical derivations presented in the following.(f (m, t) has to be distinguished carefully from fractions f k which are marked by a lower index.)ṁ ( = dm/dt) is the depositional mass growth/sublimation rate of a single crystal of mass m.Koenig (1971) has shown that ṁ can be parameterized as a function of m in the simple form: Here, a and b depend on temperature and pressure, and a depends additionally on the degree of ice supersaturation.For subsaturated cases, a<0.Spichtinger and Gierens (2009, their Fig. 5) show that this approximation is very good over four orders of magnitude in m.U M and U T are unit mass and unit time, for which we use 1 ng and 1 s, respectively.In the following equations we drop the reference to the units, but expressions like m b are always meant as (m/U M) b .With the chosen units, a certain value of a (which is a number) means that the sublimation rate of a 1 ng ice crystal is a ng/s.In this paper we use crystal masses and values of a that are typical in cirrus clouds under subsaturated conditions, i.e. the temperature is below −38 • C, the pressure below 300 hPa and the relative humidity with respect to ice ranges from say 10 to 98%.We used b = 0.37 and a width of the mass distribution of σ m = 2.85.Note that the exact values are not important for the principle considerations that follow.
We introduce the following coordinate transformation: That is, in the new coordinates the growth/sublimation rate is a simple constant.Let g(x) be the probability density function of x.Then the spectral growth equation obtains the simple form: which is solved by a characteristic: g can be any function in principle, but here we use probability density functions, of course.It is obvious that the time development of g(x, t) is merely a shift of the function as a whole along the x-axis.All central moments of order higher than one are invariant, that is, the shape of g(x, t) does not change over time.
We start (t = 0) from a certain initial mass distribution (e.g., log-normal), normalized to one for the present paper: substitute m with x (this yields another log-normal with and then replace x with x −at, which represents the temporal evolution of the mass distribution in the x-coordinates: Singularities of this formulation can appear at negative m which are irrelevant.The singularity developing at m = 0 is integrable.The temporal evolution of the initially lognormal distribution is shown in Fig. 1 for a sublimation condition (a < 0).It is evident that in the mass coordinates the distribution does not retain its initial shape; deviations from the initial shape are very pronounced.When we consider instead a growth situation (a > 0) the initial shape is much better conserved (not shown).A more general solution of the spectral growth equation is presented in Appendix A.

Mass loss and number loss fractions
We assume that ice crystals have a minimum mass of m thr (e.g. 10 −3 ng).Then we define for k ∈ {0, 1} the following integrals which give the total number and mass fractions of the ice mass distribution exceeding the threshold: Since f depends on time, so does I k .Harrington et al. (1995) consider the fractional total mass loss (φ 1 ) and number loss (φ 0 ) at time t which can be written as Note that we use the notation with numerical indices (since the expressions using moments suggest it) interchangeably with the notation using (m, n) as indices (the notation that is used in the quoted papers), i.e. the following identities hold: Time runs along the curves from (0, 0) to (1, 1).We see that initially there is mass loss combined with negligible number loss, but the number loss overhauls in the later phases of the ongoing sublimation process.This is what we expect for mass distributions with mode masses exceeding the threshold mass.Different behaviour can be expected for exponential distributions.We note that Harrington et al. (1995) also had examples with a different behaviour.
In a cloud model, we generally do neither have knowledge of the initial values I k (0) nor of the time t which must be interpreted here as the time passed by since sublimation started.(Since many processes act in a cloud model simultaneously it would not even make sense to introduce such a time variable or to track the "initial" values).Hence, in a model we usually can only consider the fractional mass and number loss per time step, which can be written as These fractions depend on the size of the time step ∆t, but still on the time since sublimation started.
Singularities of this formulation can appear at negative m which are irrelevant.The singularity developing at m = 0 is integrable.The temporal evolution of the initially log-normal distribution is shown in Fig. 1 for a sublimation condition (a<0).It is evident that in the mass coordinates the distribution does not retain its initial shape; deviations from the initial shape are very pronounced.When we consider instead a growth situation (a>0) the initial shape is much better conserved (not shown).A more general solution of the spectral growth equation is presented in Appendix A.

Mass loss and number loss fractions
We assume that ice crystals have a minimum mass of m thr (e.g. 10 −3 ng).Then we define for k ∈ {0, 1} the following integrals which give the total number and mass fractions of the ice mass distribution exceeding the threshold: Singularities of this formulation can appear at negative m which are irrelevant.The singularity developing at m = 0 is integrable.The temporal evolution of the initially lognormal distribution is shown in Fig. 1 for a sublimation condition (a < 0).It is evident that in the mass coordinates the distribution does not retain its initial shape; deviations from the initial shape are very pronounced.When we consider instead a growth situation (a > 0) the initial shape is much better conserved (not shown).A more general solution of the spectral growth equation is presented in Appendix A.

Mass loss and number loss fractions
We assume that ice crystals have a minimum mass of m thr (e.g. 10 −3 ng).Then we define for k ∈ {0, 1} the following integrals which give the total number and mass fractions of the ice mass distribution exceeding the threshold: Since f depends on time, so does I k .Harrington et al. (1995) consider the fractional total mass loss (φ 1 ) and number loss (φ 0 ) at time t which can be written as Note that we use the notation with numerical indices (since the expressions using moments suggest it) interchangeably with the notation using (m, n) as indices (the notation that is used in the quoted papers), i.e. the following identities hold: φ 0 ≡ φ n , φ 1 ≡ φ m , I 0 ≡ I n , I 1 ≡ I m (and correspondingly for the f k introduced below).φ n (φ m ) is plotted in Fig. 2. Time runs along the curves from (0, 0) to (1, 1).We see that initially there is mass loss combined with negligible number loss, but the number loss overhauls in the later phases of the ongoing sublimation process.This is what we expect for mass distributions with mode masses exceeding the threshold mass.Different behaviour can be expected for exponential distributions.We note that Harrington et al. (1995) also had examples with a different behaviour.
In a cloud model, we generally do neither have knowledge of the initial values I k (0) nor of the time t which must be interpreted here as the time passed by since sublimation started.(Since many processes act in a cloud model simultaneously it would not even make sense to introduce such a time variable or to track the "initial" values).Hence, in a model we usually can only consider the fractional mass and number loss per time step, which can be written as These fractions depend on the size of the time step ∆t, but still on the time since sublimation started.Examples for various time step values ∆t are plotted in Fig. 3: These curves link pairs {f m (t i , ∆t), f n (t i , ∆t)} for times t i = i • ∆t (i = 1, 2, 3, . ..) running until the ice is completely sublimated.(The yellow dot at (1, 1) was computed with a very large time step, such that the ice sublimated completely within this single time step).The curves f n (f m ) all have similar shape.They start near the f m -axis which means that the initial mass loss occurs alone almost without number loss, Since f depends on time, so does I k .Harrington et al. (1995) consider the fractional total mass loss (φ 1 ) and number loss (φ 0 ) at time t which can be written as Note that we use the notation with numerical indices (since the expressions using moments suggest it) interchangeably with the notation using (m, n) as indices (the notation that is used in the quoted papers), i.e. the following identities hold: φ 0 ≡φ n , φ 1 ≡φ m , I 0 ≡I n , I 1 ≡I m (and correspondingly for the f k introduced below).φ n (φ m ) is plotted in Fig. 2. Time runs along the curves from (0, 0) to (1, 1).We see that initially there is mass loss combined with negligible number loss, but the number loss overhauls in the later phases of the ongoing sublimation process.This is what we expect for mass distributions with mode masses exceeding the threshold mass.Different behaviour can be expected for exponential distributions.We note that Harrington et al. (1995) also had examples with a different behaviour.
In a cloud model, we generally do neither have knowledge of the initial values I k (0) nor of the time t which must be interpreted here as the time passed by since sublimation started.(Since many processes act in a cloud model simultaneously it would not even make sense to introduce such a time variable or to track the "initial" values).Hence, in a model we usually can only consider the fractional mass and number loss per time step, which can be written as The yellow point at (1, 1) was the result of a test with a very long time step that should guarantee that the ice sublimates completely within that step.The black line is the attempt of a fit of the later parts of the curves, fn=f 0.89 m .The initial geometric mean mass for the calculations is m 0 =1 ng. and after about 600 s they reach an "attractor" that can be approximately be fitted by f n = f α m with α = 0.89.α < 1 again signifies that later in the sublimation process the number loss catches up.The shape of the curves shows that f n is not really a function (in the mathematical sense) of f m because it is not unique even if the timestep is fixed.This makes it questionable whether a functional dependence for parameterisation of sublimation should be used in bulk models.
The timestep dependence of f k suggests a Taylor expansion which yields As Fig. 3 shows, the higher orders cannot be neglected at time steps of the order 100 s and longer.Hence the following analysis is strictly valid only for cloud resolving models with small time steps of the order seconds (where the higher order terms are negligible), but it might still give some guidance for the treatment of sublimation in large-scale models.First, we have which still retains the time dependence.We note that 1/f k (t) can be interpreted as a timescale for the change of the cor- responding integral.The derivatives of the integrals can be computed in the following way: and the partial derivative of the mass distribution function is: f n vs. f m is plotted in Fig. 4.These curves look similar to those in Fig. 3, and indeed they are equivalent to those for a unit time step of 1 s.

Time scales
The curves φ n (φ m ) in Fig. 2 look very similar for different initial mean masses, but if we would include tick marks for time along the curves they would differ for different curves.We can therefore try to unify these functions with respect to time as well by introducing a dimensionless time τ .The corresponding τ -tick marks would then become nearly congruent.For this we can first compute the time T it needs to sublimate completely (i.e. to mass zero) a crystal having exactly the mass m 0 , which can be any characteristic mass of the chosen distribution, (for instance the median, mean, or Fig. 3. Mass loss (f m ) and number loss (f n ) fractions for series of subsequent sublimation time steps.Time steps are 100, 200, 300, 400 s (solid red, green, blue, dashed red), and time runs along each curve from the point neares to the f m -axis to the point nearest to the f n -axis.The individual time steps are marked on each curve.Arrows indicate the direction of time.The yellow point at (1, 1) was the result of a test with a very long time step that should guarantee that the ice sublimates completely within that step.The black line is the attempt of a fit of the later parts of the curves, f n = f 0.89 m .The initial geometric mean mass for the calculations is m 0 = 1 ng.
These fractions depend on the size of the time step t, but still on the time since sublimation started.Examples for various time step values t are plotted in Fig. 3: These curves link pairs {f m (t i , t), f n (t i , t)} for times t i = i • t (i = 1, 2, 3, . ..) running until the ice is completely sublimated.(The yellow dot at (1, 1) was computed with a very large time step, such that the ice sublimated completely within this single time step).The curves f n (f m ) all have similar shape.They start near the f m -axis which means that the initial mass loss occurs alone almost without number loss, and after about 600 s they reach an "attractor" that can be approximately be fitted by f n = f α m with α = 0.89.α < 1 again signifies that later in the sublimation process the number loss catches up.The shape of the curves shows that f n is not really a function (in the mathematical sense) of f m because it is not unique even if the timestep is fixed.This makes it questionable whether a functional dependence for parameterisation of sublimation should be used in bulk models.
The timestep dependence of f k suggests a Taylor expansion which yields As Fig. 3 shows, the higher orders cannot be neglected at time steps of the order 100 s and longer.Hence the following The yellow point at (1, 1) was the result of a test with a very long time step that should guarantee that the ice sublimates completely within that step.The black line is the attempt of a fit of the later parts of the curves, fn=f 0.89 m .The initial geometric mean mass for the calculations is m0=1 ng. and after about 600 s they reach an "attractor" that can be approximately be fitted by f n = f α m with α = 0.89.α < 1 again signifies that later in the sublimation process the number loss catches up.The shape of the curves shows that f n is not really a function (in the mathematical sense) of f m because it is not unique even if the timestep is fixed.This makes it questionable whether a functional dependence for parameterisation of sublimation should be used in bulk models.
The timestep dependence of f k suggests a Taylor expansion which yields As Fig. 3 shows, the higher orders cannot be neglected at time steps of the order 100 s and longer.Hence the following analysis is strictly valid only for cloud resolving models with small time steps of the order seconds (where the higher order terms are negligible), but it might still give some guidance for the treatment of sublimation in large-scale models.First, we have which still retains the time dependence.We note that 1/f k (t) can be interpreted as a timescale for the change of the cor- responding integral.The derivatives of the integrals can be computed in the following way: and the partial derivative of the mass distribution function is: f n vs. f m is plotted in Fig. 4.These curves look similar to those in Fig. 3, and indeed they are equivalent to those for a unit time step of 1 s.

Time scales
The curves φ n (φ m ) in Fig. 2 look very similar for different initial mean masses, but if we would include tick marks for time along the curves they would differ for different curves.We can therefore try to unify these functions with respect to time as well by introducing a dimensionless time τ .The corresponding τ -tick marks would then become nearly congruent.For this we can first compute the time T it needs to sublimate completely (i.e. to mass zero) a crystal having exactly the mass m 0 , which can be any characteristic mass of the chosen distribution, (for instance the median, mean, or analysis is strictly valid only for cloud resolving models with small time steps of the order seconds (where the higher order terms are negligible), but it might still give some guidance for the treatment of sublimation in large-scale models.
First, we have which still retains the time dependence.We note that 1/f k (t) can be interpreted as a timescale for the change of the corresponding integral.The derivatives of the integrals can be computed in the following way: and the partial derivative of the mass distribution function is: 4.These curves look similar to those in Fig. 3, and indeed they are equivalent to those for a unit time step of 1 s.

Time scales
The curves φ n (φ m ) in Fig. 2 look very similar for different initial mean masses, but if we would include tick marks for time along the curves they would differ for different curves.We can therefore try to unify these functions with respect to time as well by introducing a dimensionless time τ .The corresponding τ -tick marks would then become nearly congruent.For this we can first compute the time T it needs to sublimate completely (i.e. to mass zero) a crystal having exactly the mass m 0 , which can be any characteristic mass of the chosen distribution, (for instance the median, mean, or mode mass).Here we take for m 0 the geometric mean mass of the chosen log-normal distribution.This time is With this time we introduce τ as t/T , i.e.
Hence, the dimensionless time variable takes into account the chosen initial characteristic mass of the mass distribution and the current sublimation rate.It does not take into account other quantities like the width of the mass distribution.We define h(τ ) can be interpreted as the ratio of two timescales for change of the I k integrals.These timescales and their ratio vary in time.h(τ ) is plotted in Fig. 5 for various initial geometric mean masses (from 1 to 1000 ng) and various sublimation rates.We can observe the following behaviour: up to a normalised time τ ≈0.1 sublimation leads almost only to a mass loss and the number density of the ice crystals stays nearly constant.Number loss commences at about τ ≈0.1 and the fractional number loss rate exceeds the fractional mass loss rate (h>1) at about τ ≈1 and thereafter, that is, once the characteristic mass had sufficient time to sublimate completely.The shape of a curve h(τ ) depends on the width (and shape) of the initial mass spectrum.The transition from mass loss to number loss regime is sudden for a narrow distribution and gradual for a broad one (not shown).However, as long as we chose a fixed value for σ m (as in the figure), all these functions are very similar and they can be fitted with a generalised log-logistic function: The formulation is such that τ 0 is the inflexion point of the fit function.The latter is plotted in Fig. 5 as well with ψ = 1.24, p = 2.4, τ 0 = 0.3.ψ is the asymptotic value of the fit for large values of τ .It exceeds unity which is an expression for our earlier finding that in the end the number loss exceeds the mass loss rate.The parameter p controls the steepness of the fit around the inflexion point, i.e. it measures K. Gierens and S. Bretl: Analytical treatment of sublimation 5 mode mass).Here we take for m 0 the geometric mean mass of the chosen log-normal distribution.This time is With this time we introduce τ as t/T , i.e.
Hence, the dimensionless time variable takes into account the chosen initial characteristic mass of the mass distribution and the current sublimation rate.It does not take into account other quantities like the width of the mass distribution.We define h(τ ) can be interpreted as the ratio of two timescales for change of the I k integrals.These timescales and their ratio vary in time.h(τ ) is plotted in Fig. 5 for various initial geometric mean masses (from 1 to 1000 ng) and various sublimation rates.We can observe the following behaviour: up to a normalised time τ ≈ 0.1 sublimation leads almost only to a mass loss and the number density of the ice crystals stays nearly constant.Number loss commences at about τ ≈ 0.1 and the fractional number loss rate exceeds the fractional mass loss rate (h > 1) at about τ ≈ 1 and thereafter, that is, once the characteristic mass had sufficient time to sublimate completely.The shape of a curve h(τ ) depends on the width (and shape) of the initial mass spectrum.The transition from mass loss to number loss regime is sudden for a narrow distribution and gradual for a broad one (not shown).However, as long as we chose a fixed value for σ m (as in the figure), all these functions are very similar and they can be fitted with a generalised log-logistic function: The formulation is such that τ 0 is the inflexion point of the fit function.The is plotted in Fig. 5 as well with ψ = 1.24, p = 2.4, τ 0 = 0.3.ψ is the asymptotic value of the fit for large values of τ .It exceeds unity which is an expression for our earlier finding that in the end the number loss exceeds the mass loss rate.The parameter p controls the steepness of the fit around the inflexion point, i.e. it measures how fast sublimation changes from the mass loss to the number loss regime.Finally, the inflexion point is found here at 0.3 which means that the transition into the number loss regime occurs at a considerably shorter time than is needed to sublimate a crystal with the chosen characteristic mass.We see that it is possible to find a certain universal behaviour of the sublimation curves that do not depend on initial mean mass or sublimation rate.The function does however change when we use different initial mass distributions, e.g. by changing to a different σ m (not shown).The main problem is still that the value of h depends on the time since sublimation started.As we do not have this information in cloud models we cannot know easily whether the process is still in the phase where mass loss is much larger than the number loss (τ < 0.1) or already in the phase where number loss catches up (τ > 1) or somewhere in between (τ close to the inflexion point).It is clear that the time corresponding to τ 0 or τ = 1 is a characteristic time of the sublimation process.It might be useful to know it for practical applications that we consider next.For the parameters used in Fig. 5 the time scales range from 40 to 250 s, that is, time steps of cloud resolving models are mostly smaller than the characteristic times, while those of climate models are often larger.

Tests of parameterisations in two-moment models
In the foregoing sections we have seen that it is necessary to know the time since sublimation started when its effect on the number and mass concentrations is to be treated correctly, and the the problem is just that this time cannot be defined generally once other processes act in a cloud as well.
In a two-moment model one needs both quantities, f n and f m .It is easier to formulate an equation for f m than for f n because there is an analytic equation for the mass change of single crystals, while there is none such for the number (it is either zero or one, which is non-analytic in the mathematical sense).A formulation for f m can be made in a relatively straightforward way once the current mass spectrum is known or assumed.We have seen in Sect. 2 how sublimation Fig. 5. h(τ ) vs. dimensionless time τ for various initial geometric mean masses and sublimation rates: m 0 = 1 ng, a = −0.04(red); m 0 = 10 ng, a = −0.09(green); m 0 = 100 ng, a = −0.2(blue); m 0 = 1000 ng, a = −0.5 (violet).These curves are so similar that only one can be seen.The dashed black line is a generalized loglogistic function that fits h(τ ).
how fast sublimation changes from the mass loss to the number loss regime.Finally, the inflexion point is found here at 0.3 which means that the transition into the number loss regime occurs at a considerably shorter time than is needed to sublimate a crystal with the chosen characteristic mass.
We see that it is possible to find a certain universal behaviour of the sublimation curves that do not depend on initial mean mass or sublimation rate.The function does however change when we use different initial mass distributions, e.g. by changing to a different σ m (not shown).The main problem is still that the value of h depends on the time since sublimation started.As we do not have this information in cloud models we cannot know easily whether the process is still in the phase where mass loss is much larger than the number loss (τ <0.1) or already in the phase where number loss catches up (τ >1) or somewhere in between (τ close to the inflexion point).
It is clear that the time corresponding to τ 0 or τ = 1 is a characteristic time of the sublimation process.It might be useful to know it for practical applications that we consider next.For the parameters used in Fig. 5 the time scales range from 40 to 250 s, that is, time steps of cloud resolving models are mostly smaller than the characteristic times, while those of climate models are often larger.

Tests of parameterisations in two-moment models
In the foregoing sections we have seen that it is necessary to know the time since sublimation started when its effect on the number and mass concentrations is to be treated correctly, and the the problem is just that this time cannot be defined generally once other processes act in a cloud as well.
In a two-moment model one needs both quantities, f n and f m .It is easier to formulate an equation for f m than for f n because there is an analytic equation for the mass change of single crystals, while there is none such for the number (it is either zero or one, which is non-analytic in the mathematical sense).A formulation for f m can be made in a relatively straightforward way once the current mass spectrum is known or assumed.We have seen in Sect. 2 how sublimation changes the shape of the mass spectrum.This is ignored in two-moment models, such that errors arise in the mass loss computation as well as in the number loss fractions (as we will see below).
There are two principle possibilities to parameterise sublimation in two-moment models, either one defines f n as a direct function of f m or one formulates f n independently of f m .The second approach is a bit risky, since it does not guarantee that f n = (0, 1) for f m = (0, 1), so that a mixed approach that gives this guarantee would perhaps be preferable.When f n is a direct function of f m , the power law f n = f α m is the simplest formulation that fulfills f n = (0, 1) for f m = (0, 1).
In the following we test a number of potential formulations.The tests are all based on the corresponding "benchmark" functions φ n (φ m ) obtained from the analytical solution derived in Sect.2, see Fig. 2. That is, we take a formulation of f n , and compute the corresponding function φ n (φ m ) up to complete sublimation and compare the result to the benchmark φ n (φ m ).The tests use the following algorithm: 5. update mass and number, update mode mass for f (m), output of various quantities; 6. end time step loop.

Power laws
A first example showing a test of f n = f α m with α = 1.1 (Spichtinger and Gierens, 2009) is given in Fig. 6.We see that the parameterisation produces too high number loss fractions in the early phases of sublimation while in the final phase the mass loss fractions are overestimated relative to the number loss.Obviously, a choice of α<1 would deteriorate the situation, and a much larger value of α would only improve the agreement for the initial phase of sublimation at the price of much worse results in the later phases.
One could try to use a large α in the initial phase and α<1 later.The problem is, however, that one cannot decide in the bulk cloud model in which phase the sublimation process 6 K. Gierens and S. Bretl: Analytical treatment of sublimation changes the shape of the mass spectrum.This is ignored in two-moment models, such that errors arise in the mass loss computation as well as in the number loss fractions (as we will see below).
There are two principle possibilities to parameterise sublimation in two-moment models, either one defines f n as a direct function of f m or one formulates f n independently of f m .The second approach is a bit risky, since it does not guarantee that f n = (0, 1) for f m = (0, 1), so that a mixed approach that gives this guarantee would perhaps be preferable.When f n is a direct function of f m , the power law f n = f α m is the simplest formulation that fulfills f n = (0, 1) for f m = (0, 1).
In the following we test a number of potential formulations.The tests are all based on the corresponding "benchmark" functions φ n (φ m ) obtained from the analytical solution derived in Sect.2, see Fig. 2. That is, we take a formulation of f n , and compute the corresponding function φ n (φ m ) up to complete sublimation and compare the result to the benchmark φ n (φ m ).The tests use the following algorithm: e.g.
5. update mass and number, update mode mass for f (m), output of various quantities; 6. end time step loop.

Power laws
A first example showing a test of f n = f α m with α = 1.1 (Spichtinger and Gierens, 2009) is given in Fig. 6.We see that the parameterisation produces too high number loss fractions in the early phases of sublimation while in the final phase the mass loss fractions are overestimated relative to the number loss.Obviously, a choice of α < 1 would deteriorate the situation, and a much larger value of α would only improve the agreement for the initial phase of sublimation at the price of much worse results in the later phases.
One could try to use a large α in the initial phase and α < 1 later.The problem is, however, that one cannot decide in the bulk cloud model in which phase the sublimation process is, as we have seen in the foregoing sections.Thus, such a possibility is not given.

Other functions
We tested other functional relations between f n and f m as well, in particular such that give a zero derivative at f m = 0 (not shown).Examples are 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fig. 6. φn vs. φm, as parameterized in Spichtinger and Gierens (2009) (red lines, for various choices of sublimation rates and initial geometric mean masses), compared with the analytical solution for various initial geometric mean masses and sublimation rates (black line).Only one black line can be seen because the differences between the analytical solutions are smaller than the thickness of the curve.
and f n = [cos((f m − 1)π) + 1]/2.This did not yield real improvements over the simple power law.In these cases, φ n ≈ 0 until most of the ice mass is sublimated (unless α is very close to one).

Using maximum sublimating crystal mass
Alternatively one can compute in the cloud model the maximum crystal mass m max that sublimates within a timestep.Assuming the Koenig approximation this is given as With this quantity at hand we can compute while f m is computed by the model in the usual way.For the calculation of the numerator integral we may use expressions for truncated moments when analytical expressions for them are available.For our special choice of a log-normal distribution we have (cf.Jawitz, 2004): It turns out that this method produces bad results for too large and too small time steps.For too large time steps (even 1 s Fig. 6. φ n vs. φ m , as parameterized in Spichtinger and Gierens (2009) (red lines, for various choices of sublimation rates and initial geometric mean masses), compared with the analytical solution for various initial geometric mean masses and sublimation rates (black line).Only one black line can be seen because the differences between the analytical solutions are smaller than the thickness of the curve.
is, as we have seen in the foregoing sections.Thus, such a possibility is not given.

Other functions
We tested other functional relations between f n and f m as well, in particular such that give a zero derivative at f m = 0 (not shown).Examples are f n = 1− 1−f α m with α>1 and f n = [cos((f m −1)π )+1]/2.This did not yield real improvements over the simple power law.In these cases, φ n ≈0 until most of the ice mass is sublimated (unless α is very close to one).

Using maximum sublimating crystal mass
Alternatively one can compute in the cloud model the maximum crystal mass m max that sublimates within a timestep.Assuming the Koenig approximation this is given as With this quantity at hand we can compute while f m is computed by the model in the usual way.For the calculation of the numerator integral we may use expressions for truncated moments when analytical expressions for them are available.For our special choice of a log-normal distribution we have (cf.Jawitz, 2004): It turns out that this method produces bad results for too large and too small time steps.For too large time steps (even 1 s turned out too large for m 0 = 1 ng and a = − 0.04) f n is non-zero at the beginning, and for too short time steps (0.1 s) we have again the problem that φ n ≈ 0 until most of the ice mass is sublimated.

Using moments
Another possibility one could think about is to formulate f n in a cloud model via the moments of the mass distribution.These are defined as (25) Spichtinger and Gierens (2009) use dq c /dt = aµ b (corrections neglected) to calculate the tendency of the ice mixing ratio from which f m is computed.In analogy, we may set dN/dt = aµ b−1 to calculate the tendency of the number mixing ratio and compute f n from it.(Although µ b−1 is a moment of negative order, it works since the corresponding integral is finite).Although this sounds a logical approach, it is not.In this case we have f m = − aµ b /µ 1 and f n = − aµ b−1 /µ 0 , which shows that h(τ ) turns out as a constant (since the time dependent mode mass cancels out in the division f n /f m and σ m is constant).This is not only in contradiction to the analytical behaviour of h(τ ) (Fig. 5), buteven worse -it also implies f n = cf m (c a constant) which violates the boundary conditions at (1, 1), unless c = 1, which it is generally not since it depends on b (pressure and temperature dependent) and σ m .

Effect of humidity fluctuations around ice saturation
Finally we test the effect of small-scale humidity fluctuations around ice saturation (caused by random fluctuations or atmospheric waves) on sublimation.Such a situation can arise in old contrails or cirrus clouds (where relative humidity had enough time to equilibrate with the ice crystal ensemble) under certain synoptic conditions, namely when the large-scale vertical wind is very weak (less than a few cm s −1 ).Since a is a non-explicit function of time in such a case, the analytical solution is no longer valid and also the more general solution of Appendix A is not applicable.Hence we choose a simple numerical procedure.We divide the lognormal mass distribution in 1000 discrete initial masses, compute the mass growth/sublimation rates for each mass with 1 s timesteps.Once a mass gets smaller than a threshold of 0.001 ng, it is assumed to be lost.The fluctuations are represented with a sinusoidal oscillation.The sublimated ice is assumed to add to the water vapour, i.e. to increase the relative humidity.
The initial conditions are chosen such that we start with ice saturation and if we would completely sublimate the ice we would reach RH i = 105%.The oscillation amplitude is assumed to be ±5%, as well.Hence we have: Note that φ m can become negative in the growth phases.We use a geometric mean mass of 100 ng.At RH i = 95% (220 K, 250 hPa) we get a = − 9.1×10 −4 , the characteristic time (i.e.τ = 1) is 5800 s.The oscillation frequency chosen was ω = 1/250 s −1 , which is a typical value for a Brunt-Väisälä oscillation in stably stratified air.A first test of the numerical method with constant subsaturation gave φ n (φ m ) as in Fig. 2 which shows that the numerical procedure works well.When we switch on the oscillation and add the sublimated ice to the relative humidity both mass and number loss are reduced strongly; after 30000 s simulation time only a few percent of the mass and less than one per mille of the number are lost.With a lower frequency of ω = 1/2500 s −1 (gravity wave) these numbers are larger, because the initial sublimation period lasts much longer, but still less than 10 percent of the number is lost after 30 000 s.The results can be explained as follows: Initially the mean relative humidity is ice saturation, but the sublimating ice adds to this mean, such that quickly a new mean value is achieved that exceeds saturation.Hence, in spite of the oscillations that always lead transiently through subsaturated states, on the average the sublimation process is halted.After the first sublimation period all crystals that were small enough for sublimation are lost, but in the following sublimation phases the remaining crystals are large enough to survive.Unfortunately, a different behaviour is displayed by a model employing the simple parameterisation f n = f α m .Here, more than 20% of the ice is lost within 30 000 s and with ω = 1/250 s −1 , while the ice mass stays approximatly constant.This is obvioulsy the effect of adjusting the mass distribution after each time step to a new log-normal distribution (with changed geometric mean mass).Hence although the small crystals are gone within the first sublimation phase, there are new small crystals in the next, merely because the fixed distribution type enforces this.This leads to crystal loss in every new sublimation phase.This means that under conditions of a cloud in equilibrium (i.e.RH i ≈100%) small scale fluctuations in the cloud model can lead to an artificial crystal loss with corresponding consequences on further microphysical evolution and optical properties.

Discussion
Having performed the simple tests of the various possible parameterisations we see now clearer the difficulties inherent in modelling ice sublimation in the framework of two-moment www.atmos-chem-phys.net/9/7481/2009/Atmos.Chem.Phys., 9, 7481-7490, 2009 models.These difficulties arise in the same way when another mass distribution is chosen (we have also tested a generalised gamma distribution).An exponential distribution which is typical for cloud droplets could give different test results; this has not been tested in this paper because we are mainly interested in improvement of cirrus and contrail models.But such tests could be done in an analogous manner.
As a general rule of thumb the analytical solutions showed that over a large fraction of the sublimation phase the number loss fraction equals roughly the mass loss fraction, as seen for instance in Figs. 2 and 3.In particular for large-scale models with long timesteps of the order of the sublimation timescale it seems therefore reasonable to set f n = f m .For cloud resolving models which need timesteps much shorter than the sublimation timescale problems can arise, in particular in the beginning of the process, when the number concentration is hardly affected although the mass concentration is already diminishing.Hence it is not easy to find a better parameterisation for a small-scale model.The simple parameterisation f n = f α m with α 1 seems to obey the simple rule, but in the oscillation case it went wrong although the more exact simulation showed that the rule is valid in this case as well.
Unfortunately, the oscillation case is not so academic like the other tests we have performed; we even took into account that the sublimating ice mass contributes to the vapour phase and thus to the relative humidity in the cloud.The oscillation case is typical for mature clouds where crystals have consumed the excess vapour completely and where random fluctuations (turbulence) and atmospheric waves of various kind affect the cloud evolution.In order to avoid artificial number concentration loss in such quasi-equilibrium cases one could either employ a larger value of α, such that larger mass loss is needed before crystal loss commences.Or one could introduce a sublimation humidity threshold of several percent below saturation.Both strategies have their disadvantages.A sublimation threshold several percent below 100% is unphysical and may have unforeseeable effects.Larger α underestimates crystal loss in situations of steady substantial sublimation (as do other functional relations, as we have seen); however, one can let α depend on the degree of subsaturation, with larger values at small and smaller values at larger saturation deficits.This has been tested as well and it produces better results than the parameterisation with constant α.
The question arises whether there are other pieces of information available at timestep level in the two-moment model.We have: two independent moments, namely number and mass concentrations, and direct functions of them (all other moments, effective sizes, and so on).More information can only be provided by the thermodynamic state at the timestep, i.e. the relative humidity and the temperature.These quantities allow to additionally compute the maximum mass that can evaporate within one timestep.However, we have seen that even with this additional information it is difficult to construct a better parameterisation of sublimation; at least our approach in the previous section was not successful although it sounded plausible.
Our analytical tests were admittedly a bit academic.The analytical solution is only valid for constant a, which would mean that sublimated ice would disappear from the system.In reality, the sublimating ice increases the vapour concentration which generally will lead to increasing a.Also the sublimation process was considered in isolation, which is reasonable for testing, but in reality other processes act simultaneously.The sublimation timescales can be very long especially when the subsaturation is low.In such a case, however, there may well be other processes with shorter timescales that then dominate the cloud evolution rendering the problems with sublimation less important.When subsaturation is larger or the mean mass smaller, sublimation time scales get shorter; then the sublimation process will quickly evolve into that regime where f n ≈f m , such that the problems with the initial phases lasts for a shorter period of time.Here it would turn out disadvantagous to have a larger α or one of the other approaches discussed in the previous section, because they lead to underestimation of the number loss over a large fraction of the whole sublimation process.
Considering these arguments we think that, in spite of its weaknesses, the simple approach f n = f α m with α slightly exceeding unity or dependent on the degree of subsaturation (or unity for large-scale models) is still one of the best choices one can make.

Conclusions
The analytic solution of the spectral form of the equation for deposition/sublimation and its application to ice sublimation gave the following insights: -The relation between fractional number loss and fractional mass loss is non-unique, that is, it is no function.However it seems advantageous to formulate it as a function in parameterisations of ice microphysics for both large-scale and small-scale models.
-Sublimation timescales are usually longer than timesteps of cloud resolving models, and often but not generally shorter than timesteps of large-scale models.The latter case is easier to parameterise than the first one.
-As a rule of thumb the analytical solutions showed that over most of the sublimation process the number loss rate equals approximatly the mass loss rate.It is therefore a good idea to set them equal in large-scale models with long timesteps (say, 15 min.and more).
-For small-scale models (e.g.cloud-resolving and largeeddy models) we have tested a number of alternatives to the standard power law, but the power law turned out to produce the best results compared to the analytical solution.
-Care has to be taken with the power law formulation when a cirrus cloud or condensation trail at ice saturation undergoes small-scale humidity fluctuations (probably irrelevant for large-scale models).Without any counter-measure the power law with α close to one leads to severe overestimation of crystal loss.It might be better to let α be a function of saturation deficit with larger values at low deficit and values closer to one at larger subsaturations.
The academic test cases of the present paper can only give insight about the problems inherent in parameterisations of ice sublimation.The success of an alternative parameterisation in "every day work" should, however, be checked against spectrally resolving microphysics models where sublimation acts in combination and in competition with all other processes.
The analytical solution points to the very core of the problem: the sublimation process requires an additional parameter for the complete description of the mass spectrum.The additional parameter must be obtained from an additional equation which is difficult to formulate in general once a variety of other processes affect the shape of the spectrum as well.To introduce such new equations is probably not the best way to overcome the problems since this would add a significant amount of complexity to the bulk models.If exact solutions are required one should resort to bin models which do not suffer from the mentioned problems.Nevertheless, bulk models are useful and needed whenever there are computer time and memory constraints, which is a common situation.Therefore, in spite of the problems, we encourage their use, provided that their limits are thoroughly checked.ẋ transformed mass growth rate x 0 transformed geometric mean mass α power law exponent for the sublimation parameterisation α fit parameter t model time step µ k moment of order k φ 0 ≡ φ n total fractional number loss φ 1 ≡ φ m total fractional mass loss σ m geometric standard deviation of mass spectrum σ x geometric standard deviation of transformed mass spectrum ω oscillation frequency

Fig. 2 .
Fig.2.Total fraction of number loss, φn, vs. total fraction of mass loss, φm, for various initial geometric mean masses from 1 to 1000 ng.Note that the curves are almost congruent.Time runs along the curves from (0, 0) to (1, 1).

Fig. 2 .
Fig.2.Total fraction of number loss, φn, vs. total fraction of mass loss, φm, for various initial geometric mean masses from 1 to 1000 ng.Note that the curves are almost congruent.Time runs along the curves from (0, 0) to (1, 1).

Fig. 2 .
Fig. 2.Total fraction of number loss, φ n , vs. total fraction of mass loss, φ m , for various initial geometric mean masses from 1 to 1000 ng.Note that the curves are almost congruent.Time runs along the curves from (0, 0) to (1, 1).

Fig. 3 .
Fig.3.Mass loss (f m ) and number loss (f n ) fractions for series of subsequent sublimation time steps.Time steps are 100, 200, 300, 400 s (solid red, green, blue, dashed red), and time runs along each curve from the point neares to the fm -axis to the point nearest to the f n -axis.The individual time steps are marked on each curve.Arrows indicate the direction of time.The yellow point at (1, 1) was the result of a test with a very long time step that should guarantee that the ice sublimates completely within that step.The black line is the attempt of a fit of the later parts of the curves, fn=f 0.89 m .The initial geometric mean mass for the calculations is m 0 =1 ng.

Fig. 4 .
Fig.4.f n vs. f m for various initial geometric mean masses from 1 to 1000 ng (dashed red, solid green, blue, red).Note that the inverse values of the axis correspond to the timescales for the change of the corresponding integrals I k (t) in seconds.Arrows indicate the direction of time.

Fig. 3 .
Fig.3.Mass loss (fm) and number loss (fn) fractions for series of subsequent sublimation time steps.Time steps are 100, 200, 300, 400 s (solid red, green, blue, dashed red), and time runs along each curve from the point neares to the fm -axis to the point nearest to the fn-axis.The individual time steps are marked on each curve.Arrows indicate the direction of time.The yellow point at (1, 1) was the result of a test with a very long time step that should guarantee that the ice sublimates completely within that step.The black line is the attempt of a fit of the later parts of the curves, fn=f 0.89 m .The initial geometric mean mass for the calculations is m0=1 ng.

Fig. 4 .
Fig.4.f n vs. f m for various initial geometric mean masses from 1 to 1000 ng (dashed red, solid green, blue, red).Note that the inverse values of the axis correspond to the timescales for the change of the corresponding integrals Ik(t) in seconds.Arrows indicate the direction of time.

Fig. 4 .
Fig.4.f n vs. f m for various initial geometric mean masses from 1 to 1000 ng (dashed red, solid green, blue, red).Note that the inverse values of the axis correspond to the timescales for the change of the corresponding integrals I k (t) in seconds.Arrows indicate the direction of time.