Clouds – Part 1

Abstract. A double-moment bulk microphysics scheme for modelling cirrus clouds including explicit impact of aerosols on different types of nucleation mechanism is described. Process rates are formulated in terms of generalised moments of the underlying a priori size distributions in order to allow simple switching between various distribution types. The scheme has been implemented into a simple box model and into the anelastic non-hydrostatic model EULAG. The new microphysics is validated against simulations with detailed microphysics for idealised process studies and for a well documented case of arctic cirrostratus. Additionally, the formation of ice crystals with realistic background aerosol concentration is modelled and the effect of ambient pressure on homogeneous nucleation is investigated in the box model. The model stands all tests and is thus suitable for cloud-resolving simulations of cirrus clouds.


Introduction
The role of clouds is crucial for our understanding of the current and the changing climate (IPCC, 2007).Cirrus clouds modulate the Earth's radiation budget significantly.It is assumed that (thin) cirrus clouds contribute to a net warming of the Earth-Atmosphere system (e.g.Chen et al., 2000), but the magnitude of this warming has not been quantified yet.Recently, the impact of thin cirrus clouds in the mid latitudes was estimated in idealised framework using vertical profiles from radiosondes (Fusina et al., 2007), but the global effect is still uncertain.The formation and evolution of cirrus clouds depends in a complex way on a variety of environmental conditions (temperature, relative humidity, wind fields) as well as on the impact of background aerosol acting as ice nuclei.
Correspondence to: P. Spichtinger (peter.spichtinger@env.ethz.ch)The interaction of various processes and their non-linear dependence on ambient conditions renders the understanding of cirrus clouds in general a difficult task.
Cirrus clouds (except anvils) are closely related to their formation regions, so-called ice-supersaturated regions (ISSRs, see e.g.Gierens et al., 1999).These are large, initially cloud free airmasses in the upper troposphere (and sometimes lowermost stratosphere) in the status of supersaturation with respect to (wrt) ice.These regions are quite frequent in the tropopause region (see e.g.Spichtinger et al., 2003a,b;Gettelman et al., 2006).From former investigations (Spichtinger et al., 2005a) it turns out that cirrus clouds often are embedded in horizontally extended ISSRs; ISSRs and their embedded clouds form a system.Largescale dynamical processes like synoptic upward motions, but also mesoscale waves and small scale turbulence play crucial roles for the formation and evolution of the system ISSR/cirrus (Spichtinger et al., 2005a,b).Local dynamics and microphysics are acting on the cloud and sub cloudscale.From this point of view there is need of a cloud resolving model which can be used for idealised studies of cirrus clouds interacting with various scales of dynamics.
While from theory and measurements it is quite understood that in cloud free air masses the relative humidity wrt ice can reach very high values up to the freezing thresholds for homogeneous freezing (i.e.140-170% RHi, depending on temperature, see Koop et al., 2000), substantial and persistent supersaturation inside cirrus clouds is more difficult to understand.Ice crystals act as a strong sink for water vapour.Thus, one expects that RHi-distributions inside cirrus are centred around saturation, as some measurements indicate (Ovarlez et al., 2002;Spichtinger et al., 2004).However, there are also many measurements from inside cirrus clouds that indicate considerable degrees of supersaturation (see Comstock et al., 2004;Lee et al., 2004;Ovarlez et al., 2002;Krämer et al., 2008;Peter et al., 2008, and M. Krämer, personal communication).These findings seem to Published by Copernicus Publications on behalf of the European Geosciences Union.P. Spichtinger and K. Gierens: Modelling cirrus clouds -Part 1a be contrary to our current understanding of microphysics inside cirrus clouds (Peter et al., 2006) and call for an explanation.One possible explanation will be given in the present paper.Another pending question is the impact of different nucleation mechanisms on the formation and evolution of cirrus clouds.While it is generally assumed that homogeneous nucleation is the dominant formation process for cold (T <−38 • C) cirrus clouds (e.g.Sassen and Dodd, 1988;Heymsfield and Sabin, 1989;Haag et al., 2003b), there are indications that heterogeneous nucleation can substantially modify the conditions for homogeneous nucleation bringing forth large change in resulting cloud properties.Therefore a model for studying the competition of different nucleation processes would be useful.
We have developed a new ice microphysics scheme for the use in box models and cloud resolving models, based on earlier work (Gierens, 2003).A novel feature of the model is the use of arbitrary many classes of ice, discriminated by their formation mechanism.The corresponding aerosol types that are involved in the formation of the various ice classes are treated as well.This new concept allows us to investigate the impact of different nucleation processes in the same airmass, in particular how air pollution (heterogeneous nucleation) affects the cloud evolution.Another novel feature is the formulation of the various process rates in terms of moments of the underlying crystal size distribution.This makes it possible to the choose between various distribution types.
The structure of this article is as follows: in the next section we will describe the basic dynamical models, i.e. the box model and the anelastic, non-hydrostatic model EULAG very briefly.In Sect. 3 the new microphysics scheme is described in detail.In Sect. 4 the validation of the model using box model simulations and 1-D simulations is shown, and the model's performance is discussed.We end with a summary and draw conclusions in Sect. 5.

Model description -dynamics
The new microphysics scheme is implemented into two different types of models: First, we implemented the ice micro-physics into a simple box model for validating the nucleation parameterisation and for fast calculations serving a principal understanding of the interaction of different processes.The box model can also be coupled to trajectories, e.g. to the output of a trajectory model (in our case LAGRANTO, Wernli and Davies, 1997).In a second step we implemented the tested microphysics into the anelastic, non-hydrostatic model EULAG (Smolarkiewicz and Margolin, 1997).The twofold approach was not only for testing the model but also to have two different tools which can be used for different applications, which have exactly the same ice microphysics parameterisations.
In the following we describe first the more complex dynamics of the EULAG model and the coupling of the dynamics to the microphysics.
Then, we describe the box model which was developed together with the microphysics in the spirit of the EULAG model, i.e. using background states (e.g. for potential temperature) as well.This choice was made to make the "transition" between the two models as smooth as possible.

EULAG model -dynamics
As a basic dynamical model we use the anelastic nonhydrostatic model EULAG (see e.g.Smolarkiewicz and Margolin, 1997).The anelastic equations for the dry dynamics can be written in perturbation form as follows (cf.Smolarkiewicz et al., 2001;Grabowski and Smolarkiewicz, 2002) Here, u is the velocity vector; p, ρ and θ denote pressure, density and potential temperature, respectively; g and f denote gravity and "Coriolis" vectors, respectively; θ and ρ are the anelastic reference state profiles for potential temperature and density; M denotes additional appropriate metric terms, depending on the coordinate system chosen.The subscript e refers to the environmental profiles, which must not necessarily be equal to the reference states.Primes denote deviations from the environmental state (e.g.θ =θ −θ e ).D Dt : =∂/∂t + u • ∇ denotes the total derivative.The perturbation pressure p is calculated using the mass continuity constraint ∇ • ρu=0.
For solving the governing equations we use the unified semi-Lagrangian-Eulerian approach described in Smolarkiewicz and Margolin (1997); Smolarkiewicz et al. (2001): Let and F denote the vectors of variables (θ , u, v, w) and their forcings, respectively.With ˜ : = n + 0.5 tF n and the generalised transport operator LE the approximation can be described as whereby i and n denote the spatial and temporal location, respectively.This results into a trapezoidal rule for the approximations.This treatment and the non-oscillatory forwardin-time (NFT) semi-Lagrangian/Eulerian approximations of the integrals were carried out as described in detail in Smolarkiewicz and Margolin (1997) and Smolarkiewicz et al. (2001).The model was used for many applications on different scales and several problems in atmospheric dynamics (e.g.stratified flow over mountains, convectively induced gravity waves etc.).One main advantage of this model is the less diffusive advection scheme MPDATA (Multidimensional Positive Definite Advection Transport Algorithm, see e.g.Smolarkiewicz and Margolin, 1998).
For including cloud physics into the model, we have to perform "moist" dynamics and a coupling of dynamics and thermodynamics.This can be done as follows (see also Grabowski and Smolarkiewicz, 2002): We define the density potential temperature (see Emanuel, 1994) including specific humidity q v and the mixing ratio of cloud ice q c as follows: with p =(1/ )−1, where =R g /R v denotes the ratio of the ideal gas constants for dry air and water vapour, respectively.Using the definition of an environmental density potential temperature θ de : =θ e + θ p q ve for the representation of the perturbation θ d =θ d −θ de the governing equations for the moist dynamics read as follows: The coupling of dynamics and thermodynamics manifests itself in two parts: 1.An additional buoyancy source from the deviations in water vapour and the cloud ice in Eq. ( 5) in the density potential temperature 2. An additional source term F θ on the right hand side of Eq. ( 6) due to diabatic processes (phase changes etc.).

Box model
The (zero-dimensional) box model represents an air parcel which is moved in the vertical direction with a velocity w=w(t), prescribed for the whole simulation time t s .Here, we assume only adiabatic processes, i.e. the vertical velocity produces an adiabatic cooling (expansion) or warming (compression) for the background temperature T e due to The governing equations for the dynamics reduce to i.e. the coupling of dynamics and thermodynamics is reduced to the additional diabatic forcing term.As indicated above, for consistency with the formulation of the dynamics of EULAG we have formulated the box model using environmental states T e , p e , θ e , i.e. the model is formulated in perturbation form: All diabatic processes will only change θ =θ−θ e , the adiabatic processes change the environment, i.e.T e , p e while θ e remains constant.This concept can also be used for the dynamical model EULAG, where we can simulate adiabatic cooling due to upward motion by changing the background physical temperature T e .

Model description -ice microphysics
In this section the newly developed ice microphysics scheme is described.First, we define the set of variables and derive the governing equations.The classes of aerosol and ice crystals will be defined, as well as their properties and mass (or size) distribution types.The moments of these distributions will be used in the microphysical equations.Finally, the currently implemented microphysical processes (nucleation, deposition, sedimentation) are described in detail.For all quantities we use SI units.Aggregation is not yet implemented in the microphysics scheme.However, aggregation is of less importance for the cold temperature regime (T <−40 • C) and/or for moderate vertical velocities (Kajikawa and Heymsfield, 1989).The development of an appropriate parameterisation consistent with the two-moment scheme is subject to ongoing research.

Prognostic equations
As we are interested in particular in the interplay of various competing nucleation modes acting in a cirrus cloud (in particular homogeneous vs. heterogeneous nucleation), we allow a non-specified number of aerosol and ice classes.Each aerosol class corresponds to an ice class that it nucleates and vice versa.Each class has a number and a mass concentration, N x and q x , respectively, which are the zeroth and first moments of a mass distribution f (m).Note that N x is meant in a mass specific sense, that is, the unit of the number concentration is kg −1 (per kg of dry air).We use this convention to be consistent with the formulation of the advection equations in flux form.This formulation will be used for the sedimentation of ice crystals, see below.However, ice crystal number densities and ice water content can be derived using the ambient density ρ: n c =N c • ρ, IWC =q c • ρ.The kind of the distribution is preselected (e.g.log-normal or gamma etc.).The moments of f (m) are with the normalisation µ 0 =N x .Prognostic equations for the number and mass concentrations of each aerosol and ice class form the basis of our bulk microphysical two-moment scheme.
Note that we prefer mass distributions whereas observers present their measurements usually as size distributions.Ice crystals in nature appear in a myriad of shapes.Crystal mass is therefore a more convenient choice for model formulation than crystal size.Therefore we formulate the prognostic equations for masses; sizes and shapes are diagnosed from them.
Ice crystal shape is important for depositional growth (Stephens, 1983), sedimentation (Heymsfield and Iaquinta, 2000) and radiative properties (Wendisch et al., 2005(Wendisch et al., , 2007) ) of ice crystals.The ice crystal shape depends clearly on temperature and ice supersaturation (e.g.Bailey and Hallet, 2004;Libbrecht, 2005) but columns seem to be a frequent habit below −40 • C in a variety of field measurements (Heymsfield and McFarquhar, 2002, their Table 4.1).Therefore we assume columnar ice crystals in our model.
One essential difference of our scheme to other schemes (e.g.Kessler, 1969;Seifert and Beheng, 2005) is that we do not differentiate ice classes according to their size.Traditionally, cirrus ice was classified as cloud ice and snow, where cloud ice consisted of ice crystals that are so small that their terminal fall speed could be set to zero, while snow was the ice fraction that had non-negligible fall speeds.This classification allowed a better treatment of sedimentation in single-moment models.As a two-moment model partly overcomes the problems with sedimentation by introducing two sedimentation fluxes, for number and mass concentration, respectively, we no longer differentiate between cloud ice and snow.Instead our ice classes correspond to various nucleation processes (or to the respective aerosol class).
The processes ice crystal nucleation, depositional growth/sublimation and sedimentation are currently implemented in the scheme.The prognostic equations for potential temperature and the microphysical variables for n classes of ice (index c) and aerosol (index a) are thus: L θ e c p T e (NUC + DEP) + D θ (10) DN c,j Dt = Dq a,j Dt = NUCA j + DEPA j + D q a,j where j is the respective class index and Here, θ denotes potential temperature, T e and θ e denote the physical and potential temperature of the environmental state, respectively.L s is the latent heat of sublimation, c p denotes the specific heat at constant pressure and ρ denotes the density of the reference state.q v is the specific humidity.D ψ terms represents sources and sinks of the model variable ψ due to processes not explicitely represented in the equations (e.g.turbulent transport, diffusion etc.).v m,j and v n,j are the mass and number weighted terminal velocities (see Sect. 3.4), respectively.The terms NUC j and NNUC j represent the sources and sinks for the ice crystal mass and number concentration, respectively, due to nucleation.Sources and sinks related to diffusional growth or sublimation are represented by the terms DEP and NDEP.We assume that every ice crystal nucleates from an aerosol particle (for both nucleation types, heterogeneous and homogeneous), hence the aerosol particle will be removed (i.e. it is now inside the crystal) from the aerosol number concentration and it will be released if the ice crystal sublimates completely.Therefore we can treat the sources and sinks for the aerosol number concentration using the same terms NNUC j and NDEP j .For treating the sources and sinks of the aerosol mass concentration we use the terms NUCA j and DEPA j .Note, that there is no sedimentation of the aerosol particles (due to their small masses).

Crystal shape and mass-size relation
For parameterisation of the various processes we need assumptions on the properties of single ice crystals, for instance their shapes.Generally, we assume that ice crystals are hexagonal columns (Bailey and Hallet, 2004), with height L and diameter D (twice the side length of the hexagon).
The aspect ratio r a : =L/D depends on crystal size, such that small crystals have r a =1 and larger crystals have r a >1.Mass-length relations from Heymsfield and Iaquinta (2000) of the form L(m)=α(m) m β(m) with piecewise constant parameter functions α(m), β(m) (Table 1) are used to derive crystal size from the prognostic quantity crystal mass.The bulk ice density is assumed as ρ b =0.81 • 10 3 kg m −3 (Heymsfield and Iaquinta, 2000).The parameter values are such that the boundary between small and large crystals is m t =2.146 • 10 −13 kg, or equivalently L t =7.416 µm.The aspect ratio can be formulated as a function of m; using the formula for the volume of a hexagonal column (V =( √ 27/8)•D 2 •L) and the mass-length relations we arrive at:  In Fig. 1 length and diameter of the ice crystals are shown as functions of crystal mass, and in Fig. 2 the aspect ratio vs. the ice crystal mass and length is shown.

Sedimentation and mass-fall speed relation
For formulating the sedimentation fluxes of the ice mass and number concentrations, we need assumptions on the terminal velocity of a single ice crystal.Here again we use an ansatz by Heymsfield and Iaquinta (2000): Using the mass-length relations derived above we formulate the terminal velocity as a function of crystal mass: with piecewise constant parameter functions γ (m), δ(m), as given in Table 2.The parameters have been derived using the coefficients for small columns by Heymsfield and Iaquinta (2000).The terminal velocities v 0 are valid for reference values of T 0,1 =233 K, p 0,1 =300 hPa, for other temperatures and pressures we apply the correction factor such that v(m, T , p)=γ (m) m δ(m) c(T , p).For very large ice crystals (L>1899 µm) we use coefficients adapted from   Barthazy and Schefold (2006), who used a reference state of T 0,2 =270 K, p 0,2 =815 hPa.This second reference state was taken into account in the derivation of our coefficients, which are then valid for the whole tropospheric temperature and pressure range.The terminal velocity of a single ice crystal vs. its mass is shown in Fig. 3.

Choice of distribution type
For the formulation of the process rates in Eqs.(10-15) it is required to choose an a priori distribution type for the crystal masses and the masses of single aerosol particles.This choice can be handled very flexible (for instance as an entry in a FORTRAN namelist) when we are able to formulate the process rates in terms of general moments of the distribution.Then specification of a distribution in the model simply means branching into the corresponding function for the general moments.Because of this flexibility it does not matter much which kind of distribution we use here to derive the process rates.However, the model formulation is simpler when we select a lognormal distribution than when we select something else.The reason for that is the special form of the mass-size and mass-fall speed relations that we have assumed and that are common to cloud microphysics formulations at least since the 1970s.These relations have the form of a power law.It can easily be shown that when two quantities are related by such a power law and one quantity is lognormally distributed, then the other quantity is lognormally distributed as well.There are other distributions that have that property as well, e.g.Weibull and generalised gamma distributions.The former of these is a generalisation of the exponential distribution and might sometimes apply, e.g. for young contrails (Gierens, 1996).The latter has more parameters then the lognormal, which has two.Since it is difficult to derive equations for all parameters, one usually has to fix them somehow in an ad hoc way.Hence, we prefer to have a distribution with few parameters that avoids such procedures as much as possible.When a gamma distribution is selected for crystal sizes, as is often done, then masses and terminal velocities are generalised gamma distributed.However, from observations there is no evidence for preferring (generalised) gamma distributions against lognormal distributions for fitting size distributions of ice crystals.For many observational data, a lognormal distribution can be fitted as well as a generalised gamma distribution (M.DeReus, personal communication), and indeed, lognormal fits were used in several studies (e.g.Schröder et al., 2000) Also the observation of multiple size modes does not call for multi-parameter distributions.We strongly believe that multiple modes are the result of mixing of ice crystal populations that have different origin.In our model such cases are covered by the use of multiple ice classes, where the combined size distribution of several classes will often show several modes.Hence, we prefer the lognormal distribution.
The lognormal distribution for the ice crystal mass can be written as with geometric mean mass m m and geometric mass standard deviation σ m .The lognormal distribution is completely specified once its zeroth, first, and second moment are given.The prognostic variables of the two-moment scheme let the second moment a free parameter that we either have to fix to a constant or to make a function of the mean mass.We use the latter possibility.For this purpose we follow Höller (1986) and define a "predominant mass": divide it by the mean mass m : =µ 1 [m]/µ 0 [m] and set the ratio r 0 of these masses constant, that is: Hence, the geometric standard deviation can be expressed as Using Eq. ( 22) the analytical expression for the moments of the lognormal distribution f (m) can be written as follows: Obviously the formulation using the ratio r 0 is much simpler than the formulation using σ m .Given the lognormal mass distribution, the corresponding lognormal distributions for the related quantities size and fall speed are obtained using the following transformation (e.g. for crystal length L): Here we have considered the parameters as constant for the sake of simplicity.Since the coefficients are actually piecewise constant functions of mass, the transformation formula above is not strictly correct.One possible correction is the use of truncated moments (see below).
In Table 3 the values of σ m and σ L depending on the ratio r 0 and the mass range are shown.

Nucleation
Two different nucleation processes are parameterised in our scheme: First, homogeneous freezing of supercooled aqueous solution droplets (Sect.3.2.1)and second, heterogeneous freezing on solid aerosol particles (Sect.3.2.2).Both mechanisms depend on relative humidity wrt ice; for calculating relative humidities from prognostic variables (temperature, specific humidity, pressure), we use formulae from Murphy and Koop (2005) for the saturation pressures of (supercooled) water and ice.

Homogeneous nucleation
The solute mass (mass of H 2 SO 4 ) in a solution droplet is equivalent to the radius r of a sphere of the pure solute.For this radius we prescribe a lognormal size distribution f (r) with the geometric mean radius r m and a geometric standard deviation σ r .
Given the solute mass, the radius r d of the solution droplet is obtained from the Köhler theory, depending on temperature and relative humidity.From the Köhler theory in its simplest form we know that at (water) saturation the equilibrium droplet radius r d is proportional to the square root of the solute mass, i.e. proportional to r 3/2 , a power law relationship.Hence the solution droplets are approximately lognormally distributed as well, which is in accordance to measurements of upper tropospheric aerosol (Minikin et al., 2003).
The probability for an aqueous H 2 O/H 2 SO 4 solution droplet of volume V d to freeze within a time period t is where J hom ( a w , T ) denotes the homogeneous nucleation rate which is parameterised according to Koop et al. (2000) in terms of temperature and a w : =a w −a i w , the difference of water activity in the solution and a i w , the activity of the water in the solution in equilibrium with ice at temperature T (i.e. a i w =e * i (T )/e * w (T ), the ratio of the saturation vapour pressures wrt ice and liquid water, respectively).a i w is thus a function of temperature alone.The water activity itself is the ratio of the equilibrium vapour pressure over the solution to the equilibrium vapour pressure over pure liquid water.
When dynamical processes (uplift) are slow enough that solution droplets can equilibrate their volume to changes in relative humidity, then the water activity equals the saturation ratio (wrt water).This is assumed here.Hence, the homogeneous nucleation rate J hom can be expressed as a function of relative humidity and the temperature.The number of ice crystals created within a time step t from an initial aerosol concentration N a can now be calculated as (where f a (r) is the aerosol size distribution normalised to 1) and the frozen ice water content can be calculated as where w w denotes the H 2 O weight fraction and ρ d the density of the solution, respectively.Exploiting the lognormal character of the dry aerosol size distribution, we can use a Gauss-Hermite integration for numerical calculation of the integral (Gierens and Ström, 1998).The integral in Eq. ( 28) is usually much less than one, so that homogeneous nucleation usually is not number limited.Note, that homogeneous freezing of solution droplets only occurs at T <−38 • C, i.e. below the supercooling limit of pure water.Formation of N c ice crystals simply implies loss of N a = − N c aerosol particles.However, a priori it is not clear how much aerosol mass is transferred to the ice in the nucleation process.Here we use the following procedure.While the aerosol number concentration decreases by a factor f n =| N a |/N a the dry aerosol mass concentration in the aerosol class decreases by a factor f m .Hence the mean dry aerosol mass is reduced by a factor (1−f m )/(1−f n ).We can compute and apply this factor once a relation between f m and f n is given.For this we make the ansatz The postulate α≥1 expresses that fact that large droplets (consisting of large aerosols) will freeze first and vanish from the aerosol pool (see e.g.Haag et al., 2003a, Fig. 8).It turns out from additional calculations, that a factor α=1.33 is a good approximation.We will use the same approach for the sublimation of ice crystals, see Sect.3.3.Note that we let the width parameter (σ m ) of the aerosol mass distribution unchanged during the nucleation event.
The shift of the mean mass of the background aerosol distribution is optional and can be switched on and off.For very large background aerosol concentrations, the shift of the mean mass is marginal and can be switched off.Then, only the number concentration is decreased, while the mean mass or size of the aerosol distribution remains constant; the new mass distribution is calculated using the new number concentration.For simulations with high vertical velocities at very cold temperatures (T <210 K) the shift of the mean mass is recommended.
The impact of shifting the modal radius of the background aerosol distribution leads to smaller radii of solution droplets during the ongoing nucleation event, resulting in a decrease of number density of nucleated ice crystals due to the abrupt change during a nucleation event.

Heterogeneous nucleation
Although the investigation of different nucleation mechanisms within the same environment is one of the key issues which will be studied using the new cirrus microphysics, we mention heterogeneous nucleation only briefly here, because we will not use it further in the present paper.Generally, the parameterisation for heterogeneous nucleation makes use of prescribed background aerosols, which act as ice nuclei.Thus, there is an explicit impact of the aerosol on the formation of ice crystals and also washout of the background aerosol trapped in the sedimenting ice crystals is described.After nucleation the aerosol particles are trapped inside the ice crystals and are released as soon as the ice crystals evaporate.
In our model we can use parameterisations of heterogeneous nucleation of any type, more or less sophisticated.In Part 1b (Spichtinger and Gierens, 2009) we will extensively use heterogeneous and homogeneous nucleation within the same environment; the parameterisation that is used is described there.

Deposition growth and sublimation
The growth Eq.(see e.g.Stephens, 1983) for a single ice crystal of mass m reads: where ρ v (T e ) and ρ s,i (T s ) denote the ambient water vapour density, derived from ambient humidity and temperature T e and the saturated (with respect to ice) water vapour density at the ice crystal surface, i.e. at surface temperature T s .The other factors in Eq. ( 30) are as follows: -Diffusivity of water vapour in air according to Pruppacher and Klett (1997) using reference values T 0 =273.15K and p 0 =101325 Pa, respectively.
-Capacitance factor C which accounts for the nonspherical crystal shape.We could have used the capacitance factors for hexagonal columns of Chiruta and Wang (2005) here, but we did not for the following reasons: First, the well-known uncertainty in the deposition coefficient (see below) has a larger effect on the results than the choice of C. Second, the formulation of ventilation factors (see below) has been derived for spherical water drops and only few experiments have been carried out for other shapes than spherical Hall and Pruppacher (1976).Third, an ice cloud contains anyway a mixture of crystal shapes and habits.Inclusion of the correct capacitance factor for one habit leads to inconsistencies with other habits.For these reasons we simply followed Hall and Pruppacher (1976).For spheres of radius r, C=r.Hexagonal columns with length L=2a and diameter D=2b can be approximated by prolate spheroids with semi axes a and b (a≥b), i.e. with an aspect ratio r a =L/D=a/b.The capacitance factor C can be determined using the electrostatic analogy (McDonald, 1963) a denotes the eccentricity of the spheroid.We use the aspect ratio r a introduced in Sect.3.1 for our calculations.
-Correction factor f 1 for the difference between the transfer of water molecules to the crystal by pure diffusion and that according to kinetic treatment of individual water molecules (important for very small crystals with sizes less than 1 µm): where A is the surface area of the ice crystal, M w denotes molecular weight of water and R is the universal gas constant.α d denotes the deposition coefficient.Currently there is no agreement on a generally accepted value of α d .Measured values range between 0.004≤α d ≤1 (see e.g.Pruppacher and Klett, 1997;Magee et al., 2006), however most models work well with α d ≥ 0.1 (Lin et al., 2002) for reasons discussed in Kay and Wood (2008).The deposition coefficient could even depend on crystal size (e.g.Gierens et al., 2003) or on ice supersaturation (Wood et al., 2001).
For our validation runs we have set α d =0.5 (see e.g.Kärcher and Lohmann, 2002a,b).
-Ventilation factor f 2 to correct for the enhanced growth of ice crystals due to enhanced water vapour flux arising from motion of the crystal relative to the environmental air (important for large crystals).We follow Hall and Pruppacher (1976) and set: Re ≤ 1 0.86 + 0.28(N where N Sc =Dη/ρ denotes the Schmidt number, η is the viscosity of air and N Re is the Reynolds number defined by characteristic dimensions of the ice crystal.
The latent heat released on the growing crystal must be diffused to the ambient air.This is described by an analogous equation: L s denotes the latent heat of sublimation, the coefficients f * 1 , f * 2 are the counterparts to those in Eq. ( 30): Here, K is the thermal conductivity of moist air, u r is the average thermal velocity of air molecules striking the ice surface and β d =1 is the thermal accommodation coefficient.
-The ventilation factor f * 2 for thermal diffusion is calculated as follows (Hall and Pruppacher, 1976): Re,l *

Q
, where N P r denotes the Prandtl number and N Re,l * Q denotes the Reynolds number for the characteristic length l * Q .
In order to compute the parameters required for a parameterisation, the two diffusion equations are solved iteratively using a fourth order Runge-Kutta scheme.For crystals of intermediate size, i.e. above the kinetic regime but still small enough such that ventilation is negligible, the ansatz by Koenig (1971) dm dt ≈am b , provides a good approximation to the numerical solution.The coefficients a, b depend on saturation ratio, temperature and pressure and can be derived using a least squares regression of exponential type f (x)=ax b .The coefficients a, b have been calculated for pressures in the range 150≤p≤600 hPa in 50 hPa bins and for the temperature range −80≤T ≤−20 • C in 1 K bins.For the Runge-Kutta integration we have assumed water saturation.The actual value of a during a simulation is the product of the tabulated value and (e−e * i )/(e−e * w ).This factor is negative for ice-subsaturated conditions, i.e. the equations then describe crystal sublimation.The functions a(p, T ) and b(p, T ) for selected pressures are displayed in Fig. 4. The form of the Koenig approximation makes it ideal for later use in the prognostic equation for q c .However, it overestimates strongly the crystal growth in the kinetic regime and underestimates the influence of the ventilation factor for large crystals.In order to overcome these problems, we introduce a correction of the following form: where m 0 =m 0 (p, T ) ∼ 10 −16 − 10 −14 kg, γ =γ (p, T ) ∼ 0.2−0.25,m l =2.2 • 10 −10 kg, and δ=0.12.
Using these corrections we are able to approximate the mass growth rates for single crystals within an error margin of less than 5% compared to the numerical solution.
The numerical solutions together with the original Koenig ansatz and the new approximations are shown in Fig. 5. Now we are going to derive the "integrated" equation for the cloud ice mixing ratio tendency: We may interchange the derivative with integration, invoke the product rule and arrive at: The second integral vanishes, since ∂m/∂t ≡ 0 (because m must be interpreted here a co-ordinate in m-space, and the coordinate system is, of course, fixed

Fig. 5.
Growth rates dm/dt of single ice crystals vs. crystal mass.Results of the numerical integration (red) and approximations using the original Koenig ansatz (green) and the Koenig ansatz with our corrections (blue) for various temperatures (T =230/220/210/200 K) and fixed pressure (p=300 hPa).
cast into another form when we make use of the "continuity equation" in mass-space, which reads Inserting this into the first integral above, we find

Partial integration yields
where the integrated part in the square brackets vanishes, because f (m) vanishes at infinity, and at the lower boundary m=0.Finally, we arrive at the following simple expression: Here, dm dt can be interpreted as the "advection velocity" in the "mass"-space due to crystal growth; hence, we can insert the modified Koenig approximation from above.Let us first treat the case where the largest crystals are still smaller than m l (i.e.f (m) ≈ 0) for m>m l ).This gives The original Koenig approximation results into an ice water mass rate of aµ b [m].Additional numerical integrations and our simulations showed that we can approximate the integral in Eq. ( 42) so that for sufficiently small crystals the tendency for q c can be cast into the following form: where χ≈20.Hence, we simply need another correction factor.It turns out that this correction will only have an impact for low temperatures and high vertical velocities (see below).
For large crystals, the ventilation correction becomes important.These crystals alone give the following contribution: which can be computed using known expressions for truncated moments (Jawitz, 2004).Finally, we arrive at the following expression: Note that dq c dt is expressed using the mean mass (first moment) and the moment of order b of the mass distribution.This formulation makes it possible to use any kind of mass distribution for which analytical expressions for the moments are known.
Growth of ice crystals does not affect their number concentration (as long as there is no aggregation), but sublimation does when the smallest ice crystals sublimate completely while larger crystals only loose some mass.We parameterise this effect in a simple way.Let the ice mass mixing ratio in a certain grid box be reduced by a fraction f q during one time step.Then we assume that the corresponding fractional reduction of number concentration is given by f N =f α q , with α>1.This relation implies that when a small mass fraction sublimates, this is mainly due to shrinking of big crystals; when a large mass fraction vanishes, then also a large fraction of crystal number concentration must vanish.In the limiting cases, f q =0 or f q =1, we also have f N =0 or f N =1, respectively, as it should be.From numerical studies we found that α=1.1 produces plausible results.This value is in agreement with Harrington et al (1995), who derived from numerical simulations a range 1≤α≤1.5.

Sedimentation of ice crystals
From our parameterisation of the terminal velocity of a single ice crystal we derive terminal velocities for the mass and number concentrations, respectively.To this end let us consider the sedimentation fluxes of mass and number concentrations: We use these definitions to define mass and number weighted terminal velocities, such that F m =ρ q c v m , F n =ρ N c v n .Hence: The terminal velocity of one single ice crystal is a function of its mass, viz.v(m)=γ (m) • m δ(m) .For simplicity, let us first assume γ (m) = γ 0 , δ(m)=δ 0 to be constants.This allows to simply express the mass and number weighted terminal velocities via the moments of f (m): Actually, the coefficients in the mass vs. fall speed relation are piecewise constant, hence we can express the integrals by using truncated moments: where µ k (m l , m u ) denotes the k th truncated moment with boundaries m l <m<m u .Here, m 0 =0, m 4 =∞, the values m 1 , m 2 , m 3 are given in Table 2.For the lognormal distribution the truncated moments can be expressed as follows (Jawitz, 2004): with the transformation For the calculation of the integral we have to evaluate the error function  6. Mean terminal velocities v m , v n for r 0 =3 and the approximations by using the expressions of general moments as in Eq. ( 50) as piecewise constant functions.
In Fig. 6 the mean terminal velocities v m , v n for a typical value of r 0 =3 (i.e.σ m =2.852) are shown.
The treatment using truncated moments is computationally expensive.Therefore we approximate the analytical solutions by using general moments for different mass intervals of the curve, i.e.
The approximations are also shown in Fig. 6.Evidently, the approximation works rather well.Certain mathematical relations between moments (Lyapounov's inequality) guarantee that the ratio v m /v n >1 for all r 0 >1 which implies that the mass mixing ratio sediments faster than the number concentration, or in other terms, that large ice crystals fall faster than smaller ones.
The ratio v m /v n depends strongly on the width of the ice crystal mass distribution (parameter r 0 ).For too large values of r 0 a kind of decoupling of the variables N c , q c could be observed: The ice crystal number density concentrates in the upper layers of a cloud, while the cloud ice falls downwards very fast.This leads to unphysical behaviour, e.g. to extraordinary large ice crystals in the lower cloud levels and fall streaks.For realistic values of r 0 in the range 1<r 0 ≤4 this does not occur.
In the EULAG model sedimentation is treated as a 1-D advection in vertical direction.Several advection schemes can be used: a simple implicit scheme or an explicit scheme (as used for the 3-D) with various orders and non-oscillatory option.In our studies we will usually use the explicit scheme of 2nd order (MPDATA) without the non-oscillatory option.P. Spichtinger and K. Gierens: Modelling cirrus clouds -Part 1a

Microphysical time step
As we will see later in the validation section, the parameterisations can be sensitive to the prescribed time step.For box model calculations this is not a problem; the computational effort is so small, that we can easily use a very small time step.However, for the application of the microphysics scheme in a 2D/3D framework of the EULAG model the computational effort is much higher.For solving this problem, we implemented a so-called microphysical time step.Here, the "dynamical" time step, i.e. the time step for the dynamics in the EULAG model is split into a number of subtime steps, i.e. t mp = t/n mp .There is an additional loop for each vertical column of n mp sub time steps; here, only the microphysics (nucleation, growth/sublimation, sedimentation) is calculated.Then, the summed forcing terms from the sub-time steps are used for the next dynamical time steps.We use the microphysical time step adaptively: First, the forcings for a normal dynamical time step are calculated.If homogeneous nucleation takes place within this time step or the change in relative humidity wrt ice is larger than 1% due to deposition/sublimation, then the time splitting is used.We will see later in the validation of the nucleation parameterisation, which number of sub-time steps is appropriate.

Validation and discussion
The microphysics scheme is validated using two different types of comparisons.First, we use our microphysics in a box model setting to compare the ice crystal number densities that are nucleated homogeneously under a large variety of conditions with values derived from a box model with detailed ice microphysics scheme (Kärcher and Lohmann, 2002a,b).Additionally, some sensitivity studies are carried out.Second, we simulate a well-documented case of arctic cirrostratus and compare our results with 1-D simulations (Lin et al., 2005;Kärcher, 2005).Kärcher and Lohmann (2002a) used a detailed box model (particle tracking, highly resolved aerosol size distribution) for testing an analytically derived relationship between the maximum possible ice crystal number density, formed by homogeneous nucleation, and vertical velocity.In our first step of validating the model, we carry out box model runs for the same conditions to compare ice crystal number densities with the results of Kärcher and Lohmann (2002a).The code of Kärcher and Lohmann (2002a) has been successfully applied to measurements of homogeneous nucleation in the large cloud chamber AIDA (Haag et al., 2003a) and to field measurements during INCA (Interhemispheric differences in cirrus properties from anthropogenic emissions, Gayet et al., 2006).So, our comparison with results of Kärcher and Lohmann (2002a) can be viewed as an indirect comparison with the AIDA and INCA results.

Comparison with detailed box model calculations
We use the following setup: We allow only homogeneous freezing and prescribe a constant vertical velocity in the range w=0.01−5.0m s −1 .The initial conditions (p, T , q v ) are adapted such that temperature and pressure reach prescribed values p±5 hPa, T ±0.1 K at the beginning of the homogeneous freezing event.
Initially we assume a very high aerosol number density so that nucleation is not constrained by the number of available nuclei.Later we will investigate cases with realistic background aerosol conditions, obtained from observations (e.g.Minikin et al., 2003).Within the box model framework we will also investigate the impact of pressure on the number density of ice crystals produced in the nucleation event.

Idealised simulations
In our first series of experiments we use only one class of (homogeneously formed) ice.The number density of the background aerosol is set to the very large value of N a =10 000 cm −3 , such that nucleation cannot exhaust the background aerosol.Therefore we may safely neglect the shift in the mean mass of the aerosol size distribution.We choose a geometric standard deviation of σ m =2.85 (r 0 =3) for the lognormal distribution of the ice crystal mass.Additional sensitivity tests have been performed with r 0 in the range 2≤r 0 ≤4; these resulted in only slight variations compared to r 0 =3.For a certain vertical velocity w we choose the time step such that the nucleation event is resolved: t=(0.05m/w) (B.Kärcher, personal communication).The number of time steps is fixed at n t =12 000.
We choose the same conditions (p, T ) as in Kärcher and Lohmann (2002a), i.e.
T =196, 216, 235 K and p=200 hPa, and a mean radius for the (dry) H 2 SO 4 aerosol of r m =25 nm.The reported values of temperature and pressure, respectively, in all boxmodel simulations refer always to the onset of homogeneous freezing.For investigating the effect of the width of the aerosol size distribution we vary the geometrical standard deviation: σ r =1.3, 1.6, 1.9.
The results are shown in Fig. 7 in comparison with the values of Kärcher and Lohmann (2002a,b).The agreement between our results and those of the more detailed model is quite satisfying, in particular for the two higher temperatures.At the lowest temperature our model still produces results similar to those of Kärcher and Lohmann (2002a,b), but only up to w=10 cm s −1 .Generally we find from additional simulations (not shown) that the maximal vertical velocity at which our model produces reliable results decreases with decreasing temperature.At higher vertical velocities our model starts to underestimate the number of aerosols freezing, and the underestimation grows with the width of the aerosol size distribution.We believe that these errors have two sources.
The first error source is the assumption of equilibrium implicit in the Koehler equation and in the Koop  Kärcher and Lohmann (2002a).The panels show the impact of the width of the aerosol distribution: top: σ r =1.3, middle: σ r =1.6, bottom: σ r =1.9.Fig. 8. Effect of a nucleation event during a time step of t=0.1 s at T =200 K and a relative humidity wrt ice RH i=161% for a quite broad aerosol distribution (r m =50 nm, σ r =1.9):The H 2 O/H 2 SO 4 droplet distribution before the nucleation event is indicated by the red line points, the green line indicates the distribution of freshly nucleated ice crystals and the blue line describes the H 2 O/H 2 SO 4 droplet distribution after the nucleation event.
parameterisation.Both the Koehler equation and the Koop parameterisation are based on the assumption of equilibrium, i.e. they are only strictly applicable when rates of water uptake on aerosol droplets are fast compared to the rate of change of saturation ratio.For high vertical velocities, especially in the low temperature range, this condition is not fulfilled.Hence, in a fast cooling environment the solution droplets are actually smaller than computed with the Koehler equation.The initial nucleation rate (when the threshold supersaturation is reached) is larger when the aerosol sizes are overestimated, first because of their increased volume, and second because of the overestimated water activity.However, the final number of ice crystals is the integral of J hom V d over time.Since the nucleation rate is initially overestimated, and since these crystals are too big, the initial crystal growth is overestimated as well, and so the initial depletion of excess vapour is overestimated.This can have two effects: First, the actual maximum supersaturation beyond the threshold is not reached in our model, hence maximum nucleation rates are underestimated.Second, the duration of the nucleation event is clipped.Overall, the time integral of J hom V d becomes underestimated, which is the effect that we see.
The solution of this problem requires introduction of prognostic equations for the aerosol dynamics, which currently is not part of our scheme.
The second error source works in a similar way, but it arises as an artifact of the assumptions required in a bulk model.When nucleation starts, it transforms the largest aerosol droplets into ice, cutting off the right tail of the droplet size distribution (see Fig. 8).For the simulations, we assume a fixed width for the aerosol distribution of σ r =1.4.Here, the impact of geometric mean size of the aerosol distribution is shown: r m =12.5 nm (black), r m =25.0 nm (cyan), r m =50.0 nm (magenta).
In the next time step this cutoff is forgotten, however, because the bulk model always assumes the same type of size distribution.This has the effect that at each time step during nucleation (except the first) the size of freezing droplets is overestimated.The newly formed crystals come out too large, consume too much excess vapour and quench nucleation prematurely, so that eventually too few crystals are produced.We can try to mitigate this error by reducing the mean aerosol size during the nucleation phase (switched off for the idealised experiments), but this does not always help, because the error depends in a non-linear way on the time step (simply because the mass distribution is a non-linear function).Later we will see that the mentioned problems become unimportant in cases with realistic aerosol background concentrations, so that fortunately it turns out unnecessary to invent a much more complicated correction for the mean aerosol size.
With the same setup as before we now study the effect of changing the geometric mean radius of the "dry" sulphuric acid: r m =12.5, 25.0, 50.0 nm.We hold the width constant at σ =1.4.Fig. 9 shows the results.
At the two higher temperatures there is hardly an effect of the aerosol core mass distribution.But at the lowest temperature and at high uplift rates, there is an effect.Our model follows the detailed microphysical results of Kärcher and Lohmann (2002a,b) the best when the acid mass in the droplets is smallest.In this case equilibrium between the solution droplets and their environment is easier to maintain compared to cases with larger acid fractions; hence use of the Koehler equations is still justified at strong uplift.We Fig. 10.Box model calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity and on the temperature.For the simulations, we assume a fixed width for the aerosol size distribution of σ r =1.4 and a geometric mean radius of r m =25 nm.Here, the impact of pressure on the formation of ice crystals is shown using values of p=200/300/400 hPa.have already stated above that at (water) saturation the equilibrium droplet radius r d is proportional to the square root of the solute mass, i.e. proportional to r 3/2 , which shows that generally droplets are smaller with smaller solute mass under otherwise identical conditions.The relaxation time (the time needed to equilibrate with changed ambient saturation ratio) of a droplet increases with its size because the number of water molecules that have to be transferred between the vapour and the droplets increases with droplet size.A simple remedy of the problem encountered above could therefore be to choose smaller values of r m or σ r than in corresponding spectral microphysics simulations.As we have seen, this choice gives better results at low temperatures and high uplift velocities while it has merely a weak effect at other conditions.Such an approach is justified as long as it is rather the ice than the aerosol that is the focus of the studies.
Finally we study the pressure dependence of the number concentration of ice crystals freezing by homogeneous nucleation.We expect an influence via the pressure dependence of the diffusion coefficient.Also results of Hoyle et al. (2005) point at the possibility that changes in pressure can change the amount of ice crystals formed in a homogeneous freezing event up to one order of magnitude.We select σ r =1.4,r m =25 nm as a setup that gave reasonably consistent results with Kärcher and Lohmann (2002a) who, however, did not investigate the pressure dependence.We choose temperatures in the range T =200, 215, 230 K and pressures in the range p=200, 300, 400 hPa.The results are presented in Fig. 10.
We can see clearly that the ambient pressure has an impact on the ice crystal number densities produced during the cooling experiment.The diffusion constant depends on pressure, i.e.D v ∝p −1 (see Eq. 31); thus diffusional growth rates decrease with increasing ambient pressure due to decreasing mean free path of the water molecules in air.Slower growth at higher pressures implies that supersaturation can stay longer above the nucleation threshold, hence more crystals nucleate.Varying the pressure from 200 to 400 hPa leads to an increase in the ice crystal number concentrations by a factor of 4-5.

Realistic background conditions
Realistic background aerosol concentrations are much lower than in the idealised simulations, for instance in the range n a =N a ρ∼100−300 cm −3 (Minikin et al., 2003).Hence, we repeat the simulations of Sect.4.1.1 with a realistic background aerosol density of n a =300 cm −3 .In this case we also use the correction for the mean aerosol mass, i.e. the aerosol size distribution is now shifted to smaller masses after a nucleation event.Under realistic conditions the number of ice crystals that form in homogenous nucleation events can be constrained by the available aerosol.
Figure 11 shows the impact of the width of the aerosol size distribution.We see that it hardly has an effect on the number of crystals formed, even at the lowest temperature and the highest uplift velocities (except for a very broad distribution with σ r =1.9).This is in sharp contrast to the idealised cases where σ r had a much larger influence at the lowest temperature.We conclude that it is here the constraint by the available aerosol rather than the equilibrium assumptions in the Koehler and Koop theories that leads to lower crystal numbers than in the Kärcher and Lohmann (2002a) studies.This implies that the problems encountered above are less important in reality than in the idealised situations, which means that our parameterisation can be used at low temperatures and high uplift velocities.
We also repeated the simulations with different pressures with the realistic aerosol distribution.While we still find the pressure effect at low uplift speeds, the effect gets weak or even vanishes at high uplift speeds due to the constraints posed by the number of available aerosol particles (see Fig. 12).

Time step issues
In all previous box model simulations we have adapted the time step such that the homogeneous nucleation event can be resolved in time.However, for more expensive 2-D/3-D simulations one would like to use longer time steps.We have tested how the model behaves with fixed time steps of t=0.1, 0.5, 1.0, 2.0 s, respectively, both in the idealised and realistic cases from above.Figure 13 shows the results.
Obviously, problems appear with too long time steps at high uplift speeds.
While there are no problems even with a time step of 1 s at small vertical velocities (w≤20−50 cm s −1 ), there are strong deviations from For the simulations, we assume a fixed geometric mean size of the aerosol distribution (r m =25.0 nm).Here, the impact of different width of the aerosol distribution is shown (black: σ r =1.3, cyan: σ r =1.6, magenta: σ r =1.9).In contrast to the simulations shown in Fig. 7, we use realistic background concentrations for the sulphuric acid aerosol (n c =300cm −3 ).This leads to less ice crystals for low temperatures and/or high vertical velocities due to a limited reservoir of solution droplets.Box model calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity and on the temperature.For the simulations, we assume a fixed width for the aerosol size distribution of σ r =1.4 and a geometric mean radius of r m =25 nm.Here, the impact of pressure on the formation of ice crystals is shown using values of p=200/300/400 hPa.In contrast to the simulations showed in Fig. 10, we use realistic background concentrations for the sulphuric acid aerosol (n c =300cm −3 ).This leads to less ice crystals for low temperatures and/or high vertical velocities due to a limited reservoir of solution droplets.the previously shown cases at higher vertical velocities (w>50 cm s −1 ) where the non-linear behaviour of both, the  For the simulations, we assume a fixed width for the aerosol distribution of σ r =1.4 and a geometric mean radius of r m =25 nm.Here, the impact of the prescribed (fixed) time step (dt=0.1/0.5/1.0/2.0) for the simulations is shown, compared to simulations using an adapted time step (simulation "ref", red line).For the simulations shown at the top panel, the background aerosol reservoir is large, while for the simulations at the bottom panel, we used a realistic background aerosol concentration (n c =300cm −3 ).nucleation process and the depositional growth introduce large deviations from the reference cases.The deviations are more severe for higher than colder temperatures.Crystal growth proceeds faster at warmer than at colder temperatures, which means that the duration of the nucleation pulse increases with decreasing temperature.When the nucleation event is not resolved, its duration will be overestimated and too much crystals will form.This effect is evidently the more severe the shorter is the nucleation pulse; i.e. the largest error occurs at the warmest temperature considered.

Formation and evolution of an arctic cirrostratus
In this section we compare the performance of our bulk microphysics scheme with the spectrally resolving schemes of Lin et al. (2005) and Kärcher (2005) for the case of an arctic cirrostratus triggered by a constant vertical updraught.

Setup
We use the following setup for our simulations: The whole 2D model domain (0≤x≤6.3km, 2≤z≤11 km) is lifted up adiabatically with a constant updraught velocity of w=0.05 m s −1 as described in Kärcher (2005).This is equivalent to a constant cooling of the background profile T e with a rate of dT /dt=dT /dz•dz/dt=−g/c p •w =−0.000489K/s.The same cooling rate was used in the box model simulations in Sect.4.1.The cooling is adiabatic (i.e.θ e is constant), and is continued for a total simulation time of t s =7 h.In Fig. 14 the initial profiles for the simulations are shown.
We use a horizontal resolution of x=100 m with a horizontal extension of 6.3 km, cyclic boundary conditions in x-direction, a vertical resolution of z=10 m and a dynamical time step of t=1 s.According to our discussion in Sect.4.1.3,there is no need of a small microphysical time step for the moderate vertical updraught in this case.For the background aerosol (H 2 SO 4 ) we use a number density of n a =N a ρ=300 cm −3 with geometric standard deviation σ r =1.4 and geometric mean radius of r m =25 nm for the lognormal distribution, as these values gave good results in Sect.4.1.

Results
For our comparisons we mostly refer to the simulation by Kärcher (2005), because he also parameterised homogeneous nucleation according to Koop et al. (2000).
In Fig. 15 the time evolution of the variables relative humidity wrt ice, ice crystal mass and number concentrations, resp., are shown.
The first nucleation event occurs at t≈60 min.The supersaturation peak of about 154% RHi triggers homogeneous nucleation.Within a few minutes a large amount of ice crystals (N c ρ∼100 L −1 ) is formed.Because of the high supersaturation the ice crystals can grow quickly and deplete a fraction of the water vapour, which reduces the relative humidity, see Fig. 16.Ice crystals grow and soon start to fall.Therefore the peak of high supersaturation at the top of the ISSR is influenced very weakly by the depletion of the water vapour.The peak is permanently maintained for the whole simulation time and is a permanent source for homogeneous nucleation at the top of the ISSR.
The combination of crystal growth and sedimentation causes two effects: On the one hand, the supersaturation is reduced by crystal growth such that the relative humidity cannot reach the threshold for homogeneous nucleation in layer would everywhere lead to nucleation, then the crystals would grow until the excess vapour would be consumed over the whole depth of the layer.Sedimentation obviously plays a crucial role for the development and the structure of the simulated cirrus cloud and for the maintenance of supersaturation within the cloud.The first nucleation event forms a large number of ice crystals.Many crystals fall downwards, resulting in a downward moving peak of high ice crystal number densities in fig.17.These results agree qualitatively well with those of Kärcher (2005), yet there are differences in details.For instance, the in-cloud supersaturation is higher in our simulation than in Kärcher's.This is simply a consequence of different assumptions of crystal shape in the two codes.While we use hexagonal columns with an aspect ratio r a > 1 for crystal lengths L ≥ 7.42µm, Kärcher (2005)  ical crystals (r a = 1) for crystal lengths up to L ∼ 25µm; spherical crystals always grow faster than columns.Hence high supersaturation can be maintained for a longer period in our simulation than in Kärcher's.Second, our treatment of sedimentation is more diffusive than that of Kärcher (2005) who uses a particle approach with the advection scheme by Walcek (2000).Numerical diffusion in the double-moment scheme (see Wacker and Seifert, 2001) leads to smoothing of vertical gradients in ice crystal number concentrations and cloud ice mixing ratio, such that peak values are smaller than in Kärcher (2005).
In spite of these differences in details we find that our bulk microphysics scheme is able to well reproduce the main features of the arctic cirrostratus case, namely: -High supersaturation at the top of the cloud with continuously ongoing homogeneous nucleation; the lower part of the cloud.On the other hand, the falling ice crystals formed at the top of the cloud are the only sink for the water vapour.Although the continuous homogeneous nucleation events permanently form new ice crystals, these are spread vertically over the whole cloud depth resulting in relatively low number densities.Thus, inside the cloud, ice supersaturation is maintained.If there were no sedimentation or if it were very weak, the cooling in the supersaturated layer would everywhere lead to nucleation, then the crystals would grow until the excess vapour would be consumed over the whole depth of the layer.Sedimentation obviously plays a crucial role for the development and the structure of the simulated cirrus cloud and for the maintenance of supersaturation within the cloud.The first nucleation event forms a large number of ice crystals.Many crystals fall downwards, resulting in a downward moving peak of high ice crystal number densities in Fig. 17  These results agree qualitatively well with those of Kärcher (2005), yet there are differences in details.For instance, the in-cloud supersaturation is higher in our simulation than in Kärcher's.This is simply a consequence of different assumptions of crystal shape in the two codes.While we use hexagonal columns with an aspect ratio r a >1 for crystal lengths L≥7.42 µm, Kärcher (2005) assumes spherical crystals (r a =1) for crystal lengths up to L∼25 µm; spherical crystals always grow faster than columns.Hence high supersaturation can be maintained for a longer period in our simulation than in Kärcher's.Second, our treatment of sedimentation is more diffusive than that of Kärcher (2005) who uses a particle approach with the advection scheme by Walcek (2000).Numerical diffusion in the double-moment scheme (see Wacker and Seifert, 2001) leads to smoothing of vertical gradients in ice crystal number concentrations and cloud ice mixing ratio, such that peak values are smaller than in Kärcher (2005).In spite of these differences in details we find that our bulk microphysics scheme is able to well reproduce the main features of the arctic cirrostratus case, namely: -High supersaturation at the top of the cloud with continuously ongoing homogeneous nucleation; -Significant supersaturation within the cloud for the whole simulation time; -Downward moving peak of high ice crystal number concentrations formed at the first nucleation event.Lin et al. (2005) noted that high vertical resolution is needed for the reproduction of these features.We have checked this issue using various resolutions in the range 5≤z≤50 m.It turned out that the number of ice crystals that are produced is underestimated with too coarse vertical resolution, which occurs when the continuous nucleation source at the cloud top is unresolved.This renders the supersaturation inside the cloud too high.Results change only marginally between simulations with resolutions of z=5 m and z=10 m.Therefore, we used z=10 m.
Regarding the occurrence of ice supersaturation inside cirrus clouds we want to remark that this feature is not limited to the very cold temperature range in the tropics (an example is presented by Gao et al., 2004).Also in the temperature range of midlatitude and arctic cirrus clouds ice supersaturation inside the clouds can be found frequently, see for instance the results of Ovarlez et al. (2002); Comstock et al. (2004); Krämer et al. (2008) and also the recent SPARC newsletter about the supersaturation issue (Peter et al., 2008).

Conclusions
We have described a new bulk microphysics parameterisation for simulation of cirrus clouds on the cloud resolving scale.In the two-moment scheme the processes nucleation (homogeneous and heterogeneous), deposition growth/sublimation and sedimentation of ice crystals are implemented.An arbitrary number of ice classes can be treated that are discriminated by their formation mechanism.Each ice class is connected to an aerosol type that freezes into the respective ice class.We tried to formulate process rates as functions of general moments of the assumed underlying ice mass distribution in order to keep the model flexible enough that various kinds of mass distribution can be chosen.Currently we use log-normal distributions for the aerosol and the ice masses.The new parameterisation was implemented into a simple box model and into the anelastic, non-hydrostatic model EULAG.
In a first validation step we compared the ice crystal number densities generated during box model simulations with a steady and constant updraught with results from a box model that was equipped with a more detailed microphysics scheme (Kärcher and Lohmann, 2002a).The agreement of our results with the results from the detailed model was very good.Additionally we used the box model to study the impact of the ambient pressure on the generated ice crystal number concentrations.Also, for realistic background aerosol conditions simulations were carried out.In this case we could find that the background aerosol acts as a limiting factor, i.e. the produced ice crystal number density is strongly reduced in the low temperature range in combination with high vertical velocities (consistent with results of Kay and Wood, 2008).This partly offsets the underestimation of ice crystal number densities of our parameterisation.Additionally, the pressure dependence of the parameterisation was investigated.
In a second validation step we simulated the case of an arctic cirrostratus using EULAG as a 1-D column model first.Also in this case the agreement with a much more detailed microphysical model (Kärcher, 2005) was very good.
These studies led to the following conclusions: -The model is validated against models with detailed microphysics and produces reliable results; -Sedimentation is of utmost importance in the evolution of the cloud structure and the in-cloud humidity field; -Persistent in-cloud supersaturation is found in our simulations consistent with Kärcher (2005).
Finally, we conclude that the model is suitable to represent cirrus clouds in cloud-resolving simulations.In the second part of our study (Spichtinger and Gierens, 2009), we will investigate the validation case of an arctic cirrostratus carefully in terms of changing the large scale dynamics (i.e.updraught velocity) and in terms of possible changes due to 2-D effects introduced by temperature fluctuations and wind shear.In future applications the model will be used for investigating the effects of the competition of different nucleation mechanisms (Spichtinger and Gierens, 2008) and for the multiscale problem of the impact of mesoscale dynamics on the formation and evolution of cirrus clouds (as in Spichtinger and Dörnbrack, 2006).

Fig. 1 .
Fig. 1.Length (L) and diameter (D) of hexagonal ice crystals as functions of the ice crystal mass m.

Fig. 2 .
Fig. 2. Aspect ratio r a of hexagonal ice crystals versus ice crystal mass m (upper panel) and length L (lower panel), respectively.

PFig. 3 .
Fig. 3. Terminal velocity v(m) of a single ice crystal as function of the crystal mass m.

Fig. 4 .
Fig. 4. Coefficients a(p, T ) (top) and b(p, T ) (bottom) for Koenig's ansatz dm/dt≈a(p, T )m b(p,T ) and the approximation of the numerical solution of Eq. (30) for different pressures.

Fig. 7 .
Fig. 7. Box model calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity.The calculations for different temperatures are indicated by the colours (red: T =235 K, green: T =216 K, blue: T =196 K).Line points indicate calculations with our box model, points denotes values fromKärcher and Lohmann (2002a).The panels show the impact of the width of the aerosol distribution: top: σ r =1.3, middle: σ r =1.6, bottom: σ r =1.9.

PFig. 9 .
Fig.9.Box calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity and on the temperature.For comparisons, the values from detailed microphysics calculations inKärcher and Lohmann (2002a) are shown.For the simulations, we assume a fixed width for the aerosol distribution of σ r =1.4.Here, the impact of geometric mean size of the aerosol distribution is shown: r m =12.5 nm (black), r m =25.0 nm (cyan), r m =50.0 nm (magenta).

Fig. 11 .
Fig.11.Box model calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity and on the temperature.For comparisons, the values from detailed microphysics calculations inKärcher and Lohmann (2002a) are shown.For the simulations, we assume a fixed geometric mean size of the aerosol distribution (r m =25.0 nm).Here, the impact of different width of the aerosol distribution is shown (black: σ r =1.3, cyan: σ r =1.6, magenta: σ r =1.9).In contrast to the simulations shown in Fig.7, we use realistic background concentrations for the sulphuric acid aerosol (n c =300cm −3 ).This leads to less ice crystals for low temperatures and/or high vertical velocities due to a limited reservoir of solution droplets.
Fig. 12.Box model calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity and on the temperature.For the simulations, we assume a fixed width for the aerosol size distribution of σ r =1.4 and a geometric mean radius of r m =25 nm.Here, the impact of pressure on the formation of ice crystals is shown using values of p=200/300/400 hPa.In contrast to the simulations showed in Fig.10, we use realistic background concentrations for the sulphuric acid aerosol (n c =300cm −3 ).This leads to less ice crystals for low temperatures and/or high vertical velocities due to a limited reservoir of solution droplets.

Fig. 13 .
Fig.13.Box model calculations for the homogeneously formed ice crystal number densities depending on the vertical velocity and on the temperature.For comparison, the values from detailed microphysics calculations inKärcher and Lohmann (2002a) are shown.For the simulations, we assume a fixed width for the aerosol distribution of σ r =1.4 and a geometric mean radius of r m =25 nm.Here, the impact of the prescribed (fixed) time step (dt=0.1/0.5/1.0/2.0) for the simulations is shown, compared to simulations using an adapted time step (simulation "ref", red line).For the simulations shown at the top panel, the background aerosol reservoir is large, while for the simulations at the bottom panel, we used a realistic background aerosol concentration (n c =300cm −3 ).

Fig. 16 .
Fig. 16.Temporal evolution of the vertical RHi profiles for simulation time 60 ≤ t ≤ 90 min.The notch of the profile due to ice crystal growth is represented clearly, while the supersaturation at the cloud top remains.

Fig. 16 .
Fig. 16.Temporal evolution of the vertical RHi profiles for simulation time 60≤t≤90 min.The notch of the profile due to ice crystal growth is represented clearly, while the supersaturation at the cloud top remains.

PFig. 17 .
Fig. 17.Downward moving peak of high ice crystal number densities formed by the first homogeneous nucleation event at t≈60 min in the reference simulation (i.e.constant uplift of w=0.05 m s −1 ).At the top of the cloud there is continuous nucleation, indicated by the high peak values of ice crystal number density.
of spheroid a(p, T ) factor for Koenig's approximation (abbrev.by a) a w water activity in the solution a i w water activity in the solution in equilib.with ice A Ice crystal surface A = √ a 2 − b 2 b short half axis of spheroid b(p, T ) exponent for Koenig's approximation (abbrev.by b) C Capacitance factor c p heat capacity for constant pressure c(T , p) correction factor for sedimentation

Table 1 .
Values for α, β in the general mass-length relation; here, m t =2.146 • 10 −13 kg denotes the transition between aspect ratio r a =1 and r a >1, this ice crystal mass is equivalent to a ice crystal length of L t =7.416µm

Table 3 .
Geometric standard deviations σ m and σ L , depending on the ratio r 0 .The transition mass m t =2.146 • 10 −13 kg is equivalent to a length of L t =7.416µm.
(p, T ) mass for appr. in Eq. (38) (abbrev.m 0 ) m i i=1, 2, 3 transition masses for terminal velocity m m geometric mean mass in lognormal distribution m pre predominant mass m t transition mass between spherical crystals and columns M w molecular weight of water vapour N a aerosol number concentration per kg dry air n a aerosol number density (n a =N a ρ) per m 3 dry air N c ice crystal number concentration per kg dry air n c ice crystal number density (n c =N c ρ) per m 3 dry air N P r