(4.112)
and the vertical velocity ω in the pressure coordinate system is defined as
(4.113)
For any arbitrary vertical coordinate η, the “horizontal” gradient ∇ηΨ of a function Ψ is related to the horizontal gradient of this function ∇z Ψ expressed in the geometric coordinate framework z by (see (4.110))
(4.114)
We deduce when η = p and when applying the hydrostatic approximation that
(4.115)
since ∇p p = 0. The continuity equation expressed in pressure coordinates takes a particularly elegant expression. Separating in the continuity equation (4.59) the contributions provided by the horizontal (vz) and vertical (w) winds
(4.116)
and using the hydrostatic equilibrium relation to replace ρa by (–1/g ∂p/∂z) we find the remarkable result
(4.117)
where vp (u, v) is the “horizontal” wind (components along isobaric surfaces) and ω = dp/dt is the previously defined vertical velocity in the pressure coordinate system. This form of the continuity equation contains no reference to the air density and is purely diagnostic (no time derivative). Its simplicity is one of the advantages of the pressure coordinate system. In spherical coordinates, the continuity equation is written as
(4.118)
By integrating (4.117), we find the isobaric tendency equation that provides the vertical velocity ω at any pressure level p by integrating the isobaric divergence of the horizontal wind:
(4.119)
In pressure coordinates, the state of the atmosphere is often represented by the geopotential Φ on isobaric surfaces (rather than the pressure on equal altitude levels). We note that
(4.120)
In pressure coordinates, the momentum equation is expressed as
(4.121)
or in spherical coordinates
(4.122)
(4.123)
where u and v represent here the wind along isobaric surfaces.
The general form of the energy equation
(4.124)
remains unchanged in pressure coordinates. In this and in the momentum equation, the total derivative of a function Ψ under the shallow atmosphere approximation is expressed as
(4.125)
where vp(u, v) represents the “horizontal” velocity vector on an isobaric surface. In spherical coordinates, this expression becomes
(4.126)
Here, the derivatives ∂/∂λ and ∂/∂φ are calculated along isobaric surfaces rather than along a constant altitude surface.
4.6.2 Log-Pressure Altitude Coordinate System
The dynamical equations can also be simplified in the log-pressure coordinate system by defining a log-pressure altitude (Chapter 2)
(4.127)
where ℋ is a constant effective scale height. This vertical coordinate is often used for stratospheric models. We have
(4.128)
where H(z) = kT/mg is the atmospheric scale height (m). The total time derivative of function Ψ is now expressed as
(4.129)
or
(4.130)
where vz is the wind vector on a constant pressure–altitude surface and
(4.131)
Components u and v stand here for the wind components along isobaric surfaces or equivalently for given pressure altitude (Z) surfaces. By choosing ℋ to be equal to 7 km, the value of Z (which corresponds to a given pressure level) is approximately equal to the value of the geometric altitude z.
When the log-pressure representation is used, the hydrostatic equation takes the form
(4.132)
and the continuity equation becomes
(4.133)
Newton’s second law is expressed as in (4.122) and (4.123).
4.6.3 Terrain-Following Coordinate Systems
A difficulty in using geometric height or pressure as the vertical coordinate is the treatment of surface topography. To address this issue, one can introduce a normalized pressure coordinate (commonly called sigma coordinate)
(4.134)
where ps is the pressure at the Earth’s surface (which varies with the topography) and ptop is the pressure at the highest level of the model. The sigma coordinate is illustrated in Figure 4.6. The value of σ varies from zero at the top of the model (p = ptop) to unity at the surface (p = ps), where a simple boundary condition dσ/dt = 0 is applied. In the original definition of the sigma coordinate by Phillips (1957), the top of the model was assumed to be the top of the atmosphere (ptop = 0) with therefore σ = p/ps. The advantage of the sigma coordinate system is that it conforms to the natural terrain and therefore eliminates the problem of intersection with the ground when the terrain is not flat. It is particularly well suited for the boundary layer.
Figure 4.6 Sigma coordinate levels above a region with variable topography.
Courtesy of Martin Schultz, Forschungszentrum Jülich.
In the sigma coordinate system, the total derivative of a function Ψ is expressed as
(4.135)
where vσ is the “horizontal” velocity along a sigma surface and dσ/dt is the vertical velocity. Alternatively, we can write
(4.136)
where the derivatives ∂Ψ/∂λ and ∂Ψ/∂φ are taken along a constant sigma surface. If π = ps – ptop, the horizontal equations of motion become
(4.137)
the hydrostatic equation
(4.138)
the continuity equation
(4.139)
and the thermodynamic energy equation
(4.140)
Even though sigma surfaces do not intersect the ground, the use of the sigma coordinate system requires some precautions. Large errors can occur in calculating pressure gradients in regions of complex topography because the constant-sigma surfaces are steeply sloped. To address this problem, Mesinger (1984) introduced the step-mountain coordinate, commonly called eta coordinate, as
(4.141)
where pref(z) is a reference pressure defined as a function of the geometric height z (e.g., pressure in the standard atmosphere with pref (0) = 1013 hPa), and zs is the local terrain elevation. The scaling factor applied to the sigma coordinate ensures that the η surfaces are quasi horizontal.
The influence of topography on the mean flow decreases with altitude, so that the sigma-coordinate system is less desirable in the upper troposphere and above the tropopause. Models frequently use hybrid σ–p coordinate systems (Figure 4.7) that follow terrain in the lower troposphere and transition gradually to follow pressure in the stratosphere. The pressure pk at model vertical level k (k = 1, K) is given by
pk = Akp0 + Bkps
(4.142)
where coefficients Ak and Bk have values that depend only on k. Parameter p0 is chosen to be equal to the pressure at sea level (1013 hPa). The surface pressure ps varies along the topography of the Earth’s surface and may also vary with time. At the surface (k = 1), A1 = 0 and B1 =1, while at the top of the model domain (k = K and pK = ptop), one imposes AK = ptop/p0 and BK = 0. The value of Bk (equal to the value of σ at the surface) decreases with height, typically down to zero in the upper troposphere or in the stratosphere.
Figure 4.7 Representation of a hybrid sigma–pressure vertical coordinate system with coefficients Ak = [0.0, 0.0, 0.0, 0.02, 0.1, 0.15, 0.18, 0.16, 0.14, 0.1, 0.05] and Bk = [1.0, 0.95, 0.90, 0.80, 0.70, 0.45, 0.25, 0.12, 0.04, 0.01, 0.0, 0.0].
Courtesy of Martin Schultz, Forschungszentrum Jülich.
4.6.4 Isentropic Coordinate System
In models that focus on the middle atmosphere, the vertical coordinate is sometimes chosen to be the potential temperature θ (Figure 4.8). The advantage is that, under adiabatic conditions, the flow follows isentropic surfaces and is therefore simple to analyze. In this isentropic coordinate system, the total derivative of a function Ψ is written as
(4.143)
where vθ (u,v) is the “horizontal” wind vector on an isentropic surface, and in spherical coordinates
(
4.144)
In this case the “vertical velocity” dθ/dt is directly proportional to the net diabatic heating rate Q. The “horizontal” momentum equations are expressed as
(4.145)
(4.146)
where Ψ = cpT + Φ is now the Montgomery streamfunction.
Figure 4.8 Isentropic coordinate system.
From the National Oceanic and Atmospheric Administration (NOAA).
The hydrostatic relation takes the form
(4.147)
and the continuity equation is written as
(4.148)
where
(4.149)
is the pseudo-density.
Air motions under diabatic conditions are “horizontal” in the isentropic coordinate system. The surfaces of constant species mixing ratio (isopleths) tend to align themselves with the isentropic surfaces and the 3-D advection becomes essentially a 2-D problem. Isentropic coordinates are therefore convenient to analyze tracer motions over timescales of 1–2 weeks since air parcels remain close to their isentropes over such a period of time. Numerical noise is also reduced.
Although attractive for the stratosphere, the isentropic coordinate system is problematic in the lower troposphere because the flow is diabatic, θ is not a monotonous declining function of z, and the isentropes intersect the Earth’s surface (see Figure 4.8). A hybrid system with isentropic coordinate above the boundary layer and sigma coordinate near the Earth’s surface is sometimes adopted.
4.7 Lower-Dimensional Models
The computational cost of a model can be decreased considerably by reducing its dimensionality. Although a 3-D representation of the atmosphere is most realistic and often necessary, there are many cases where simpler 2-D (zonal mean), 1-D(column), and 0-D (box) models can provide valuable insights (Figure 4.9). We discuss here the foundations of these lower-dimensional models and their applications.
Figure 4.9 Conceptual representation of 0-D, 1-D, 2-D, and 3-D models.
From Irina Sokolik, personal communication.
4.7.1 Two-Dimensional Models
Global 2-D (latitude–altitude) models of the middle and upper atmosphere have been used extensively to simulate the meridional distribution of ozone and other chemical species. The motivation for their use is that zonal (longitudinal) gradients are generally weak, so that resolving the longitudinal dimension would add little information. In 2-D models the dependent variables Ψ such as concentrations, temperature, wind velocity, etc. are separated according to their zonally averaged value
(4.150)
and the departure Ψ′(λ, φ, z, t) from this mean value. Thus
(4.151)
By introducing this type of variable separation into the equations presented in Sections 4.5 and 4.6, we obtain the following continuity equation
(4.152)
the zonal mean momentum equation
(4.153)
the zonal mean thermodynamic equation
(4.154)
and the zonal mean continuity equation for chemical species
(4.155)
Here, ρ0(Z) is a standard vertical profile of the air density, Z the log-pressure altitude, a the Earth’s radius, f the Coriolis factor, Q the net heating rate, and S the chemical source term. The correlation (eddy) terms such as or represent the effects of atmospheric waves on the zonal mean quantities. Closure relations expressing the eddy terms as a function of the mean quantities must be added to this system. In most cases, the eddy terms are parameterized as a function of empirical eddy diffusion coefficients. One can show that, for steady and conservative waves, the eddy flux divergence and the zonal mean advection terms cancel (non-transport theorem), so that errors in the specified eddy diffusion coefficients can have a large effect on the solution of the system, and specifically on the calculated distribution of chemical species.
To avoid the numerical problems associated with the quasi eddy-mean flow cancellation (small difference between two relatively large terms), it is useful to introduce the transformed Eulerian mean (TEM) velocities v* and w* (Boyd, 1976; Andrews and McIntyre, 1976):
(4.156)
(4.157)
It can be shown that, in the presence of conservative and steady waves, the TEM meridional velocity components are equal to zero. The quasi-compensation between large zonal mean and eddy terms that characterizes the classic Eulerian mean equations is therefore replaced by the concept of residual velocities that describes a meridional residual circulation. In the TEM framework, the continuity, momentum, and thermodynamic equations are expressed as
(4.158)
(4.159)
(4.160)
where
(4.161)
represents the diabatic heating term expressed in K per unit time. With the TEM transformation, all eddy contributions are included in the forcing terms and . Wave momentum forcing is assumed to be proportional to the divergence of the Eliassen–Palm (EP) flux E
(4.162)
For planetary waves under quasi-geostrophic scaling, the meridional and vertical components of this vector are proportional to the eddy momentum and heat fluxes, respectively. With a good approximation,
(4.163)
(4.164)
The divergence is regarded as the wave stress that accelerates or decelerates the mean flow. It vanishes when the waves are steady and conservative. The contribution of eddies to the mean temperature tendency is given by
(4.165)
This term vanishes exactly for steady and conservative waves (Andrews and McIntyre, 1978), and is otherwise generally small and often ignored. In this case, to a good approximation, the vertical TEM velocity is directly proportional to the diabatic heating rate
(4.166)
In the presence of gravity waves, the eddy term reduces to
(4.167)
Finally, the zonally averaged continuity equation for chemical species is written in the TEM framework as
(4.168)
where
(4.169)
where the net eddy flux components and in the TEM framework are expressed as a function of the zonally average eddy mass and heat fluxes by (Garcia and Solomon, 1983; Andrews et al. 1987)
Again, the term vanishes if the waves are steady and conservative, and the tracer has no sources/sinks. Otherwise, it is generally small, and accounts for chemical eddy transport, i.e., the net meridional and vertical transport that occurs when reactive gases displaced by atmospheric waves encounter different photochemical environments. In this case, the eddy correlations and must be parameterized by closure relations, for example by diffusion terms
(4.170)
(4.171)
where the coefficients Kij can be expressed as a function of the chemical lifetimes of the trace gases and of the time needed for an air parcel to move through a dynamical disturbance.
4.7.2 One-Dimensional Models
One-dimensional (1-D) models are valuable conceptual tools for vertical transport and chemistry in an atmospheric column. Uniformity is assumed in the horizontal plane; in other words, horizontal flux divergence is assumed to be negligible. Vertical transport is computed using an eddy (turbulent) diffusion parameterization. One-dimensional models may be designed for limited domains (such as the continental boundary layer over a source region) or for the global domain (vertical transport and chemistry in the atmospheric column). By averaging the continuity equation and denoting horizontal means as brackets ⟨ ⟩, we obtain the following equation for the density 〈ρi〉 of a chemical species i:
(4.172)
where the net vertical flux is parameterized as
(4.173)
Here, ρa is the mean air density and Kz is the vertical eddy diffusion coefficient. This coefficient accounts for the effects of motions at all scales. The vertical flux is assumed to be proportional to the vertical gradient in the mean mixing ratio. Values of Kz are specified as input. In 1-D models for the boundary layer, semi-empirical formulations of Kz are available as a function of the surface fluxes of momentum
and heat that force a diurnal cycle of mixed layer growth and decay (see Seinfeld and Pandis, 2006). In global 1-D models, empirical values of Kz are derived from observations of chemical tracers (see, for example, Liu et al., 1982).
In the thermosphere, vertical exchanges are dominated by molecular diffusion rather than eddy mixing. The 1-D continuity equation is then formulated in terms of the zonally and meridionally averaged number density ⟨ni⟩ as
(4.174)
Here the vertical flux ⟨Φi⟩ is expressed as
(4.175)
where Di is the molecular diffusion coefficient, αT is the thermal diffusion factor (–0.40 for helium, –0.25 for hydrogen, 0 for heavier species) and Hi = kT/mig is the scale height of species i with molecular or atomic mass mi. The other symbols have their usual meaning. An approximate expression for the molecular diffusion coefficient [here in cm2 s–1] is provided by Banks and Kockarts (1973):
(4.176)
where mi and the mean molecular mass m are expressed in atomic mass units (amu), the total air number density na is expressed in [molecules cm–3], and T is in [K]. The value of the molecular diffusion coefficient increases rapidly with height as the air density decreases, which can lead to some numerical difficulties when solving the continuity equations.
Modeling of Atmospheric Chemistry Page 18