(4.22)
The formulation can be readily extended to include a local source pi(r′, t′) and linear loss rate constant ℓi [s–1]. We then have
(4.23)
where pi and ℓi are applied over each trajectory arriving at point r at time t. Nonlinear chemistry cannot be accommodated in this framework because individual trajectories do not interact.
The Green function provides a general statement of source–receptor relationships and is used in a wide range of atmospheric chemistry applications. It is of specific use in receptor-oriented problems where we seek the contributions from a 2-D or 3-D source field to concentrations at a particular time and location. This forms the basis for simple inversion techniques relating observed atmospheric concentrations to surface fluxes (Enting, 2000). Approaches using a Green function are also used to characterize the transport history of air parcels. In the case of a single source at location r′ affecting the entire modeling domain, we can express the mixing ratio anywhere in the domain with a Green function in unit of inverse time:
(4.24)
The first moment of this Green function,
(4.25)
defines the mean age of air for transport from r′ to r. It has been used in particular to characterize transport times in the stratosphere for air originating at the tropical tropopause (Hall and Plumb, 1994).
4.2.9 Initial and Boundary Conditions
The Eulerian form of the continuity equation is a 4-D partial differential equation (PDE) in time and space, and solution requires specification of initial and spatial boundary conditions. The Lagrangian form expressed for a moving air parcel is a 1-D ordinary differential equation (ODE) in time, but spatial boundary conditions are still needed to describe the flow of air parcels.
Initial conditions describe the chemical concentrations over the 3-D domain at the beginning of the simulation and account therefore for the former evolution of these variables. They can be provided by a previous simulation using the same model, by a simulation from another model, by assimilated observations (Chapter 11), or by some mean climatological state. Initial conditions are often not well characterized and may not be consistent for the different species described in the model. The model simulation then carries that legacy of inconsistency during the initial phase of the simulation. Good practice is to initialize the model by conducting a simulation for some time period prior to the period of interest, sufficiently long that the inconsistency of initial conditions is dissipated. This is called spinning up the model. The length of the spin-up period depends on the characteristic timescales over which the species of interest respond. These timescales are discussed in Section 4.4.
Vertical boundary conditions are defined by imposed fluxes or concentrations at the boundaries of the domain. The lower boundary is usually the Earth’s surface, so that emission and dry deposition fluxes are appropriate boundary conditions. The dry deposition flux is usually computed on the basis of the concentration at the lowest model level with a specified deposition velocity, as given by (4.18). We show in Chapter 9 that this represents in fact a zero-concentration boundary condition at the surface, and can be made into a two-way exchange when the surface concentration is not zero. Further discussion of surface boundary conditions is presented in Chapter 9. A pure advection problem in which the model domain extends to the top of the atmosphere would require no other boundary condition since the continuity equation is first-order in space. Representation of vertical transport as turbulent diffusion (4.14) makes the continuity equation second-order in space and requires an additional boundary condition. Also and in practice, models operate over a limited vertical domain and there is some inflow through the top into the model domain. Therefore an upper flux boundary condition at the model top is needed. The model top may be sufficiently high that a zero-flux boundary condition is acceptable, or concentrations may be imposed at the model top based on information from observations or another model.
Lateral boundary conditions depend on whether the model is global or of limited domain. In a global model, the lateral boundary conditions are periodic as concentrations at longitude λ and latitude φ must be reproduced at (λ + 2π, φ + 2π). In a limited-domain model, lateral boundary conditions must be provided as fluxes or concentrations at the edges of the domain. Lateral boundary conditions are actually needed only for the upwind boundary of the domain, but since the wind varies in direction one needs in practice to prescribe boundary conditions at all edges of the domain. The superfluous boundary conditions at the downwind edges may lead to model noise if allowed to propagate into the model domain. It is standard to specify lateral boundary conditions as concentrations just outside the model domain, to be entrained into the model domain by the wind at the model boundaries. This ensures that only the upwind boundary conditions propagate into the model at any given time.
4.3 Continuity Equation for Aerosols
The continuity equation can be applied to different size classes and chemical components of aerosols in the same way as for gases. The transport terms are the same, since the particles are sufficiently small that they are advected by the wind in the same way as gases, except that one should add a sink term from gravitational settling in the case of very large particles (>10 μm in the troposphere, > 1 μm in the stratosphere). The chemical term in the continuity equation is of the same form as for gases if we only simulate the total mass concentrations of aerosol chemical components (without regard to their size distributions).
If we are interested in describing the evolution of the aerosol size distribution, however, we need to introduce additional terms in the continuity equation to account for aerosol nucleation, condensational growth, coagulation, and cloud interactions. This ensemble of processes is called aerosol microphysics. For accounting purposes it is more convenient to use volume rather than radius as the independent variable to characterize the aerosol size distribution. We thus define a volume distribution function nN(V) such that nN(V) dV represents the number concentration of particles in the volume range [V, V + dV]. The relationship between nN(V) and nV(r) defined in Chapter 3 can be derived by equaling the number of particles in the volume range [V, V + dV] with that in the equivalent size range [r, r + dr]:
(4.26)
The local evolution of the aerosol size distribution can be written as
(4.27)
where nucleation describes the formation of new particles from the gas phase, condensation/evaporation describes the exchange of mass between the gas phase and the particles, and coagulation describes the formation of larger particles from the collision of smaller particles by Brownian diffusion. Representation of aerosol–cloud interactions may require additional terms.
The nucleation term is dependent on the supersaturation of species in the gas phase and on the stability of successively larger molecular clusters formed from these gas-phase molecules. Nucleation acts as a source of new particles at the low end of the size distribution (~10–3 μm). The condensation/evaporation term can be calculated from knowledge of the condensation growth rate I(V) = dV/dt of particles of volume V. Consider a volume element [V, V + dV]. The particles growing into that volume element over time dt are those that were in the volume element [V, V – I(V)dt]; their number concentration is nN(V) I(V) dt. Similarly, the number concentration of particles growing out of the volume element is n(V + dV) I(V + dV) dt. Thus the change in the volume distribution function is
(4.28)
The coagulation term is defined by the frequency of collisions between particles (collision usually results in coagulation). We characterize the collision frequency by a coagulation coefficient β(V, V′) [m3 particle–1s–1] representing the rate constant at which particles of volume V collide with particles of volume V′. In this manner, the rate of production of particles of volume (V + V′) by collision of particles of volume V and V′ is given by β(V, V′) nN(V) nN(V′) (dV)2. To determine the change with time in the number concentration of particles in the volume element [V, V + dV]
due to coagulation processes we consider collisions with particles over the entire range of the size distribution and account for both production and loss of particles out of the volume element:
(4.29)
where the ½ coefficient on the first term of the right-hand side is to avoid double-counting. Physical formulations for I(V), β(V, V′), and aerosol nucleation are presented by Seinfeld and Pandis (2006).
Standard numerical algorithms for aerosol microphysics approximate the size distribution by a series of square functions called sections or “bins.” These are called sectional models. Accurate representation of aerosol microphysics generally requires more than 30 bins for each chemical component of the aerosol, so including multiple components can become computationally cumbersome. An alternative is to make some assumption about the shape of the size distribution and track the evolution of the low-order moments of that distribution. See Section 5.6 for additional information on the methods used in aerosol models.
4.4 Atmospheric Lifetime and Characteristic Timescales
4.4.1 Atmospheric Lifetime
The atmospheric lifetime of a chemical species is defined as the average time that the species remains in the atmosphere before it is removed by one of its sinks. The term residence time is equivalently used, most often when the sink involves deposition to the surface. If a species i has a local mass concentration ρi [kg m–3] and loss rate Li [kg m–3 s–1], the local atmospheric lifetime is given by
(4.30)
We are often interested in the mean atmospheric lifetime averaged over some atmospheric domain V and time period T. This is obtained by summation of the concentrations and loss rates:
(4.31)
For example, the global mean atmospheric lifetime is defined by summing over the whole atmosphere and a full annual cycle.
The atmospheric lifetime of a species is a very useful thing to know because it characterizes the spatial and temporal variability of its concentration. A species with a short lifetime will have strong concentration gradients defined by fluctuations of its sources and sinks, while a species with a long lifetime will be well mixed. A species with a short lifetime will respond rapidly to changes in sources or sinks, while a species with a long lifetime will respond more slowly.
The loss rate of a species is usually proportional to its concentration (first-order loss). We then write Li = li ρi, where li [s–1] is a loss rate coefficient, and derive τi = 1/li. If there is no compensating production, the evolution of the concentration is expressed by
(4.32)
with solution
ρi(t) = ρi(0) exp [−li t] = ρi(0) exp [−t/τi]
(4.33)
Thus we see that τi defines a characteristic e-folding timescale for decay: ρi(τ) = ρi(0)/e = 0.37 ρi(0). τi is sometimes called the e-folding lifetime to distinguish it from the half-life t1/2 used in the radiochemistry literature to describe 50% loss of an initial amount of radioactive material. For a first-order loss, τ and t1/2 are related by τ = t1/2/ln2 = 1.44 t1/2.
Different processes can contribute to the removal of a species from the atmosphere, including chemical loss, transport, wet scavenging, and surface uptake, as described by the different terms in the continuity equation (Section 4.2). If several processes contribute to the sink and all are first-order, the overall loss coefficient li is the sum of the loss coefficients for the individual processes. Thus we have
(4.34)
where li,q is the loss coefficient associated with process q. One can similarly define a lifetime τi,q = 1/li,q associated with process q. The overall lifetime τ resulting from all loss processes is
(4.35)
Thus the loss coefficients for individual processes add in series while the corresponding atmospheric lifetimes add in parallel.
We now examine the formulation of lifetimes for individual processes. Consider first the chemical loss of species i by reaction with species j with a rate constant k. The corresponding chemical lifetime is
(4.36)
For a more general case where species i reacts with several species j (rate constants kj) and also photolyzes with a photolysis frequency J, the chemical lifetime is given by
(4.37)
Chemical lifetimes in the atmosphere range from less than a second for highly reactive radicals to practically infinite for noble gases.
Consider now the timescales for transport of species i. For one-dimensional advection along direction x with a velocity u, the continuity equation is
(4.38)
Let us view advection as representing a sink for species i, such as ventilation from a source region. In that case we can define a loss coefficient associated with advection as per (4.30): ladv = ∂(u ρi)/∂x / ρi. We approximate the divergence term as ∂(u ρi)/∂x ≈ u ρi/D, where D represents a characteristic distance over which the transport flux varies. The lifetime associated with advection is then given by
(4.39)
For a source region of dimension D, (4.39) defines the lifetime against ventilation out of the source region by the mean wind. In a more general case, advection acts as both a source and a sink for the concentration at a given location. τadv in (4.39) can be viewed more generally and usefully as a transport timescale for species i to be transported over a distance D.
For small-scale turbulent mixing parameterized with a turbulent diffusion coefficient K, the continuity equation is
(4.40)
Again one can make an order-of-magnitude approximation of the diffusion term as Kρi/D2, and the corresponding mixing timescale is
(4.41)
Transport timescales can be derived from knowledge of the atmospheric circulation and also from observations of chemical tracers. Typical values for the troposphere are given in Figure 4.4. Horizontal transport is fastest in the longitudinal direction because of strong winds; thus air can circumnavigate the world in a given latitudinal band on a timescale of a month. Meridional mixing within a hemisphere takes longer, on a timescale of about three months, and exchange across hemispheres requires the order of a year because of the lack of a strong thermal forcing gradient across the Equator. Vertical transport in the troposphere is driven largely by buoyant mixing with characteristic timescales of a day for mixing of the planetary boundary layer (PBL, surface to 1–3 km) and a month for the full troposphere. Tropospheric air has a residence time of 5–10 years against transport to the stratosphere, while stratospheric air has on average a residence time of 1–2 years against transport to the troposphere.
Figure 4.4 Typical time constants associated with global atmospheric transport.
From Jacob (1999).
The lifetime of species i against dry deposition (surface uptake) for a well-mixed atmospheric layer extending from the surface to altitude h is expressed as a function of the deposition velocity wD,i by
(4.42)
For a gas removed efficiently at the surface, wD,i is of the order of 1 cm s–1 so that the lifetime against dry deposition in the PBL is of the order of a day.
The lifetime of water-soluble gases and aerosols against wet scavenging can be roughly estimated from the frequency of precipitation. It is typically of the order of a week in the lower troposphere, longer in the upper troposphere. The estimate is complicated by the coupling between atmospheric transport and precipitation and by differences in scavenging efficiencies for different forms of precipitation.
4.4.2 Relaxation Timescales in Response to a Perturbation
An important application of the concept of atmospheric lifetime is to determine the time required for the atmosphere to relax in response to a perturbation. Consider a species i initially at steady state to which an instantaneous perturbation Δρi(0) is applied at time t = 0. If the loss is linear with coefficient li, then the perturbation decays as Δρi(t) = Δρi(0)exp[–lit] and the corresponding timescale is defined by the species lifetime τi = 1/li. The situation is more complicated if the loss is nonlinear due to chemical coupling between species. To analyze the behavior of a nonlinear c
oupled system, let us consider the general kinetics equation
(4.43)
in which ρ (ρ1, ρ2, … ρN)T represents the concentration vector for N chemical species and A(ρ) is an N × N matrix operator containing the kinetic expressions for the chemical production and destruction rates of each species. If a small perturbation Δρ is applied to the system, the response of the perturbed state is determined by the linearized equation (Prather, 2007)
(4.44)
where the (N × N) Jacobian matrix J of operator A is the first term in the Taylor expansion of A(ρ + Δρ). We neglect the higher-order terms O(Δρ2) in what follows. If the perturbation is an eigenvector of J with eigenvalue λ (Box 4.1), then
(4.45)
and the solution is
Δρ(t) = Δρ(0) exp (λt)
(4.46)
If the N species are linearly independent, then J has full rank with N linearly independent eigenvectors aj each with an eigenvalue λj (j = 1, N). In the general case, an initial perturbation Δρ(0) can be decomposed on the basis of eigenvectors with coefficients αj:
(4.47)
so that the time evolution of the perturbation is given by
(4.48)
Box 4.1 Eigenvalues and Eigenvectors
The vector a = (a1, …aN)T is an eigenvector of the square matrix A (N × N) if it satisfies the equation
Aa = λa
Modeling of Atmospheric Chemistry Page 15