(5.129)
where I(V) = dV/dt is the condensation growth rate. Equation (5.129) is called the condensation equation. It is mathematically similar to the advection equation and numerical algorithms face the same difficulties, discussed in Chapter 7. Numerical diffusion associated with the algorithm can lead to erroneous damping of the peak values of the distribution. Numerical dispersion leads to wakes around peak values, particularly near fronts in the size distribution (Seigneur et al., 1986).
The condensation growth rate for species i is proportional to the difference between the bulk vapor pressure pi (far from the particle) and the equilibrium vapor pressure peq,i. It is expressed by (Seinfeld and Pandis, 2006):
(5.130)
where Dg(V) [cm2 s–1] is the Brownian diffusion coefficient for the particle, and mi [kg mol–1] is the molecular mass of species i. Function f(Kn, α) is a correction factor for the non-continuum regime (Section 5.5.2) that depends on the Knudsen number Kn and the mass accommodation coefficient α.
Coagulation
Coagulation is the process by which two particles that collide by Brownian motion stick together to form a new, larger particle. It shifts the size distribution toward larger sizes and reduces the number of smaller particles.
The coagulation rate Ji,j [m–3 s–1] resulting from the collisions between two particles i and j is proportional to the number density Ni and Nj [m–3] of these two particles:
Ji , j = βi , jNiNj
(5.131)
where βi,j [m3 s–1] is the binary coagulation coefficient. In the continuum regime with Kn ≪ 1 (particle diameter considerably larger than the mean free path), βi,j is given by
βi , j = 2π(Di + Dj)(Dp , i + Dp , j)
(5.132)
where Dp,i is the diameter of particle i and Di is the Brownian diffusion coefficient given by the Stokes–Einstein formula
(5.133)
where μ is the viscosity of air. Thus, in this case,
(5.134)
In the free molecular regime with Kn ≫ 1 (very small particles with diameters considerably smaller than the mean free path), the coagulation coefficient becomes
(5.135)
where vi = (8kT/πmp,i)1/2 is the mean thermal velocity with mp,i the mass of particle i. In the transition regime, one generally adopts the continuum regime formula with a correction factor (Seinfeld and Pandis, 2006). Coagulation is most effective when the sizes of the two particles that collide are very different (i.e., collision between a small particle with high thermal velocity and a large particle that provides a large area for collision).
The rate of change in the aerosol size distribution resulting from coagulation is given by
(5.136)
where V0 denotes the volume of the smallest particles considered in the size distribution (associated typically with nucleation, cf. (5.128)). Factor ½ is introduced to avoid double counting.
5.6.2 Representation of the Size Distribution
Computation of the microphysical terms in the continuity equation for aerosols requires that the size distribution be approximated with a limited number of parameters. This is illustrated in Figure 5.17 with discrete, spline, sectional, modal, and monodisperse approximations. The discrete representation provides a fine resolution of the size distribution function with a value for each discrete size. The spline defines a single continuous size distribution function over the whole size range. The sectional representation partitions the size distribution into discrete size intervals called “bins,” with fixed values of the distribution functions within these intervals. The modal representation superimposes several continuous functions covering different ranges of the size distribution, one for each mode. The monodisperse representation assigns a single size to each mode. We elaborate here on the most popular methods.
Figure 5.17 Numerical approximations of the aerosol size distribution function.
From Whitby et al. (1991).
Sectional method
In the sectional method (Gelbard et al., 1980; Adams and Seinfeld, 2002), the aerosol size domain is divided into a discrete number of bins K within which the size distribution functions are assumed to be constant. Thus the size distribution function is represented by K parameters Nk(k = 1, K) representing the number concentration of particles in bin k bounded by volumes [Vk, Vk+1]:
(5.137)
Applying the continuity equation (5.127) to these parameters yields a system of K coupled ODEs:
(5.138)
Nucleation (J) provides a source of particles to the smallest bin (k = 1). Condensation growth (p) and evaporation (γ) provide source/sink terms for adjacent bins. Coagulation couples the whole size distribution. This system of coupled ODEs is mathematically equivalent to the system representing chemical production/loss terms for different species in a chemical mechanism. Numerical algorithms for solving such systems are presented in Chapter 6.
Moments and modal
The detailed aerosol size distribution generally contains more information than required to address a specific problem. Rather than calculating its time evolution by solving the conventional continuity equation (5.127), it is often sufficient and computationally more efficient to estimate the low-order moments of the size distribution (Friedlander, 1977). The moment of order k is defined by
(5.139)
where n(Dp) dDp represents the number of particles (assumed to be spherical) per unit volume of air in the diameter size range [Dp, Dp + dDp]. The evolution equations for moment Mk of the aerosol distribution (called moment dynamics equations or MDEs) are derived by starting from the aerosol continuity equation (5.127) with n(Dp) as the independent variable, multiplying each term by , and integrating each term over all particle diameters (Whitby and McMurry, 1997). The equations can be solved if all terms can be expressed with only moments as the dependent variables. This requires either an assumption on the mathematical form of the size distribution (see modal method below) or approximations to the terms (closure relations) to force them to be expressed in terms of moments. The latter approach describes the method of moments (MOM) and makes no a-priori assumptions on the form of the size distribution. However, there may be large errors associated with the closure relations.
In the modal method (Whitby, 1978), the aerosol size distribution is specified as the superimposition of a limited number K of functional forms, each representing a particular aerosol mode. For example, a given mode may be defined by a log-normal size distribution with imposed median diameter Dm,k and geometric standard deviation σg,k (5.126). The contribution of each mode to the overall size distribution is then defined by the number concentration of particles in that mode. Figure 5.18 gives an example. The microphysical terms of the continuity equation are integrated over each mode and transfer particles between modes. The evolution of the aerosol size distribution can be calculated using low-order moments by noting that these moments characterize the parameters of the distribution. For example, a log-normal distribution is fully defined by its first three moments (k = 0, 1, 2 in (5.139)) characterizing the total number of particles, the mean diameter, and the standard deviation of the distribution. The modal method is computationally much faster than the sectional method but relies on the suitability of decomposing the actual aerosol size distribution along specified modes.
Figure 5.18 Modal distributions of different aerosol components (sulfate, black carbon, organic matter, sea salt, dust). The aerosol size distribution for each component is represented by the number concentrations in seven log-normal modes. The smallest (nucleation) mode is exclusively sulfate. r is particle radius. Based on the model of Stier et al. (2005).
From Stier and Feichter, personal communication.
References
Adams P. J. and Seinfeld J. H. (2002) Predicting global aerosol size distributions in general circulation models, J. Geophys. Res., 107, 4370, doi:10.1029/2001JD001010.
Brasseur G. P. and Solomon S. (2005) Aeronomy of the Middle Atmosphere: Chemistry and Physics of the Stratosphere
and Mesosphere, 3rd edition, Springer, New York.
Chabrillat S. and Kockarts G. (1997) Simple parameterization of the absorption of the solar Lyman-α line, Geophys. Res. Lett., 24 (21), 2659–2662, doi: 10.1029/97GL52690, correction: Geophys. Res. Lett. 25 (1), 79, doi: 10.1029/97GL03569.
Chandrasekhar S. (1950) Radiative Transfer, Oxford University Press, Oxford (Reprinted by Dover Publications, 1960).
Donahue N. M., Robinson A. L., Stanier C. O., and Pandis S. N. (2006) Coupled partitioning, dilution, and chemical aging of semivolatile organics, Environ. Sci. Technol., 40, 2635–2643.
Elsasser W. M. (1938) Mean absorption and equivalent absorption coefficient of a band spectrum, Phys. Rev., 54, 126–129.
Fang T. M., Wofsy S. C., and Dalgarno A. (1974) Capacity distribution functions and absorption in Schumann–Runge bands of molecular oxygen, Planet Space Sci., 22, 413–425.
Friedlander S. K. (1977) Smoke, Dust and Haze: Fundamentals of Aerosol Behavior, Wiley, New York.
Fu Q. and Liou K. N. (1992) On the correlated k-distribution method for radiative transfer in nonhomogeneous atmosphere, J. Atmos. Sci., 49, 2139–2156.
Gelbard F., Tambour Y., and Seinfeld J. J. (1980) Sectional representation for simulating aerosol dynamics, J. Colloid Interface Sci., 76, 357–375.
Gijs A., Koppers A., and Murtagh D. P. (1997) Model studies of the influence of O2 photodissociation parameterizations in the Schumann–Runge bands on ozone related photolysis in the upper atmosphere, Ann. Geophys., 14, 68–79.
Goody R. M. (1952) A statistical model for water vapor absorption, Quart. J. Roy. Met. Soc., 78, 165–169.
Goody R. (1995) Principles of Atmospheric Physics and Chemistry, Oxford University Press, Oxford.
Goody R., West R., Chen L., and Crisp D. (1989) The correlated-k method for radiation calculations in nonhomogeneous atmosphere, J. Quant. Spectrosc. Radiat. Transfer, 42, 539–550.
Heintzenberg J., Raes F., and Schwartz S. E. (2003) Tropospheric aerosols. In Atmospheric Chemistry in a Changing World (Brasseur G. P., Prinn R. G., and Pszenny A. P., eds), Springer, New York.
Jacob D. J. (1986) The chemistry of OH in remote clouds and its role in the production of formic acid and peroxymonosulfate, J. Geophys. Res., 91, 9807–9826.
Jacob D. J. (2000) Heterogeneous chemistry and tropospheric ozone. Atmos. Environ., 34, 2131–2159.
Joseph J. H., Wiscombe W. J., and Weinman J. A. (1976) The delta-Eddington approximation for radiative flux transfer, J. Atmos. Sci., 33, 2452–2459.
Kockarts G. (1994) Penetration of solar radiation in the Schumann–Range bands of molecular oxygen: A robust approximation, Ann. Geophys., 12 (12), 1207–1217, doi: 10.1007/BF03191317.
Lenoble J. (1977) Standard Procedures to Compute Atmospheric Radiative Transfer in a Scattering Atmosphere, Vol. I, International Association of Meteorology and Atmospheric Physics (IAMAP), Boulder, CO.
Liou K. N. (2002) An Introduction to Atmospheric Radiation, Vol. 84, 2nd edition Academic Press, New York.
López-Puertas M. and Taylor F. W. (2001) Non-LTE Radiative Transfer in the Atmosphere, World Scientific Publishing, Singapore.
Lorenz L. V. (1890) Lysbevaegelsen i og uden for en af plane Lysbolger belyst Kugle, Det Kongelige Danske Videnskabernes Selskabs Skrifter, 1, 1–62.
Martin R.V., Jacob D.J., Yantosca R.M., Chin M., and Ginoux P. (2003) Global and regional decreases in tropospheric oxidants from photochemical effects of aerosols, J. Geophys. Res., 108, 4097.
Martin S. T. (2000) Phase transitions of aqueous atmospheric particles, Chem. Rev., 100, 3403–3453.
Mie G. (1908) Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen, Ann. Phys., 330, 377–445.
Minschwaner K. and Siskind D. E. (1993) A new calculation of nitric oxide photolysis in the stratosphere, mesosphere, and lower atmosphere, J. Geophys. Res., 98 (111), 20401–20412, doi: 10.1029/93JD02007.
Minschwaner K., Salawitch R. J., and McElory M. B. (1993) Absorption of solar radiation by O2: Implications for O3 and lifetimes of N2O, CFCl3, and CF2Cl2, J. Geophys. Res., 98, 10543–10561, doi: 10.1029/93JD00223.
Petty G. W. (2006) A First Course in Atmospheric Radiation, 2nd edition, Sundog Publications, Madison, WI.
Rayleigh L. (1871) On the light from the sky, its polarization and colour, Phil. Mag., 41, 107–120.
Santillana M., Le Sager P., Jacob D. J., and Brenner M. P. (2010) An adaptive reduction algorithm for efficient chemical calculations in global atmospheric chemistry models, Atmos. Environ., 44, 4426–4431.
Schwartz S. E. (1986) Mass transport considerations pertinent to aqueous-phase reactions of gases in liquid-water clouds. In Chemistry of Multiphase Atmospheric Systems (Jaeschke W., ed.), Springer-Verlag, Berlin.
Schwartz S. E. and Freiberg J. E. (1981) Oxidation of SO2 in aqueous droplets: Mass-transport limitation in laboratory studies and the ambient atmosphere, Atmos. Environ., 15, 1129–1144.
Seigneur C., Hudischewskj A. B., Seinfeld J. H., et al. (1986) Simulations of aerosol dynamics: A comparative review of mathematical models, Aerosol Sci. Technol., 5 (2), 205–222.
Seinfeld J. H. and Pandis S. N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, Wiley, New York.
Shaw J. (1953) Solar radiation, Ohio J. Sci., 53, 258.
Smith F. L. III and Smith C. (1972) Numerical evaluation of Chapman’s grazing incidence integral Ch (X, χ), J. Geophys. Res., 77, 19, 3592–3597, doi: 10.1029/JA077i019p03592.
Stamnes K., Tsay S. C., Wiscombe W., and Jayawerra K. (1988) Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl. Opt., 27, 2502–2509.
Stier P. Feichter J., Kinne S., et al. (2005) The aerosol–climate model ECHAM5-HAM, Atmos. Chem. Phys., 5, 1125–1156.
Whitby K. T. (1978) The physical characteristics of sulfur aerosols, Atmos. Environ., 12, 135–159.
Whitby E. R. and McMurry P. H. (1997) Modal aerosol dynamics modeling, Aerosol Sci.Technol., 27, 673–688.
Whitby E. R., McMurry P. H., Shankar U., and Binkowski F. S. (1991) Modal Aerosol Dynamics Modeling, Atmospheric Research and Exposure Assessment Laboratory, Research Triangle Park, NC.
Woods T. N., Prinz D. K., Rottman G. J., et al. (1996) Validation of the UARS solar ultraviolet irradiances: Comparison with the ATLAS 1 and 2 measurements, J. Geophys. Res., 101 (D6), 9541–9569, doi: 10.1029/96JD00225.
Zhu X., Yee J.-H., Lloyd S. A. and Storbel D. F. (1999) Numerical modelling of chemical–dynamical coupling in the upper stratosphere and mesosphere, J. Geophys. Res., 104, 23995–24011, doi: 10.1029/1999JD900476.
6
Numerical Methods for Chemical Systems
6.1 Introduction
Solving the 3-D continuity equations for chemical species in atmospheric models requires splitting of the transport and chemistry operators. We present in this chapter an overview of numerical methods for the chemistry operator, which solves the evolution of the system driven by chemical kinetics independently of transport. Complexity arises from the large number of coupled species in standard mechanisms for atmospheric chemistry, with time constants ranging over many orders of magnitude. The associated computational requirements are very high and this is a major challenge for the inclusion of atmospheric chemistry in Earth system models.
For a chemical mechanism involving K chemically interacting species, the task of the chemistry operator is to solve the following initial value problem over time step Δt,
(6.1)
where ψ = (ψ1, ψ2, …, ψK)T is the vector of concentrations for the K species with known initial value ψ(to), and s is a vector of chemical production and loss rates. Each component of s is a sum of terms describing the rates of individual reactions. Equation (6.1) describes a system of coupled ordinary differential equations (ODEs) with time as the only coordinate. There is no spatial dependence since the chemical evolution is a function of local concentrations only. Although we refer to solution of (
6.1) as the “chemistry operator,” s may also include non-chemical terms such as emission, precipitation scavenging, and dry deposition rates. Any local process independent of transport (and hence with no spatial dependence) can be included in the chemistry operator. When solving for aerosol microphysics, ψ may represent particle concentrations of different sizes with s including particle formation and growth terms (Section 4.3).
Nonlinearity arises in (6.1) because the rates of bimolecular and three-body reactions involve products of concentrations. This nonlinearity can be highlighted by rewriting (6.1) as:
(6.2)
Here A is a (K × K) diagonal matrix of unimolecular reaction rate coefficients in the mechanism, Qi is a (K × K) upper triangular matrix of bimolecular and three-body rate coefficients for reactions producing or consuming species i, ei is the ith column of the identity matrix of dimension K (zeros except for 1 in row i), and f = (f1, f2, …, fK)T is an independent forcing term. Nonlinearity is introduced by ΨTQiΨ, which is a summation of terms of form qijkΨjΨk producing or consuming species i. The form qijkΨjΨk applies to both bimolecular and three-body reaction rates since the concentration of the “third body” in a three-body reaction is not computed from the mechanism but is instead independently specified (see Chapter 5).
The general approach for obtaining ψ(to + Δt) from knowledge of ψ(to) is to use a finite difference approximation of the temporal derivative in (6.1). Differences lie in the way that s(ψ, t) is estimated. Fast ODE solvers use explicit methods where s is calculated on the basis of the known concentrations at to and previous time steps. In these solvers, as we will see, the time step must be smaller than the shortest time constant in the system in order to maintain stability. This is a major obstacle for atmospheric chemistry applications because radical species central to the chemical mechanisms have very short lifetimes. The systems of ODEs describing atmospheric chemistry mechanisms are stiff (Box 6.1), that is, the time constants for the different species coupled through the mechanism vary over many orders of magnitude. The stiffness is defined by the stiffness ratio R = τL/τS where τL and τS are the longest and shortest time constants in the system, corresponding roughly to the longest and shortest lifetimes of species in the mechanism (see Section 4.4). Solution with an explicit solver requires time steps Δt ~ τS, but we are generally interested in solutions over timescales ~τL. Thus the number of time steps required with an explicit solver is of order R. A typical mechanism for atmospheric chemistry may have stiffness R ~ 108, making for a formidable computational problem on a 3-D model grid.
Modeling of Atmospheric Chemistry Page 32