Adaptive method of lines for multi-component aerosol condensational growth and cloud droplet activation

Introduction Conclusions References


Introduction
Cloud droplets form on cloud condensation nuclei, CCN (Aitken, 1875, and the works of his contemporaries).CCN are the atmospheric aerosol that contribute to the heterogeneous nucleation of the water vapour.As a result, the microphysical properties of clouds (i.e. the size spectrum of cloud droplets) strongly depend on physicochemical properties of aerosol particles (see e.g., the review of McFiggans et al., 2006).Owing to the fact that the radiative properties of clouds and their ability to precipitate are both dependant on characteristics of the size spectrum of cloud droplets (see e.g., the review of Stevens and Feingold, 2009), the description of the activation process is of importance to the studies on aerosol-cloud-climate interactions (e.g., Kulmala et al., 2009).This paper discusses a numerical method for studying the kinetics of aerosol growth by vapour condensation in the context of cloud-droplet activation.Presented method is of potential applicability in other branches of aerosol/hydrosol research.Figures

Back Close
Full Aerosol growth by condensation of water vapour is conceptualized as a twofold diffusion process.First, there is the diffusion of water molecules towards the particle, and the condensation of water vapour on the drop surface.Concurrently, the latent-heat release triggers the diffusion of heat from the droplet, hence the system is akin to the wet-bulb thermometer described by Maxwell (1878).The chemical composition of the growing particle and its curvature determine the vapour pressure over the droplet surface (K öhler, 1936, and his earlier works), the driving factor of the diffusion process.Additional complexities are introduced when considering an ensemble of particles of different sizes (Howell, 1949;Kraus and Smith, 1949), with different shapes and soluble material content (Kornfeld, 1970), nuclei of different solubility (Fitzgerald, 1974), and finally nuclei composed of different chemical species (Lee and Pruppacher, 1977) all competing simultaneously for the available water vapour.Yet another level of detail is reached when considering the feedback with the dynamics of the air (turbulence, entrainment and mixing), the interactions between the individual droplets (i.e.collision/coalescence) or more sophisticated description of cloud chemistry, all being however distant from the scope of this study.Numerical solutions of the problem had been repeatedly challenged using the fastest computers in the early years of electronic computing, starting with the studies of Mordy (1959) performed using the Swedish BESK computer, and the study of Neiburger and Chien (1960) carried out on the American SWAC.
The key complexity of the problem, from the mathematical standpoint, is the necessity to solve the partial differential equation describing the evolution of aerosol size spectrum.Analytical solutions to the problem (or semi-analytical estimates) can be derived only under significant approximations (see e.g., Squires, 1952;Twomey, 1959;Gelbard, 1990;Khvorostyanov and Curry, 2009, and references therein).Numerical solutions based on finite-difference approach using a fixed-grid representation of the size spectrum may be used (also called Eulerian or fixed-bin approaches).Such methods are characterised by computational difficulties (related to the numerical diffusion and the need for employing positive-defined algorithms), and the complexity of the Figures

Back Close
Full description of chemical properties of aerosol (see e.g., Chap. 5 in Williams and Loyalka, 1991).The partial differential equation may be reduced to a system of ordinary differential equations (ODE) by spatial discretization of the initial size spectrum into sections with constant-in-time number of particles.The boundaries of these sections are then allowed to evolve with time, hence the name moving-sectional technique/Lagrangian approach/full-moving structure (Jacobson, 1999, Sect. 14.5 therein) or method of characteristics (Gelbard, 1990), or method of lines (Debry et al., 2007).This approach was introduced by Howell (1949), and has since been used in modelling CCN activation (e.g., Mordy, 1959;Neiburger and Chien, 1960;Fitzgerald, 1974;H änel, 1987;Feingold and Heymsfield, 1992;Nenes et al., 2001b;Grabowski et al., 2010).Such formulation of the model eliminates the computational difficulties inherent to numerical treatment of partial differential equations, and allows accurate description of the evolution of aerosol composition because particle properties such as mass of solute are retained throughout the computation.
The moving-sectional technique does pose computational difficulties because of the ODE system stiffness (see e.g., the discussion by Gelbard, 1990).The physical process of competition for the water vapour between the droplets of different sizes and chemical composition can be correctly resolved in a numerical solution only under sufficiently short time-steps.The approximations inherent to the bin-representation of the aerosol size spectrum introduce significant uncertainties to the solution if low-resolution discretization is used (Korhonen et al., 2005).The issues of stiffness of the system and of appropriate time-step choice are addressed herein by using an adaptive Gear-type ODE solver.The issue of sensitivity of the results to the bin-number choice, being one of the factors causing discrepancies among results of different models (Kreidenweis et al., 2003) is the very focus of this study.The problem is tackled by introducing a simple adaptive size-spectrum discretization refinement.The resultant solution procedure is analogous to what is referred to as the adaptive method of lines (Wouwer et al., 2001).Introduction

Conclusions References
Tables Figures

Back Close
Full The following section of the paper covers model formulation and implementation.Sample results are discussed in the third section.First, a set of model runs without adaptivity is discussed revealing the sensibility of the results to parameters controlling the discretization of the initial aerosol size spectrum (Sects.3.1 and 3.2).Second, the properties and efficacy of the adaptive size-spectrum discretization method are presented (Sect.3.3).Third, a comparison of the results with the predictions of the classic Twomey's (1959) analytical formulae is given (Sect.3.4).The fourth section concludes the paper with a brief outlook on the possible applications of the discussed method within different domains of atmospheric research.

Model formulation and implementation
The model provides a description of the cloud-drop formation process at the initial stage of cloud-formation triggered by vertical air motion (i.e. in convective or orographic clouds).The model utilizes the method of lines (MOL) for numerically solving the socalled dynamic equation of aerosol growth by condensation.A list of symbols used in the text is given in Table 3.

Air-parcel system and its governing equations
The model describes an air parcel with the temperature T , pressure p, and specific humidity q v being vertically displaced with constant velocity w.The parcel does not exchange heat or mass with the environment and hence its heat budget (neglecting terms related to mass and heat content of droplets) is expressed as:

Conclusions References
Tables Figures

Back Close
Full respectively; c pv and c pd depict specific heat capacities at constant pressure for water vapour and dry air, respectively.The parcel's pressure adapts instantaneously to the environmental value defined through the hydrostatic equilibrium.It follows that the pressure changes are driven by the vertical velocity of the parcel: where g is the gravitational acceleration.Physical characteristics of the aerosol suspended in the parcel are described by a set of probability density functions (PDFs) defining the size spectra of dry aerosol, each PDF corresponding to different chemical composition.The dry spectra n w are expressed in terms of particle number per unit mass of moist air, and are referred to as specific concentrations.The model formulation is aimed at describing the growth of aerosol at high humidity, and corresponding high relative water contents within particles of the wet spectra.The total mass of particles is assumed to be proportional to the third moment of the wet PDFs: where w dr w is the total specific concentration for a given chemical component.For each wet spectrum a corresponding droplet temperature field T Water is considered to be present in liquid and vapour phases only.The diffusion of both water vapour and heat around the droplets is treated assuming stationary state Introduction

Conclusions References
Tables Figures

Back Close
Full with time-varying boundary conditions what results in the following form of the Fick's first law and Fourier's law (cf.e.g., Chap.7 in Rogers and Yau, 1989): where M=4πr 3 w ρ l /3 is the mass of a drop of density ρ l , ρ v is the density of air in the parcel, D and K are diffusion coefficients discussed below.The two diffusion equations are coupled via a heat budget equation: where c l is the specific heat capacity of liquid water.Equations ( 4) and ( 5) evaluated with constant (or temperature-dependant) diffusion coefficients D and K represent the so-called continuum r égime in which the efficiency of diffusion is limited solely by the diffusivity of water vapour in air.Continuum r égime is applicable for particles of sizes larger than the free path of water molecules in air (i.e.r w 0.1 µm).Evolution of size of much smaller droplets is described using the molecular r égime, and the transition between the two r égimes may be approximately allowed for by introducing a drop radius-dependency to the diffusion coefficient (Langmuir, 1944, discussion of Eq. 52 therein).Here, we apply the transition r égime correction factors of Fuchs and Sutugin (1970) in the form suggested for cloud modelling by Laaksonen et al. (2005):

Conclusions References
Tables Figures

Back Close
Full where D 0 is the diffusivity of water vapour in air, K 0 is the thermal conductivity of air, λ D is the mean free path of water vapour in air, and λ K is a length scale for heat transfer in air (the ratio of a mean free path to a characteristic length is called the Knudsen number).The two last are defined as: and Loyalka, 1991, Eqs. 6.6 and 6.33 therein).
The condensation/evapouration of water on/from the aerosol is the only source of changes in specific humidity q v and thus: N [c]  (11) where the (1−q v ) term stems from the employment of specific quantities.The density of water vapour at drop surface ρ v | drop surface in Eq. ( 4) is assumed to be equal to the equilibrium (saturation) value taking into account the Gibbs-Kelvin and Raoult effects related to surface curvature and drop composition, respectively (the socalled K öhler curve): where p v∞ is the equilibrium vapour pressure over plane surface of pure water, σ is the surface tension coefficient, and the water activity a (Raoult term) is parameterized in Introduction

Conclusions References
Tables Figures

Back Close
Full accordance with experimental data using the single-parameter representation of Petters and Kreidenweis (2007, Eq. 6 therein): where κ [c] is the single parameter representing chemical composition in spectrum c.
The dry aerosol is assumed to be entirely soluble.A constant surface-tension coefficient σ=0.072J/m 2 is used (value for pure water at 298 K) in accordance with the parameterization of Petters and Kreidenweis (2007).
The specific heat capacities of dry air c pd , and water vapour c pv are assumed constant (perfect gas), hence the latent heat of vapourisation is given by: and the equilibrium pressure over a plane surface of pure water (solution of the Clausius-Clapeyron equation) is given by: where T 0 , p 0 and l v0 correspond to the triple-point values, and T denotes here the temperature for which the formulae are evaluated.The specific heat capacity of liquid water c l is assumed constant.

MOL's ODE system
As outlined in the first section, the mathematical formulation of the model constitutes a system of ODEs which approximate the partial differential equation for describing the evolution of aerosol size spectrum.This approach is based on an approximation of the Introduction

Conclusions References
Tables Figures

Back Close
Full  [c,b] (per unit mass of air), and variable-intime left and right boundaries r wr (t).The third moment of the particle size spectrum is thus approximated using: This allows expressing the time derivative of the total mass of wet aerosol without employing partial differentials using: An assumption that r is inherent in this formulation (cf. the discussion and a related proof of particle ordering principle in Gelbard, 1990).The same discretization is applied to dry radii r d , and particle temperatures T w introducing r While still constant from the standpoint of the ODE solver, N b , N [c,b] and r The system of ODEs defines the derivatives of all the variables with respect to time as: The drop-growth is described here using separate equations for drop-temperature and drop-radius evolution (see e.g., Vesala et al., 1997), as opposed to a single explicit equation (the so-called Mason equation, Howell, 1949;Tsuji, 1950;Mason, 1957).

Initial conditions
The initial dry-aerosol size spectra and their chemical composition, initial temperature, pressure, humidity and the vertical velocity profile are the key parameters for each model run.Introduction

Conclusions References
Tables Figures

Back Close
Full PDFs of tropospheric aerosol are generally well represented using a trimodal lognormal distribution (Whitby, 1978).The initial spectra of dry aerosol are thus defined for each model run using (Williams and Loyalka, 1991, Eq. 1.22 therein): where m is the geometric standard deviation of a mode, and r [c] m is the mode radius (all defined separately for each chemical component).N m is generally reported as concentration (defined in terms of a unit volume, as opposed to unit mass as it is in case of specific concentration n d ) hence the division by the air density ρ in Eq. ( 20).The bins are initially equally spaced in logarithm of radius over the 1 nm-100 µm range what ensures accurate representation of the log-normal shape of the initial spectra.
The initial spectra of wet aerosol are calculated by assuming an equilibrium state defined by: The initial values of T , p and q v are chosen to describe the air properties slightly beneath the cloud base, i.e. at high relative humidity (RH) but below saturation.High initial relative humidity assures consistency with the dilute-solution approximation inherent in the Raoult's law.Introduction

Conclusions References
Tables Figures

Back Close
Full

Numerical solution
Integration of Eq. ( 19) without further approximations requires a numerical algorithm capable of solving a stiff ODE system (the stiffness of the system is confirmed e.g. by the significant inefficiency of the Adams-Moulton method with functional iteration).
In addition, calculation of the initial wet spectra using an equilibrium condition defined by Eqs ( 21) is not possible analytically.The latter problem is solved by employing an iterative root-bracketing procedure (using the Brent-Dekker method, cf.Sect.2.6).The initial search interval for root-bracketing is defined using the dry radius for the lower bound, and the equilibrium value with the Gibbs-Kelvin effect neglected for the upper bound.
The ODE system is solved using the variable-order, variable-step backward differentiation formula (BDF) as implemented in the CVODE solver (Cohen and Hindmarsh, 1996, cf. Sect. 2.6 herein).In each time-step (numbered by i ) the ODE system dx/d t = F (x,t) defined by Eq. ( 19) is integrated by finding a root of f (x i ): using the Newton-Raphson iteration (here h varies between 1 and 5, and defines the order of the method at the current time-step, and α i ,k & β i are variable coefficients).
The iteration (with steps numbered by j ) is of the following form: The Jacobian matrix ∂F/∂x is updated at least once every twentieth time-step.The Jacobian is approximated by multiple evaluations of the right-hand side of Eq. ( 19) using the current, and slightly perturbed values of x.The CVODE solver adaptively chooses time-step according to the desired accuracy to be achieved in each step of integration (relative tolerance of 1×10 −8 and zero

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full absolute tolerance was used in the calculations presented here).In addition, the solver is instructed to shorten the time-step if an unphysical condition is detected during evaluation of the system right-hand side F (x,t) (e.g. if a positively-defined quantity has a negative value or if the left edge of a bin outruns its right edge).
The right-hand side of the ODE system is evaluated in the order presented in Eq. ( 19).The integration is performed until a desired time (altitude) is passed, and an interpolated set of values for this time (altitude) is calculated and returned.
After each successful integration step, a spectrum-refinement procedure is carried out, if necessary, and in that case the integration over the last time-step is repeated.The conditions for triggering the spectrum-refinement are defined and discussed in the following subsection.Each refinement causes a change of the length of x, which in turn results in reducing the current order h to unity, and re-evaluation of the Jacobian matrix.
A solution of the problem in context of cloud-microphysical studies is a set of properties of the cloud-droplet spectrum such as cloud-droplet number concentration (CDNC), liquid water content (LWC), average droplet radius ( r ), the standard deviation of the size spectrum (σ r ) and effective radius (r eff ).The k-th moment of the droplet spectrum is evaluated, taking into account that the left-most (and right-most) bins may fall only partially in the assumed cloud-droplet region, using: where

The adaptive size-spectrum discretization scheme
The key task of the adaptive discretization scheme is to control the accuracy of the discretization by splitting any overly wide bins into several smaller bins while the ODE system is integrated.Using the MOL nomenclature (Wouwer et al., 2001, Sect. 1.3 therein) it is thus an h-refinement of the spatial grid.Insertion of new bins takes place at discrete time levels only, namely when the ODE solver reaches its next internal time-step.Consequently, the grid adaptation procedure is decoupled from the ODE integration and as such the static griding is employed.The method utilizes the socalled full restart as the ODE solver is reinitialised after each refinement event, and restarts with the lowest order.The restart and the reduction of the method order is a consequence of changing the length of the model state vector x.
The criterion for triggering the division of a given bin is defined as follows.Bin splitting happens whenever the ratio of the logarithm of the width of a given bin to its initial value reaches a certain limit.The limit was set at 2 for most of the simulations discussed in the text (save for the case of urban aerosol where 1.8 was used, cf.Sect.3.3), meaning that no bin was allowed to double its width when viewed in logarithmic scale (the previous step of the solver is repeated each time new bins are added).Smaller limits cause the algorithm to carry out the refinement earlier in time and hence increase the computational effort.Larger limits delay the refinement and at the same time increase the uncertainty inherent to the interpolation performed to calculate the drop temperatures and dry/wet radii at the edges of the newly created bins.No refinement takes place when the drop concentration in a given bin is below a certain level hereinafter referred to as the tolerance.Introduction

Conclusions References
Tables Figures

Back Close
Full Any refinement event causes a full-restart of the ODE solver which, in turn, causes the solver to reduce the order, tighten its time-step, re-evaluate the Jacobian matrix, and consequently slow down the computations.For this reason, the number of bins into which a bin is split is chosen in a way ensuring no further splitting of any of the newly created bins will occur.The number of new bins added is thus chosen by dividing N [c,b]     by the tolerance.The new bins are spaced linearly in radius, and hence the specific concentration in each of the new bins is equal.The values of drop temperature and drop dry/wet radii at the newly introduced bin boundaries are linearly interpolated from the values at the boundaries of the bin being split.
The refinement criterion and the refinement procedure are local to a given bin and thus different bins may be split at different times during the integration.In principle, the presented approach replaces the single parameter of the discretization, the number of bins, with two new parameters: the initial number of bins and the tolerance.The former determines the uncertainty related to the discretization of the initial condition (initial dry radii).The later controls the uncertainty associated with the final spectrum, and has a physical meaning, being the highest number of droplets allowed to reside in a bin laying in between the two modes of the final aerosol spectrum (i.e.modes of activated and non-activated particles).

Implementation of the model
The program is implemented in an object-oriented manner in C++.The implementation is build upon the Boost.Units framework for dimensional analysis of the code at compilation time (Schabel and Watanabe, 2008) hence a physical unit is assigned to every variable/expression in the code with no performance penalty.The GNU Scientific Library (Galassi et al., 2009) is used for root-bracketing using the Brent-Dekker method (cf.Sect.2.4 herein and Sect.33.8 therein) and for root-polishing using the Introduction

Conclusions References
Tables Figures

Back Close
Full  et al. (2006) while two out of five models compared in Kreidenweis et al. (2003) and the model of Debry et al. (2007) were based on VODE.
The source code of the program is available as an electronic supplement of this paper and is released under the GNU General Public License.All model parameters are passed using UNIX-like command-line options.Toggling such settings of the model as transition r égime corrections, the Gibbs-Kelvin effect or the adaptivity do not require recompilation.Similarly, altering the numerical solution parameters such as the tolerances or switching from BDF to Adams-Moulton method for ODE integration is possible using command-line options.A web-based interface (written in HTML/PHP) allowing alteration of all model options, visualisation of single simulation results, as well as control of parallel runs of simulations with different input parameters is included.

Example run without adaptivity: sea-salt/sulfate competition in marine aerosol activation
Ghan et al. ( 1998) presented an analysis of relative importance of sea-salt and sulfate particles during CCN activation.The study was based on a series of simulations using a moving-sectional model of CCN activation.We use here a similar set-up for discussing an example model run without adaptivity and relatively low spectral resolution (45 bins over the 1 nm-100 µm radius range).
In the set-up the initial aerosol content of the air parcel is described by a sum of a single-mode log-normal distribution of ammonium sulfate aerosol, and a tri-modal log-normal distribution of sodium chloride aerosol.The parameters of all modes are listed in Table 1 and  A steady updraft w=0.25 m/s is driving the changes in the parcel's temperature, pressure, and relative humidity which are initially set at 280 K, 1000 hPa and 99%, respectively.The integration is performed for 250 s.
A summary of evolution of the aerosol spectra is presented in Fig. 1 (see figure caption for details).The shift towards larger sizes of the initial wet spectra (green lines in Fig. 1a) with respect to the initial dry spectra (red lines) reflects the growth of particles up to RH=99% calculated assuming equilibrium condition (cf.O'Dowd et al., 1997, Fig. 4 therein).Further growth of the aerosol is calculated resolving the diffusion kinetics, and leads to activation of CCN and forming cloud-droplets (peaks on right-hand side of the final blue spectra in Fig. 1b).The final spectra of both sulphate and sea-salt particles are distinctly bi-modal.The large-particle mode corresponds to cloud droplets, and is visible in both sulphate and sea-salt distribution.The smallparticle mode corresponds to unactivated particles.The largest droplets were formed on sea-salt aerosol.For both aerosol components, the region of spectra between the two modes of the final PDFs is represented with a single bin spanning over an order of magnitude in radius.This very section(s) is hereinafter referred to as the critical section(s) following Korhonen et al. (2005).
Figure 2 presents the time-evolution of particle sizes, particle temperatures (colour scale representing T w −T ), and the time/height profile of relative humidity.Wet radii corresponding to the two chemical components considered are plotted using circles and diamonds.The differences in wet radii for the same dry size (i.e. the differences arising from the differences in water activity for sea-salt and sulphate solutions) are most pronounced in the left-hand side of the plot (i.e.solution droplets of less than 10 nm in radius).Elsewhere, the two symbols (i.e.circles and diamond) overlap.The symbols are drawn every fiftieth time-step of the solver, and what follows the time-step adjustments made by the solver may be observed.
Four phases of the process can be distinguished (Howell, 1949).First, the particles grow slowly approximately until the relative humidity reaches 100%.From then on, the relative humidity does not fall below 100% and hence all subsequent phases occur at Introduction

Conclusions References
Tables Figures

Back Close
Full supersaturation with respect to the plane surface of pure water.The second phase is characterised by a rapid growth of particles in the 0.5-5 µm range, and continuous increase of the relative humidity until reaching a maximum value of 100.2%.In the third phase, occurring roughly from 120 s until 150 s, the rate of growth of particles stabilizes throughout the whole spectrum, with the exception of the region separating the unactivated and activated modes.There, one of the bin-boundaries undergoes a different evolution for sulphate and sea-salt aerosol, initially having equal volume of dry material (highlighted with the thick light-gray bifurcated line in the background in Fig. 2).The boundary of the very bin in the sulphate spectrum becomes the smallest one in the activated part of the spectrum, while the one of the sea-salt aerosol is the largest one in the unactivated part of the spectrum.The relative humidity decreases during the third phase.Finally, the fourth phase of the growth represents a quasiequilibrium state with activated droplets steadily growing at almost constant relative humidity.
The characteristics of the cloud-droplet spectra (spectrum moments calculated over the 1 µm-25 µm radius range) are presented in Fig.  and the integration process.As pointed out in Kreidenweis et al. (2003) the droplet number concentrations obtained with models relying on the moving-sectional technique are sensitive to the number of bins used for the discretization of the aerosol spectra.
To assess the degree to which the choice of the number of bins used to represent the initial spectrum influences the results, a set of simulations run with different number of bins is compared.The variability of the results due to changes in discretization parameter is compared with the variability associated with the aerosol physical and chemical properties, and the vertical air velocity.Four tri-modal aerosol distributions corresponding to properties of aerosol of different air-mass types are taken into account.Their parameters are listed in Table 2 and are based on the marine, clean continental, average background and urban cases reported in Whitby (1978).For each distribution, a nuclei mode, an accumulation mode, and a coarse mode is defined.The spectra are assumed to represent the sizes of dry particles.This is the set of cases used in Nenes et al. (2001b) for studying the influence of aerosol characteristics on cloud radiative properties using a moving-sectional model of CCN activation.
For each of the considered initial aerosol size distributions the calculations are carried out separately assuming pure ammonium sulphate and pure sodium chloride aerosol, and for three different updraft speeds, namely 0.25, 1 and 4 m/s.Calculations are performed until reaching 125 m.Each set-up (defined by one of the four sizespectra, one of the two chemical compositions, and one of the three vertical velocities) is tested with 54 different numbers of bins ranging from 30 to 300, with a step of 5.All other model parameters are defined as in the example run described in Sect.3.1.
Figure 4 presents the results from all 1296 model runs.Six quantities are compared, out of which five represent the final values reached at the altitude of 125 m above the initial level: the drop concentration CDNC, the mean radius of the droplets r , the standard deviation of the radius σ r , the effective radius r eff and the liquid water content LWC.The sixth quantity is the maximum supersaturation S max =max(RH−1) encountered during the model run, expressed in percent.It is evident that the model results

Conclusions References
Tables Figures

Back Close
Full are sensitive to the number of bins used to represent the initial spectra of aerosol.The maximum supersaturation is the only quantity not exhibiting significant variation.The variability associated with the changes in the number of bins is in some cases stronger that due to different chemical composition or updraft velocity (e.g.σ r in the marine case).Gradual increase of the number of bins does lead to the convergence of the result to an asymptotic value; however, arguably 300 bins is not a sufficient number in some cases.Clearly, the results obtained with less than 100-150 bins are characterised by significant uncertainty.
The observed range of values of CDNC is roughly in line with observations reported in Kreidenweis et al. (2003) and Korhonen et al. (2005).While the differences in the model formulations, and the differences in model set-ups rule out direct comparison of the results, it is worth to mention several of observations from the two mentioned studies.In both cases the focus was given to the drop number concentration.Kreidenweis et al. (2003, Par. 22 and Fig. 8 therein) compared results from four different models, running them with different number of bins ranging from 20 to 100.The bins were spaced equally in logarithm of radius and in all cases covered at least two decades of the 1 nm-10 µm size range.The variability of drop number concentrations was observed for all of the models.The discrepancies between the results obtained with the same model but with different number of bins reached ca.10%.Korhonen et al. (2005) compared results obtained using five different activation schemes run at low spectral resolution (10 bins within the 10 nm-1.5 µm range) with a high-resolution run in which 500 bins were used.The variations were analysed separately for marine, rural, and urban aerosol, and for a variety of vertical velocities (all below 1 m/s).Predictions of the most accurate low-resolution models differed from the high-resolution calculations by roughly, up to 1%, 10% and 50% for marine, rural, and urban aerosol, respectively (Korhonen et al., 2005, Fig. 3 therein).
The problem of the sensitivity of the moving-sectional model results to the spectrum discretization resolution was also discussed in Takeda and Kuba (1982, Sect. 2.5 therein) where the authors mention large influence of the number of particles in a

Conclusions References
Tables Figures

Back Close
Full critical section on the total number of activated cloud droplets, as well as broadening of the cloud droplet size spectrum with increasing number of bins used.Ghan et al.
(1998) used bins of equal number of particles, as opposed to the present logarithmic layout, effectively increasing the discretization resolution near the critical section(s).

Example results with adaptivity
The same series of simulations as presented in the previous section was repeated with adaptivity.The tolerances were set to min(0.5%•Nm , 50 µg −1 ) using N m of the accumulation mode (see Table 2) giving 0.3 µg −1 , 4 µg −1 , 11.5 µg −1 and 50 µg −1 for marine, clean continental, average background and urban cases, respectively.The results of multiple simulations with adaptivity are summarized in Fig. 5 intended for direct comparison with Fig. 4. The previously observed variability of the results with varying number of bins is visibly suppressed confirming the efficacy of the proposed method.The quality of the results (i.e. the stability with regard to the number of initial bins) obtained with low initial number of bins is significantly enhanced, in particular in the case of LWC.This suggests that introduction of the adaptivity helps the model to correctly resolve both droplet mass-and droplet number-related quantities, even though the model is formulated using particle number distributions.Both the effective radius r eff and clouddroplet number concentration CDNC do not exhibit virtually any variability (both being the very focus quantities of the studies on cloud-radiative properties, cf.e.g., Brenguier et al., 2000, and 2005) is related to the poor representation of the spectrum shape near the critical section(s); ii) the adaptive discretization scheme correctly identifies the critical section(s); iii) the proposed refinement strategy helps to limit the uncertainty of the results.
Figure 6 summarizes the relation between the initial and final number of bins, and thus provides a brief estimate of the memory/CPU footprint of the adaptive refinement scheme.The key observation is that the final number of bins corresponding to 1293 Introduction

Conclusions References
Tables Figures

Back Close
Full ca.50 bins on the abscissa (i.e.roughly the number of bins needed to obtain convergence with adaptivity) is, in general, lower than the number of bins needed to obtain convergence without adaptivity.Consequently, the method offers a reduction of computational effort needed to obtain results of a given precision.An in-detail comparison of one of the simulations presented in Figs. 4 and 5 is presented in Fig. 7a (without adaptivity) and 7b (with adaptivity).The vertical velocity was set at 1 m/s, the aerosol size spectrum was initialized with the average background distribution, and the composition was assumed to be pure ammonium sulphate.To enhance the readability of the figure by lowering the number of bins, the initial number of bins was set at 55, and the tolerance was decreased to 46 µg −1 .During the simulation with adaptivity 14 bins were added to the spectrum.There are two most evident differences among the plots in Fig. 7a and b.First, the final concentration calculated with adaptivity is about 15% higher, and consequently the predicted average drop radius is smaller.Second, the shape of the spectra at different heights is asymmetric for the case without adaptivity.
Enabling adaptivity, and thus increasing the resolution of the discretization in the very region of the spectrum where the drop-growth is most vigorous, allows the model to resolve the shape of the left-hand side slope of the droplet spectrum.This suggests that the asymmetry of the drop spectrum near cloud base, predicted by the model without adaptivity, is a numerical artifact.Finding support for this observation in experimental data is difficult as most of the cloud droplet spectrometers nowadays are optical instruments which, for physical reasons, cannot sense the smallest (i.e.<1.5 µm in radius) cloud droplets (see e.g., Fig. 10a in Lawson and Blyth, 1998, presenting the cloud spectra 50 m above cloud base as measured by three different instruments).
The ability of the model to faithfully resolve the spectrum shape at the very onset of droplet growth is of potential importance for the development of parameterizations of the activation process for use in cloud-resolving hydrodynamical models (see e.g., the discussion in Grabowski and Wang, 2009, where the authors underline the need for parameterizing the initial width of the spectrum of activated droplets).Introduction

Conclusions References
Tables Figures

Back Close
Full Comparison of the right-hand side panels of Fig. 7a and b suggests that the adaptivity contributes as well to a better representation of the transition-r égime corrections in the model.The increase of the number of bins near 100 nm results in probing the transition-r égime interpolation formulae at denser intervals.

Comparison with the Twomey's analytical formula
In his seminal paper of 1959, Twomey presented a derivation of an analytical solution of the CCN activation problem.His solution offers two robust expressions for upper bounds of maximum supersaturation and cloud drop concentration.The latter have since been widely used in cloud modelling (cf.e.g., the discussion of Grabowski et al., 2010, and references therein).
In the formulae of Twomey (1959) the supersaturation is expressed in terms of the dew-point elevation: where the dew-point temperature T d =p −1 vs (p v ) is the inverse function of the equilibrium vapour pressure, and p v = p v (p,q v ) is the partial pressure of water vapour.In the Twomey's solution, the upper bound of the maximum dew-point elevation is a function of vertical air velocity w, and two parameters describing the aerosol k T and c T : where CDNC is expressed in cm −3 .
Twomey (1959) has given three sets of example value pairs for k T and c T , including two for marine aerosol and one for remote continental aerosol.Twomey's relations are compared here with results of multiple simulations using the present model.Aerosol spectra from (Whitby, 1978) are used for defining the dry aerosol spectra.The marine and clean continental distributions defined in Table 2 are compared with Twomey's marine (A) (k T =1/3, c T =125 cm −3 ) and remote continental (k T =2/5, c T =2000 cm −3 ) cases, respectively.The model was initialized at temperature of 10 • C, pressure of 800 hPa, and relative humidity of 99.5% to closely resemble Twomey's set-up.The initial number of bins was set to 40, simulations were run with adaptivity, with the tolerance set at 5 µg −1 .As with previous set-ups, all calculations were done both for pure ammonium sulphate and pure sodium chloride for the purpose of comparing the degree to which the results are sensitive to chemical composition as opposed to other parameters.For both chemical species (i.e.NaCl and (NH 4 ) 2 SO 4 ) and both size spectra (i.e.marine and continental) the model was run at 25 different updraft velocities w between 0.01 and 5 m/s (with an increment of 0.2 m/s).In Fig. 8 the model results are compared with the values obtained using the Twomey's formulaefor the upper bounds of drop concentration and maximum dew-point elevation.The plotted values of drop concentration correspond to the concentrations at the altitude of 200 m above the initial altitude thus in the fourth stage of the activation process when the concentration is nearly constant.The maximum dew-point elevation corresponds to the difference between the dew-point and the parcel temperature at the timestep of the maximum relative humidity (supersaturation).The dew-point temperature is computed by iteratively (Steffensen method, cf.Sect 2.6) inverting Eq. ( 15).noteworthy; however, no further quantitative analysis seem feasible taking into account the degree of approximation employed in the Twomey's solution (including neglect of the Gibbs-Kelvin, Raoult and transition-r égime effects on the condensation process) and its upper-bound nature (see e.g., the discussion in Johnson, 1981) as well as the possible discrepancies between the aerosol described in Whitby's and Twomey's works

Summary and outlook
A simple yet robust air-parcel model of cloud drop formation utilizing the adaptive method of lines was formulated.Implementation of the model was discussed and its source code is available as an electronic supplement of the paper.Example results for single-and multi-component aerosol, with-and without adaptivity, for different vertical velocities, and different initial aerosol spectra were presented.
The key advantages of the introduced spectrum refinement method are i) the suppression of bogus sensitivity of the model results to the discretization parameters when working in low resolution; ii) the introduction of physically-based uncertainty-related parameters controlling the spectrum discretization independently of the resolution of model input; iii) the enhancement in resolution of the shape of the left edge of the model-predicted cloud drop spectrum; iv) the reduction of computational effort needed to obtain results of a given precision.The discussion of results obtained without adaptivity constitutes an assessment of the uncertainty of results obtained using classical MOL (see Fig. 5).
Presented method is potentially applicable to numerous fields of aerosol-related research where moving-sectional models of aerosol growth are used, including i) studies of closure among the CCN activation theory and in-situ measurements of cloud properties (e.g., Snider et al., 2003); ii) design and calibration of aerosol instruments (e.g., Nenes et al., 2001a); iii) research on the nucleation of aerosol in the atmosphere (e.g., Anttila et al., 2004); iv) estimations of the impact of aerosol properties and the kinetics of its growth on cloud albedo (e.g., Nenes et al., 2001b) and the related assessments Introduction

Conclusions References
Tables Figures

Back Close
Full of geoengineering techniques for global warming mitigation (e.g., Bower et al., 2006); v) investigations on cloud processing of aerosol (e.g., Romakkaniemi et al., 2006); vi) construction and evaluation of parameterizations to be employed in hydrodynamical models run at resolutions where CCN activation is a sub-scale process (e.g., Feingold and Heymsfield, 1992).An example of the last application is fitting parameters c T and k T of the Twomey's formula for CDNC to the model-predicted values (cf.Sect.3.4 and Fig. 8 herein).
Preliminary validation of the model completes a proof-of-concept stage of its development.Further validations are underway, and are aimed at comparing model results with measurements, and at enhancing model formulation.The implemented logic for efficiently handling variable number of bins shall be a valuable component for implementing an accurate treatment of the drop coagulation process using the movingsectional approach.These developments will be reported in future publications.[m] standard deviation of radius of cloud drops κ [c]  [1] chemical-composition parameter for spec.c κ [0] =1.28

Supplementary material related to this article is available online at
[1] CCN-derived value of κ for NaCl κ [1]   [1/kg] total specific conc. of aerosol in spec.c N [c,b]  [1/kg] specific conc. of aerosol in bin b of spec.c N b , N Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | d ) are invariant with time, are defined as a function of dry radius r d , separately for each chemical component c.Each dry spectrum has a corresponding wet spectrum n [c] w (r w ,t) defined analogously but varying with time.Both n Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | PDFs with piecewise constant functions.Here, each bin (numbered by b) is defined by a constant-in-time number of particles N ), respectively.The MOL is thus applied to calculate the evolution of r w (r d ,t) and T w (r d ,t) by discretizing r d and letting t to be the only independent variable.If the size spectra are discretized into adjacent bins then T b] wr = r [c,b+1] wl .Consequently, the state vector of the model assuming discretization of the size-spectrum into N b adjacent bins and N c chemical components consists ofDiscussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | r a = max(r min ,r [c,b] wl ), r b = min(r max ,r [c,b] wr ); r min was set at 1 µm and r max at 25 µm.Consequently: Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fortran-written VODE solver and a close relative of the Fortran-written LSODE solver.LSODE was used e.g. in Ghan et al. (1998), Nenes et al. (2001b) and Romakkaniemi are based on measurement-data fits summarised in O'Dowd et al. (1997, values for horizontal wind speed of 15 m/s were used).Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3. The shapes of both spectra (sea-salt and sulphate solution droplets) are plotted at six intermediate altitudes (histograms constructed with thick and thin blue lines, respectively).The cloud droplet specific number concentration (turquoise line) increases till about 140 s (35 m), and stays constant afterwards.The mean radius (thick orange line) increases continuously with height from the onset of drop formation till the end of the model run.The standard deviation of the drop radius (thin orange lines) stays approximately constant.The drops formed on sulphate aerosol dominate in the small-size part of the spectrum, while the sea-salt aerosol took part in formation of the largest drops.

3. 2
Sensitivity to the discretization resolution (without adaptivity) Apart from the physical parameters of the model (such as vertical velocity profile), and the initial values of model variables (such as initial temperature and the dry aerosol spectra), the model results depend on parameters controlling the discretization scheme 1290 Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | references therein).Comparison of Figs. 4 and 5 confirms that i) the sensitivity of the model-predicted values of CDNC to number of bins observed by Kreidenweis et al. (2003) and Korhonen et al. ( Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 27) where B denotes the Euler beta function, the 1.63×10 −3 coefficient corresponds to temperature of 10 • C and pressure of 800 hPa, the vertical velocity w is expressed in centimetres per second, c T is expressed in cm −3 and E is expressed in K.The upper bound of the drop concentration is expressed as a function of k T , c T and E max : CDNC < c T • (E max ) k T (Discussion Paper | Discussion Paper | Discussion Paper | The comparison clearly shows that the model predictions do fall below the upperbound solution within the whole range of vertical velocities.The distinction between the marine and continental aerosol regime is comparable.The shapes of all modelpredicted relations resemble the ones plotted using Twomey's formulae.The difference between the model-predicted and the limit maximum dew-point elevations is Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | =0.61 [1] CCN-derived value of κ for (NH 4 ) 2 SO 4 λ D [m] mean free path of water vapour in air λ K [m] length scale for heat transfer in air a kg/K] specific heat capacity of liquid water c p (q v ) [J/kg/K] specific heat capacity of moist air c pd =1005 [J/kg/K] specific heat capacity of dry air c pv =1850 [J/kg/K] specific heat capacity of water vapour c T [m −3 ] aerosol composition par. in Twomey's form.D(λ D ,r w ) [m 2 /s] diffusivity of water vapour in air (trans.r ég.) D 0 =2.21×10 5 [m 2 /s] diffusivity of water vapour in air (conti.r ég.) E , E max [K] dew-point elevation (maximum) F (x,t) n/a ODE system right-hand side f (x i ) n/a the function used in Newton-Raphson iteration g=9iteration enumerator K (λ K ,r w ) [J/m/s/K] thermal conductivity of air (trans.r égime) K 0 =2.4×10 −2 [J/m/s/K] thermal conductivity of air (conti.r égime) k [0,1,2, . . .] spectrum moment enumerator k T [1] aerosol composition par. in Twomey's form.LWC [kg/m 3 ] liquid water content l v (T ) [J/kg] latent heat of evapouration l v0 =2.5×10 6 [J/kg] latent heat of evapouration at the triple point Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | kg/m] dry aerosol specific conc.density n [c] w [1/kg/m] wet aerosol specific conc.density p [Pa] air-parcel pressure p 0 =611.73[Pa] water vapour pressure at the triple point p v [Pa] partial pressure of water vapour p v∞ (T ) [Pa] equil.vapour press.(plane surf., pure water) p vs (T,...at left/right bin boundary) w [m/s] air-parcel vertical velocity x n/a ODE set variables vector Z k [m k /kg] helper symbol in Eqs.(24) and (25)

Table 1 .
Aerosol distributions used in the example model run discussed in Sect.3.1.Data after O'Dowd et al. (1997), as used in Ghan et al. (1998).

Table 3 .
List of symbols.