for some scalar λ called an eigenvalue of A. A non-degenerate matrix (determinant |A| ≠ 0) has N linearly independent eigenvectors forming a base and N eigenvalues (one for each eigenvector).
The above equation can be rearranged as
(A – λI)a = 0
where I is the identity matrix. The equation must have an infinite number of solutions since any scalar multiplier of a is a solution. It follows that the matrix A – λI must be degenerate, so that its determinant must be zero:
|A – λI| = 0
When expanded, this equation takes the form of a polynomial equation of Nth degree in λ. The N roots of this equation are the eigenvalues of the system. Replacing in A a = λ a then yields the corresponding eigenvector for each eigenvalue.
As an illustrative example, consider the matrix
The eigenvalues are the roots of or
λ2 – 4 λ + 3 = 0
The roots of this quadratic equation are λ = 1 and λ = 3.
The eigenvectors are obtained by solving the linear equation
for the two values of λ. We find for λ = 1 and for λ = 3.
We see that the eigenvalues define the characteristic timescales for the responses to the perturbation, and the eigenvectors define the modes over which the timescales apply. If all eigenvalues have negative real components, then Δρ will relax back to steady state over a suite of timescales [–1/λ1, …‒1/λN] that define the characteristic timescales of the system. This can be used in particular to identify the longest timescale for response. If any of the eigenvalues has a positive real component then an initial perturbation to that mode will grow exponentially with time and the system is unstable (it is explosive). If any of the eigenvalues has a non-zero imaginary component then a perturbation to that mode will induce an oscillation that may grow or decay depending on the sign of the real component. Eigenvalues for realistic mechanisms used in atmospheric chemistry models are all real and negative, so that the mechanisms are stable against perturbations. This should not be surprising as chemical systems generally follow Le Chatelier’s principle (“any perturbation to equilibrium prompts an opposing response to restore equilibrium”). Oscillatory or unstable chemical behavior may occasionally occur under unusual conditions but these tend to be rapidly dissipated by model mixing.
4.5 Conservation Equations for Atmospheric Dynamics
Solution of the chemical continuity equation (4.1) requires information on winds and turbulence to compute the transport terms. This information must be provided by a meteorological model. We discussed in Chapter 1 how chemical transport models can operate either “online,” integrated within the meteorological model, or “offline,” using archived output from the meteorological model. Here we describe the conservation equations for atmospheric dynamics that form the basis of meteorological models.
The dynamical properties of the atmosphere are determined by the fundamental principles of mass, momentum, and energy conservation. Mass cannot be created nor destroyed, momentum can be changed only through the application of a force, and internal energy can be altered only through the existence of a heat source or sink, or by performance of work. The properties are described at any spatial location and time by six dependent variables: the pressure p [Pa], the density ρa [kg m–3], the absolute temperature T [K], and the wind vector v [m s–1] with its three components (u, v, and w).
Pressure, density, and temperature are related by the ideal gas law:
p = ρaRT
(4.49)
where R is the gas constant (287 J kg–1 K–1). Solving the dynamical system requires five additional equations. These describe the evolution of mass, momentum (in the three spatial directions), and energy for the compressible fluid.
Since the Earth is rotating around its north–south axis, the most appropriate reference frame to describe air motions for an observer located on the Earth is not an inertial frame attached to the center of the planet, but rather a rotating coordinate system attached to the surface of the Earth. The geometric coordinates used in this case are denoted (x, y, z). Variables x and y represent geometric distances along parallels and meridians, respectively. Variable z is the geometric altitude. When spherical coordinates are used, the equations are expressed as a function of longitude λ and latitude φ rather than longitudinal and meridional distances x and y. We have
(4.50)
(4.51)
(4.52)
where r represents the geometric distance from the center of the Earth. The three wind components, and specifically the curvilinear velocities along a latitude circle (u) and along a meridian (v), and the linear velocity along the vertical (w) are expressed as
(4.53)
(4.54)
(4.55)
The total derivative of a quantity Ψ is expressed as
(4.56)
In local geometric coordinates (x, y, z), the total derivative is written as
(4.57)
In spherical coordinates (λ, φ, z), an equivalent formulation of the total derivative is
(4.58)
In many applications, the thickness of the atmosphere z is assumed to be small compared to the Earth’s radius a. In this shallow atmosphere approximation, one assumes z ≪ a and r = a. For models extending to the upper atmosphere, this approximation may not be valid.
4.5.1 Mass
Mass conservation for air is expressed by the continuity equation. In flux form, it is written as
(4.59)
This equation can be expanded as
(4.60)
with the sum of the first two terms being equal to the total derivative of ρa. Thus we write
(4.61)
which shows that the change in the density of a fluid parcel is proportional to the velocity divergence. Under the shallow atmosphere approximation, (4.59) and (4.61) are expressed in spherical coordinates as
(4.62)
and
(4.63)
If the fluid is incompressible (dρa/dt = 0), the continuity equation (4.61) becomes simply ∇ ⋅ v = 0. The velocity field is said to be non-divergent and the fluid density is conserved along the flow. Air is compressible and consequently this condition is not verified in the atmosphere.
4.5.2 Momentum
Newton’s second law applied to a continuum medium (i.e., with infinitely divisible fluid parcels) states that the acceleration of a small fluid element results from (1) body forces f that affect the whole fluid element (gravity, and electromagnetic forces if the fluid is ionized) and (2) stress forces that act on the surface of the fluid element and represent interactions with the rest of the fluid. These stress forces are expressed as the gradient of a second-order stress tensor σ. The balance of momentum in an inertial frame of reference (attached to the center of the Earth) is therefore provided by the differential equation (McWilliams, 2006; Neufeld and Hernandez-Garcia, 2010):
(4.64)
where ρa is the mass density of the fluid and V is the velocity in the inertial frame of reference. The only body force considered here is gravity (f = ρa ga). Further, one distinguishes between an isotropic term (diagonal elements of the stress tensor representing the pressure of the fluid: σii = –p) and an anisotropic component expressed by a tensor τ that accounts for interactions (viscosity) between the different fluid layers that move relative to each other. The resulting equation, known as the Cauchy momentum equation, is
(4.65)
To solve this vector equation, an expression must relate the stress terms to the velocity of the flow. In the hypothesis of a Newtonian flow, the stress τ is assumed to be proportional to the gradient of the velocity perpendicular to the shear. The proportionality coefficient μ [Pa s or kg m–1 s–1] is known as the dynamic viscosity and is a scalar when the fluid is isotropic (otherwise, it is a tensor). If the fluid is incompressible with a constant scalar viscosity μ, the elements τi,j of the stress tensor are expressed as
(4.66)
if xj is the jth spatial coordinate and vi is the fluid’s velo
city in the direction of axis i. The stress tensor is therefore related to the velocity vector by (Neufeld and Hernandez-Garcia, 2010):
τ = μ(∇ ⋅ V + ∇ ⋅ VT)
(4.67)
The assumption of incompressibility (∇ ⋅ V = 0) is appropriate for compressible fluids such as air if the velocity of the flow corresponds to Mach numbers smaller than about 0.3. Under this assumption, one obtains the Navier–Stokes equation written here in an inertial frame of reference:
(4.68)
In the more general case of the compressible flows encountered in the atmosphere, the friction term takes a more complicated form, and the Navier-Stokes equation becomes
with the viscosity μ assumed to have a constant value.
Here, we have expressed the material (total) derivative as the sum of the local acceleration ∂V/∂t and the convective acceleration (V ⋅ ∇)V. This equation provides the three components of the velocity field V(r, t) at a given point in space r and time t, and is used in fluid dynamics to address many different questions at various scales. The convective acceleration produced by spatial gradients in the velocity introduces a nonlinear component in the equation that can be the source of chaotic behavior known as turbulence (Box 4.2). The viscosity term operates as a linear diffusion term for momentum. The solution of this equation requires that appropriate conditions be prescribed at the boundary of the domain. The traditional condition adopted in many fluid dynamics problems is zero flow across solid boundaries and zero fluid velocity at the boundaries (“no-slip” condition).
Box 4.2 The Navier–Stokes Equation, the Reynolds Number, and Turbulent Flows
The nature of fluid flow can be assessed by comparing the importance of the inertial (or convective) term (v ⋅ ∇)v with the viscosity (diffusive) term ν∇2v in the Navier–Stokes equation (4.68). Here ν = μ/ρa represents the kinematic viscosity [m2 s–1]. To address this question, we define reduced quantities
where L and U are characteristic length and velocity scales, and apply this transformation to the Navier–Stokes equation. Omitting the prime signs for clarity, we obtain the non-dimensional equation
where
is the dimensionless Reynolds number, a measure of the ratio between the inertial force U(U/L) and viscous forces (νU/L2). Parameter ν = μ/ρ [m2 s–1] denotes the kinematic viscosity with a value of 1.5 × 10–5 m2 s–1 at 20 °C. Reynolds (1883) showed that as the value of Re increases, the motions become progressively more complex with a gradual transition from laminar to turbulent conditions. For low Reynolds numbers (less than 10), the steady-state flow results from a balance between pressure gradient and viscous forces. For larger Reynolds numbers, the nonlinear inertial term becomes dominant and chaotic eddies, vortices, and other instabilities are produced. Under this turbulent regime, viscosity is too small to dissipate the large-scale motions and the kinetic energy of the flow “cascades” progressively to smaller scales that are eventually dissipated by viscosity when their size becomes sufficiently small.
Atmospheric flows are characterized by very large Reynolds numbers. For a typical length scale of 1 km and a fluid velocity of 10 m s–1, the Reynolds number is 7 × 108. The atmosphere can therefore be treated as a frictionless medium (inviscid fluid), prone to multi-scale turbulent motions. This is illustrated in Box 4.2 Figure 1 (b) with large-scale turbulence in the Jovian atmosphere. The solution of the momentum equation in this turbulent regime is very sensitive to small perturbations and to inaccuracies in the initial and boundary conditions, so that the predictability of the flow (i.e., weather) is limited to a few days. When describing small-scale motions such as those encountered at the interface between two fluid elements, the Reynolds number is considerably smaller; viscous dissipation leads to mixing between adjacent fluid elements, and the full Navier–Stokes equation needs to be considered. An example of turbulent mixing of initially separated chemical tracers is shown in Box 4.2 Figure 1 (a).
Box 4.2 Figure 1 (a): Concentration of two initially separated tracers mixing and reacting in an isotropic turbulent flow of a wind tunnel. The distribution is derived by the 512 × 512 × 1024 grid point direct numerical simulation (DNS) of de Bruyn Kops et al. (2001). The image shows the wide range of length scales involved in the turbulent mixing process. See www.efluids.com/efluids/pages/gallery.htm. (b): Large-scale turbulent flow on Jupiter with irregular motions and the presence of vortices.
From the National Aeronautics and Space Administration (NASA).
Deterministic representation of turbulent flow requires a direct numerical simulation (DNS) method in which the whole range of spatial and temporal scales is resolved with a ~1 mm-resolution grid and very small time steps. An alternative approach commonly adopted for PBL turbulence is the LES method in which large-scale eddy motions are resolved explicitly while the effects of small-scale eddies are parameterized. If we separate the small-scale, rapidly fluctuating components (eddy term denoted by prime) of the dependent variables (velocity, pressure, and density) from their large-scale, slowly varying components (mean term denoted by overbar) and ignore the Coriolis term, the three components of the Reynolds averaged Navier–Stokes (RANS) equation are written under the assumption of incompressibility as
where represents the components of the acceleration resulting from gravity and other body forces. The symmetric tensor , called the Reynolds stress, characterizes the action of turbulent motions on the large-scale flow. To solve the equation, a “closure relation” that relates this eddy term to the mean flow must be specified. See Chapter 8 for more details.
For the rotating Earth, the velocity V of an air parcel at a location r and expressed in an inertial frame is equal to its velocity v relative to the Earth (the velocity that is measured by an observer located on the Earth) plus the velocity owing to the rotation of the Earth:
V = v + [Ω × r]
(4.69)
where Ω is the Earth’s angular velocity vector (directed from the south to the north pole and with an amplitude of 7.292 × 10–5 s–1). We deduce that the absolute acceleration is
(4.70)
Here 2[Ω × v] is the Coriolis acceleration, and [Ω × [Ω × r]] = –Ω2 R is the centripetal acceleration where R is the position vector perpendicular from the Earth’s axis of rotation (Figure 4.5). When expressed relative to the rotating Earth, the Navier–Stokes equation (also called the equation of motion) becomes
(4.71)
where g = ga + Ω2R is the apparent gravitational acceleration (gravitational acceleration corrected by the centripetal force) and Fdiss denotes the dissipation term resulting from viscosity. In most atmospheric applications, the effect of molecular viscosity can be ignored and the corresponding equation, in which term Fdiss is neglected, is referred to as the Euler equation. As we proceed, we will keep in the momentum equation a term F (vector with three components Fλ, Fφ, and Fz) to account for any possible momentum dissipation.
Figure 4.5 Coordinate systems: Inertial frame attached to the center of the Earth and rotating frame attached to the Earth’s surface. A point in the atmosphere is determined by its longitude λ, its latitude φ, and the distance r from the center of the Earth. The altitude z is given by r – a, if a is the Earth’s radius. The Earth rotation rate is Ω.
From Brasseur and Solomon (2005).
When expressed in spherical coordinates (λ, φ, z), the three components of the equation of motion become
(4.72)
(4.73)
(4.74)
It is important to note in the last equation (vertical projection of the momentum equation) that two terms (pressure gradient and gravity) are dominant, and this equation is often replaced by the hydrostatic approximation (see Section 2.6.1):
(4.75)
When neglecting the smaller terms and adopting the shallow atmosphere assumption (with r = a), we obtain
(4.76)
(4.77)
where f = 2 Ω sin φ is the Coriolis factor, and a is the Earth’s radius. The so-called metric terms in (4.72) t
o (4.74), uw/r, [u2 + v2]/r, uv (tan φ)/r, u2 (tan φ)/r, arise from the spherical coordinate system. A scale analysis shows that these terms can be neglected when applied to mid-latitude systems. In this case, the equation for large-scale quasi-horizontal flows is expressed in the vector form
(4.78)
where vh (u,v) represents the horizontal wind vector (i.e., wind along a surface of equal geometric height),
vh = ui + vj
(4.79)
F(Fλ, Fφ) is again the friction force, and ∇h the horizontal gradient operator:
(4.80)
Here, i, j, and k are the unity vectors in the zonal (x), meridional (y), and vertical (z) directions.
By applying the k ⋅ [∇h×] and ∇⋅ vector operators on the horizontal equations of motion, one obtains the vorticity and divergence equations, which are often used in atmospheric general circulation and weather prediction models as a replacement for the momentum equations. These are expressed by
(4.81)
and
(4.82)
where the vertical projection of the curl of the horizontal velocity
(4.83)
is called the relative vorticity and
(4.84)
is the horizontal divergence. If we define the absolute vorticity as
Modeling of Atmospheric Chemistry Page 16