MAID: a model to simulate UT/LS aerosols and ice clouds

The comprehensive model MAID (model for aerosol and ice dynamics) was developed to simulate condensation and freezing in aerosol particles residing in the UT/LS (upper troposphere/lower stratosphere). The exact balancing of trace gas components is a particular emphasis of MAID. MAID is applied to and verified by experiments in the aerosol chamber AIDA, and, moreover, it is adapted to Lagrangian atmospheric cirrus cloud simulations. Here, the model is introduced, and as an example for model applications the significant influence of homogeneous or heterogeneous freezing on ice cloud microphysics and the water and nitric acid partitioning in cirrus clouds is shown.


Nomenclature
concentration of condensable trace gas i in the carrier gas x = ∞: far away from particle, x = sur: at particle surface C 0 initial particle number concentration D i diffusion coefficient of condensable trace gas i d g mean geometric diameter K (r i , r j ) coagulation frequency between particles of sizes r i , r j m i mass of component i in particle of radius ṙ m i mass growth rate for component i m nv,k non-volatile mass in size bin k M i molecular mass of component i n ( R, r, t) concentration of particles of size r at location R at time t n (r k , t) concentration of particles in size bin k at time t p i partial pressure of condensable trace gas i in the carrier gas S(r k , t) source rate of particles of size r k at time t t time T : absolute temperature v i average molecular velocity of condensable trace gas i x (...) mass fraction of component (...) α Dep (r k , t) sum of first-order deposition coefficients for particles of size r k at time t δ ik Kronecker symbol σ g geometric standard deviation

Introduction
Cirrus clouds play an important role in the Earth's radiation budget. Their frequency of occurrence and cirrus microphysics is initially determined by the freezing mechanism. For instance, decreasing the freezing threshold of ice clouds will lead to an increase in high cloud cover, which in turn feeds back to the radiation balance of the atmosphere (Gettelman and Kinnison 2007). Subsequent supersaturations inside cirrus influence the cloud's microphysics and also the vertical redistribution of water vapour through the sedimentation of ice crystals. Changes of the cirrus freezing threshold and in-cloud supersaturations can be caused by alterations in the aerosol population. To improve our understanding of the aerosol influence on ice cloud nucleation, microphysics and supersaturations, we developed the process model MAID (model for aerosol and ice dynamics). MAID is especially suitable to study the evolution of ice clouds, because it exactly balances the implemented trace gas components and is based on elementary aerosol physics without restrictions regarding size distributions and particle composition. The model is applied to and verified by ice nucleation experiments at the aerosol chamber AIDA, and, moreover, it is adapted to Lagrangian atmospheric cirrus cloud simulations. Here, the model is introduced, and examples for the different applications are given.

Model description
Since simulation of single particles can be performed only in idealized situations and small systems containing a few thousand particles only, aerosol systems have to be simulated by a density function. This approximation is the basic concept of the general dynamic equation of aerosols (GDE, e.g. Gelbard and Seinfeld 1979) which describes the behaviour of an ensemble of aerosol particles influenced by external forces and by interactions with the surrounding carrier gas and with each other. In fact, the GDE can be considered as a type of continuity equation for the density or concentration of aerosol particles n ( R, r, t). The terms in the differential equation for the time derivative of n( R, r, t) represent the influences due to changes within the particle size distribution (e.g. coagulation, condensation, etc) and transport processes due to carrier gas convection and the particles' motions caused by external forces (gravity, electrical, and phoretic forces) including the 'pseudoforce' diffusional transport due to density gradients. The actual importance and formulation of the different terms depend on the conditions the GDE has to be applied for. Therefore, MAID takes into account the physical and chemical processes considered as important for the simulation of UT/LS aerosols such as condensation and evaporation, uptake of trace gases, freezing of droplets, ice formation initiated by solid particles, deposition and coagulation of particles as an option.
Two versions of the model are available at the moment, a 'chamber version' and a 'Lagrangian version'. The two versions mainly differ in the way the thermodynamic boundary conditions such as temperature, pressure, partial pressures of trace gases and so on are treated. The details will be discussed below, but the main focus of this paper is directed to the Lagrangian applications.

Aerosol and ice treatment
2.1.1. Aerosol particles. The density function mentioned above depends on the size, composition, and the shape of the particles as well as on their location. The space dependence can and has to be replaced by appropriate boundary conditions if it can be assumed that the simulated box is homogeneously mixed due to convective turbulent mixing. This assumption is valid for many situations either for closed systems (e.g. chambers) or for small volumes in atmospheric applications, e.g. control volumes in a Lagrange calculation. Motions of particles caused by interaction with the surrounding fluid or by external forces (e.g. gravity, electrical fields) can only indirectly be treated within the framework of such a spaceindependent solution either by deposition coefficients which are the results of special solutions for the boundary layers in cases of closed containers or by transport coefficients in other cases such as Lagrangian air parcels. Complete spacedependent solutions can only be achieved by direct coupling of the GDE to fluid dynamics simulations, which considerably increases the computational effort.
Numerically this space-independent density function can either be approximated by a sectional or bin-wise representation or by special functions, e.g. log-normal distributions. In the first case the GDE can be transformed into a set of coupled first-order differential equations (Bunz and Dlugi 1991), whereas in the second case a set of differential equations for the moments of the size distribution can be formed (e.g. Ackermann et al 1998). On the basis of the moments the distribution parameters (e.g. variance, geometric mean diameter, etc) of the size distribution function can be calculated. It should be mentioned that not only a single size distribution function but also a sum of distribution functions can be used, increasing the flexibility and accuracy of the method but also the number of differential equations and the necessary computation time. The size distribution keeps this shape with distribution parameters varying as a function of time, whereas on the basis of the sectional method any arbitrary distributions can be simulated and be used as particle source functions. This additional flexibility is the reason why this method was preferred in MAID. Due to the fact that the sizes of the particles to be simulated differ by orders of magnitude, the logarithm of radius or volume is used as the distribution parameter. To minimize numerical diffusion the volume of a low or (even better) a non-volatile aerosol component is used instead of the total volume. In most cases H 2 SO 4 is taken due to its extreme low vapour pressure at our temperatures. The masses m nv in the following equations represent only this non-volatile fraction of the particles, whereas the total mass as well as the particle radius calculated on the basis of the total mass vary as a function the bin number (a moving grid method). To calculate the behaviour of the volatile fractions, separate differential equations have to be solved, which are discussed in section 2.1.2. Therefore, based on such a sectional representation the GDE can be transformed into a system of coupled first-order differential equations, and is given by The different terms represent (from left to right) the explicit time derivative of the concentration of particles characterized by size r k at time t, production rate of such particles by coagulation of smaller particles, loss of particles due to coagulation with other particles, any arbitrary production term of such particles, and first-order loss term (discussed below). Due to the fact that r k represents only the non-volatile fraction of the aerosol, no term for condensation appears in this form of the equation. This effect is treated in accompanying equations which have simultaneously to be solved, as already mentioned above. The coefficients β are necessary to interpolate the particles being formed due to coagulation into the sectional representation chosen to represent the particle size distribution. They are formulated in a way that mass as well as particle number balances are guaranteed.
Depending on the ambient conditions the particles may contain a number of different components such as organics (Murphy et al 2006) and inorganics such as H 2 SO 4 , H 2 O, HNO 3 , HBr, HCl, and non-soluble fractions such as soot or mineral dust. At the moment due to lack of detailed information only solid cores and inorganics are implemented in MAID, but due to the flexible structure of the program other components can be implemented quite easily. As discussed in more detail in the next sections, the uptake (condensation) of volatile aerosol components can be calculated as well as the ice nucleation phenomena as a function of thermodynamic conditions and particle composition. The total particle volume is computed by summing up the volatile contributions to the non-volatile core. The exact particle size is then determined for each bin. Liquid aerosols and the frozen particles are represented by two distinct distributions, interacting on the one hand via condensation and evaporation processes and on the other hand via coagulation if this process is taken into account (high number concentration cases). If a frozen particle and a liquid particle coagulate it is assumed that the resulting particle is frozen as well. Of course, it is checked if this new particle represents a thermodynamically stable state according to its composition and the surrounding temperature and relative humidity.
The coefficient α can be used to take deposition on surfaces, exchange between adjacent volume elements, and similar processes into account. The selection of the physical effects necessary to evaluate α depends on the boundary conditions. In the cases presented in this work only the change of the particle concentration due to changes of volume by pressure and temperature variations and by the adiabatic pumping of the AIDA chamber is considered, since in the case of the Lagrangian calculations a steady state is simulated and in the case of the AIDA calculations the time is quite short in comparison to the usual deposition processes (sedimentation, diffusion, phoretic effects). But is can easily be implemented if necessary.

Condensation and evaporation.
The mass flux of a volatile component to or from a particle is proportional to the difference between the partial pressure (or the concentration) of the corresponding trace gas and its vapour pressure (or saturation concentration) at the particle surface times a kinetic pre-factor depending on the particle size, the mean free path length of the gas, diffusivity, and the accommodation coefficient of the condensing species. In the model the respective formulation of Dahneke (1983) is implemented: for i = 1, NC (= number of components in the particle). (2) The equation describes the growth of a particle characterized by radius r and composition m 1 . . . m NC . Of course, only equations for the volatile components (mainly H 2 O and HNO 3 in our cases) have to be taken into account. The kinetic factor is based on the solution of the partial differential equation describing stationary diffusion around a body, in this case a spherical aerosol particle. This assumption is strictly valid as long as droplets are considered and approximately valid even for ice crystals with a compact structure. The boundary conditions applied are the bulk concentration of the trace gas far away (in fact infinity) from the particle and the composition and temperature-dependent equilibrium concentration on the particle surface. For particles, which are smaller or of the order of the mean free path length, these boundary conditions have to be modified to account for non-continuum effects. As the result, an additional factor depending on the mean free path length of the carrier gas, the particle size, and the accommodation coefficient α interpolates between the continuum regime (large particles) and the free-molecule regime (small particles). Whereas the transport in the continuum regime is determined by diffusion only, accommodation is important for the free-molecule regime dominated by collisions between the particle surface and trace gas molecules. The influence of α increases with decreasing particle size and decreasing α. Parametric calculations for particles in the size range between 10 nm and 1 μm showed that the accommodation coefficient α becomes important if α is below 0.1 (Carver and Harrington 2006). Measurements of α show that it is in the order of 0.5 for the condensation of H 2 O on electrolytic solutions (Gershenzon et al 2002), whereas values for the deposition on ice surfaces range from 0.001 to almost 1, depending on the temperature and experimental conditions (e.g. Magee et al 2006 and references quoted therein). This could have an effect on the degrees of saturation which can be maintained in ice phase clouds due to the delayed uptake of water by the ice crystals. But simulations performed with MAID for the CR-AVE 2006 field campaign (Costa Rica-Aura Validation Experiment, Gensch et al 2008) indicate that the accommodation coefficient for water vapour on cirrus cloud particles is close to one at temperatures around 180-185 K.
For temperatures above −20 • C the growth rate formulations (in particular for H 2 O) have to be modified to include the combined effect of heat and mass transfer due to the higher mass flux rates and the corresponding higher release rates of latent heat due to the condensation process. In all cases the influence of the condensation and evaporation on the H 2 O partial pressure has to be quantified, together with other loss or gain terms such as transport to and from walls or across the surfaces of the volume element, since the H 2 O partial pressure and the saturation level derived from it are the most important parameters to determine particle growth as well as freezing phenomena (section 2.1.2). With this method an excellent mass balance of any trace component can be achieved as a function of time.
The vapour pressure at the particle surface consists of the product of the vapour pressure of the solution as a function of composition and temperature and the Kelvin term taking into account the increase of vapour pressure due to the curvature of the particle surface in the case of a convex shape of the particle surface. In our model, parameterizations of the vapour pressures of solutions of H 2 O/H 2 SO 4 (Tabazadeh et al 1997), H 2 O/H 2 SO 4 /HNO 3 /HCl/HBr (Luo et al 1995), and H 2 O/H 2 SO 4 /HNO 3 /NH 3 (Lin and Tabazadeh 2001) as functions of temperature and composition are implemented.
It should be mentioned that chemical reactions taking place in and on the particles have to be treated in a similar way as differential equations modifying the mass and chemical form of certain particle components. If necessary, they can easily be implemented into the model.

Ice formation.
Ice formation can occur either homogeneously as a spontaneous process in liquid aerosols supersaturated with respect to ice or heterogeneously on a solid kernel by deposition or immersion freezing. The nucleation rate J (cm −3 s −1 ) for homogeneous freezing is calculated using the model of Koop et al (2000), which takes advantage of the fact that the freezing temperature correlates very nicely to the water activity of the droplets (figure 1 in Koop et al 2000).
It should be mentioned that the activity of the droplets is calculated in MAID on the basis of the actual particle composition for any individual particle size. Therefore, the growth kinetics influences the actual water activity which has to be taken as input to calculate the freezing probability. In all cases characterized by fast changes of temperature and H 2 O partial pressure the water activity for large drops differs from the gas phase activity since the relative growth rate (the inverse of the characteristic condensation time) scales with r −2 , as can be easily seen from equation (2). This is the case for fast rising air parcels in the atmosphere or during an expansion experiment in the AIDA chamber. The water content and the corresponding water activity of large drops will be below the equilibrium value due to their slower relative mass uptake rate slowing down their freezing probability considerably in comparison to calculations based on the gas phase water activity (figure 1). The figure shows a somehow arbitrary snapshot of the water activities of sulfuric acid droplets as a function of the size and the corresponding freezing rates (nucleation rate multiplied by the droplet volume) calculated on the basis of the Koop model during a simulated expansion in the aerosol chamber AIDA with a cooling rate of about 2.5 K min −1 . The water saturation changes (S w = 89% or S ice = 150% at the moment of the snapshot) quickly during such scenarios. In real nature, such conditions can be found in deep convection clouds. The actual position of the maximum of the freezing rate depends mainly on the cooling rate or the rate of increase of the water vapour saturation.
For the saturation water vapour pressures the descriptions of Marti and Mauersberger (1993) and Tabazadeh et al (1997) are used for ice and for supercooled liquid water, respectively.
For the heterogeneous freezing process the 'shifted activity' parameterization derived by Kärcher and Lohmann (2003) is implemented in the model, accounting for the reduced critical freezing saturation ratio (freezing thresholds) of heterogeneously freezing aerosol particles. The freezing thresholds of heterogeneous ice formation are derived from AIDA ice nucleation measurements of soot particles coated with sulfuric acid (Möhler et al 2005), which we believe to be important ice nuclei in the UT region. Using a pre-defined threshold for the heterogeneous ice nuclei number density, ice may nucleate until all these insoluble cores are consumed. At that moment, heterogeneous ice nucleation is suppressed. If the air parcel is cooled down even further, two alternatives for the later cloud evolution are open: either the water vapour depletion by deposition on the ice can compensate the increase of relative humidity causing a substantial growth of the few existing ice crystals, or the relative humidity increases continuously and exceeds the threshold for a subsequent homogeneous ice nucleation, producing high number densities of ice crystals.
The actual freezing rates of any droplets in a bin are calculated on the basis of the corresponding nucleation rates. The fractions of the aerosol components are conserved at each time step; they only pass from the liquid to the ice phase. The retention fraction of the volatile species described is considered in MAID to be near unity. Consequently, the molecules dissolved in the liquid phase are entirely retained in the ice particles upon freezing. The ice crystals are supposed to sublimate if the ambient temperature is higher than their dew point. In that case, any contained soluble trace species will evaporate at a rate consistent to the mass transfer law, depending on the partial vapour pressures in the system.
The ice crystal growth due to water vapour deposition is treated in the same manner as for liquid aerosols. Nitric acid uptake does not affect the crystal size because of very low molar ratios, HNO 3 : H 2 O < 10 −4 , but this process is described in MAID using the HNO 3 trapping parameterization derived by Kärcher and Voigt (2006).

The chamber version of MAID
The treatment of the aerosol behaviour is identical for the two versions. Therefore, calculations performed for chamber experiments can be directly compared to atmospheric conditions. The two versions differ mainly in the way in which the thermodynamic boundary conditions are handled. This version allows the user to calculate the chamber thermodynamics as a function of time, based on the initial conditions and the chamber specific boundary conditions such as expansion rate, surface temperature, and degree of ice coverage. The expansions are used to simulate atmospheric situations characterized by decreasing temperature and increasing saturation levels and to investigate particle growth and freezing. In MAID the heat fluxes caused by the adiabatic cooling and by the convective heating of the chamber volume by the walls are balanced, and the temperature is calculated accordingly. The partial pressures of trace gases (mainly H 2 O) are determined on the basis of their initial conditions, sources, and sinks such as particle condensation/evaporation and the convective transport of sublimating ice from the chamber walls which are warmer than the chamber atmosphere during the expansion. If measured data of these thermodynamic variables are available the program can be forced to use these data instead of the calculated ones.

The Lagrangian atmospheric version of MAID
The Lagrangian version of MAID is distinguished from the AIDA version by using the temperature and pressure from atmospheric air parcel trajectories calculated from ECMWF wind fields. The partial pressures of the trace gases treated in MAID can either be calculated by balancing their fluxes or can be forced by tabulated data as is the case for the other thermodynamic variables. At the moment only those trace gas fluxes caused by condensation/evaporation of the particles are taken into account as well as the changes due to variations of atmospheric pressure and temperature. But, if the necessary data are available, transfer rates across the surfaces of the simulated volume element could easily be implemented.

Model applications
3.1. AIDA chamber experiment As mentioned above, particle growth and possible subsequent freezing is initiated by adiabatic pumping of the 84 m 3 chamber. In figure 2, the water partial pressure is shown for such an expansion experiment with H 2 SO 4 particles (initial values: C 0 = 12 000 P cm −3 , d g = 0.3 μm, σ g = 1.5, x(H 2 SO 4 ) = 34 wt%) in the time immediately prior to ice formation. Comparisons between measurements and calculations demonstrate that the model is able to balance water vapour sources and sinks (mainly the dilution due to the expansion and the uptake by the growing particles) and is able to identify wall effects which are difficult to quantify by measurements alone. There are always some small inhomogeneities concerning the wall temperature and the corresponding degree of ice coverage. Water saturation increases due to the temperature drop (approximately 4-7 K) during the expansion. The saturation levels and temperature at the onset of ice formation can be determined on the basis of measurements and model calculation (figure 3). The figure shows the onset saturation levels at different temperatures for of a number of experiments with sulfuric acid aerosols. Due to temporal and spatial fluctuations in the chamber on the one hand and the assumption of spatial homogeneity in the model on the other hand, there are some deviations between modelled and measured results (figure 4). But the onset and fraction of activated droplets are described quite well by the model, having in mind the difficult experimental conditions. Some improvements can be expected in the future and MAID is directly coupled to a detailed modelling of the fluid dynamics and the temperature and trace gas distribution in the chamber.
Due to the low vapour pressure of ice in comparison to liquid water the ice crystals take up water vapour very quickly, and the saturation decreases after ice formation, stopping further ice formation rapidly.

Atmospheric ice clouds
As an example of atmospheric ice cloud simulations with MAID, we perform a sensitivity study along an atmospheric trajectory calculated using the wind fields provided by ECMWF. Mesoscale temperature fluctuations with a peak to peak amplitude of 1.6 K (Gary 2006) are superimposed on the synoptic scale temperature history. The resulting temperature along the trajectory is plotted in figure 5, panel (a). Two possible cirrus scenarios are simulated along the trajectory: one where solely homogeneous freezing of aerosol particles forms ice crystals and a second where both heterogeneous and homogeneous freezing is allowed. The model simulations are initialized based on aircraft measurements performed at the endpoint of the trajectories during the field campaign POLSTAR in 1997 (Meilinger et al 1999, Krämer et al 2003. The observed log-normal aerosol size distribution showed a total particle number of N ptcl = 156 cm −3 , a mean size of R ptcl = 0.05 μm, and a width of sigma ptcl = 2.25. 10.3 ppmv of total (gas phase + ice crystals) water was measured and used as model input. Sulfuric acid was calculated to be 0.5 ppbv and nitric acid as 0.1445 ppbv for the observed conditions. The number of heterogeneous ice nuclei (IN) is set to 0.01 cm −3 following Blake and Kato (1995), who sampled black carbon soot aerosols from the upper troposphere and lower stratosphere and derived latitudinal number densities. The number density N ice of the arising ice crystals is shown in panel (b) of figure 5 (blue: homogeneous; red: heterogeneous case), and the relative humidity with respect to ice Rh ice is shown in panel (c). The first ice crystals appear in the heterogeneously formed cirrus (het) shortly after the lower heterogeneous freezing threshold was reached (red vertical line in figure 5). A second, homogeneous ice forming event occurs for the heterogeneous case almost at the same time when the purely homogeneously formed cirrus (hom) appears, namely when the homogeneous freezing threshold is reached due to further cooling. Rh ice (het) does not significantly drop until the second nucleation step but fluctuates, showing a negative correlation to the temperature fluctuations. Rh ice (hom) and then Rh ice (het) decreases after this ice nucleation event, fluctuating around saturation. The amplitude of Rh ice (het) fluctuations is a little larger than that of Rh ice (hom). The parameter controlling Rh ice is the integral particle size N ice R ice , which is smaller in the heterogeneous case (panel (c)), so the water vapour consumption by the ice crystals is lower and more water remains in the gas phase.
The uptake of nitric acid on liquid aerosol HNO 3,liq (panel (e)) is sensitive to the evolution of Rh ice . Before the homogeneous ice nucleation begins, up to 80% nitric acid accommodates from the gas phase into the liquid particles. After the homogeneously formed cirrus appears, a significant part of HNO 3,liq escapes from the interstitial particles again in the gas phase. Due to the sharp drop of Rh ice (hom) under 100%, the release of HNO 3 from the interstitial particles is more accentuated in the homogeneous case, compared to the heterogeneous case. The HNO 3,gas is trapped efficiently by the ice crystals. The uptake of nitric acid on ice HNO 3,ice (panel (f)) is determined by the integral particle size N ice R ice , similarly to Rh ice : HNO 3 is slowly trapped by the ice crystals after the first heterogeneous ice nucleation and rapidly increases after the homogeneous event, with a higher uptake of the available HNO 3 in the solely homogeneous case.  Results of Lagrangian simulations and comparisons to the corresponding field measurements are presented and discussed in more detail in the accompanying paper of (Gensch et al 2008).

Conclusions
The model MAID is applied to both laboratory experiments at the aerosol chamber AIDA and atmospheric trajectories. It is demonstrated that the model is able to exactly reproduce the water fluxes in the aerosol chamber AIDA and that the fluxes from the walls are not negligible. From a sensitivity study of atmospheric ice cloud formation it is evident that the upper tropospheric aerosol particle composition significantly impacts the ice cloud microphysics as well as the partitioning of water and nitric acid. Due to the lower freezing thresholds of heterogeneous freezing regarding temperature as well as water vapour saturation levels, ice clouds can develop already at higher temperatures and lower saturations in comparison to scenarios characterized by homogeneous ice nucleation only. Further, heterogeneously formed ice clouds can contain less water and nitric acid in the ice phase (see also the companion paper of Gensch et al 2008). Thus, injecting IN such as soot or mineral dust particles into the upper troposphere is able to influence the cirrus coverage and the radiative properties of the clouds, as well as the vertical redistribution of water and nitric acid by sedimentation and convective transport of ice crystals.