(4.259)
For example, the spatial resolution representative of a model with a T63 truncation is 1.8 degrees or 210 km along a latitude circle. It is 1.1 degrees or 125 km for a model with a T106 truncation.
Spectral methods have several advantages: The derivatives are accurately determined because they are calculated analytically with no related numerical diffusion. The method does not produce aliasing of the quadratic nonlinear terms and hence no nonlinear numerical instability, the use of staggered grids is avoided, and there is no pole problem. There are also disadvantages: the calculation of nonlinear terms is computationally expensive and the number of arithmetic operations increases faster with spatial resolution than in grid point models. Therefore, spectral methods are not well suited for implementation in massively parallel computing architectures. Spectral methods are also unsuitable for regional models where boundary conditions must be specified at the edges of the modeling domain. Other disadvantages include non-conservation of mass and energy (Warner, 2011), spurious waves in the vicinity of large discontinuities (Gibbs phenomenon, see Box 4.9), and occurrence of negative values for the dependent variables.
4.10 Finite Element Method
The finite element method (Courant, 1943) is, to a certain extent, analogous to the spectral method. In this numerical technique, the solution of the partial differential equations is also expressed by a finite sum of spatially varying basis functions, but rather than being global as in the spectral method, they are expressed by low-order polynomials that are different from zero only in localized regions. Again, the PDEs are transformed into ordinary differential equations that are numerically integrated by standard techniques. To describe the method, we consider the PDE (in the 1-D space)
L[Ψ] = F(x)
(4.260)
to be solved over the model domain [a,b] with specified boundary conditions Ψ(a) and Ψ(b). Here, represents a differential operator. For example, in the case of the linear 1-D advection equation with a velocity u, this operator can be expressed by ∂/∂t + u ∂/∂x, and F(x) = 0. In the finite element method, the solution Ψ(x, t) of the PDE is approximated by a finite series of specified basis orthogonal functions Φk(x) that are non-zero only in a small part of the total domain called finite element (see below and Figure 4.23):
(4.261)
The unknown coefficients ak are a function of time t. As in the case of the spectral method, an advantage of this formulation is that the space derivative can be calculated exactly from the known basis functions.
Figure 4.23 Example of four basis functions in a finite element j as a function of the spatial dimension x.
Personal communication of Paul Ullrich, University of California Davis.
The error resulting from the fact that the approximation in (4.261) retains only K terms is
(4.262)
In the algorithm introduced by Galerkin, the residual eK is required to be orthogonal to each basis function Φj(x), so that for all values of j = 1, K
(4.263)
By substitution, we derive a system of K algebraic equations for the unknown coefficients ak(t)
(4.264)
This system of K coupled ordinary equations can be solved to obtain the coefficients ak(t). For this purpose, the time derivatives included in operator are usually approximated by finite differences. The basis functions inside each element are often expressed by a piecewise linear function called Chapeau function (Figure 4.24, a) equal to unity at the jth node and zero to all other nodes. Its mathematical formulation is
(4.265)
(4.266)
Φj(x) = 0 elsewhere.
(4.267)
An approximation of the true function is provided by a linear combination of the Chapeau functions (Figure 4.24, b). Higher-order shapes (polynomial, curvilinear elements) can also be adopted as basis functions.
Figure 4.24 (a) Chapeau function for finite element j. (b) Linear combination of basis functions (black line) to define a piecewise (in green) linear function that approximates the solution Ψ(x).
Contrary to the spectral method that describes the physical quantities by expansion functions over the entire domain, the finite-element method represents the fields for individual elements by a combination of a small number of basis functions (Figure 4.24, b). Much of the computation is thus local to a single element. In this case, the accuracy is not anymore achieved by using a large number of basis functions but by increasing the number of elements inside the entire model domain. These elements can be made smaller where a high-resolution formulation is required. The finite element technique is attractive for complex grid geometries including unstructured meshes or when the desired precision varies over the entire domain. It has been developed and often applied for aeronautic or civil engineering applications, but is not often used in chemical transport modeling. It could become an attractive approach for unstructured or adaptive grids. An application of the method to the 1-D advection equation is given in Section 7.8. The spectral element method (Canuto et al., 1984; Patera, 1984), which is used in some modern atmospheric general circulation models (Fournier et al., 2004), is a high-order finite element technique that combines the geometric flexibility of finite elements with the high accuracy of spectral methods.
4.11 Lagrangian Approaches
The methods described in the previous sections are Eulerian since the dependent variables are calculated relative to a numerical grid attached to the rotating Earth. In the Lagrangian approach, the dependent variables such as chemical concentrations are calculated following the trajectory of infinitesimally small air parcels (called particles) displaced by air motions. Different types of Lagrangian models have been developed (Figure 4.25) and we give here a brief overview.
Figure 4.25 Schematic representation illustrating different Lagrangian model formulations. (a) Trajectory of a single air particle that retains its identity as it moves along a single line determined by the mean wind velocity; turbulent mixing is neglected. This formulation is appropriate when the flow is laminar. (b) Trajectories of multiple particles aggregated in a small volume (box). Each individual particle is displaced along the flow. The deformation of the box provides information on the dispersion of the particles resulting from inhomogenieties in the flow. (c) Gaussian puff advected by the wind and affected by turbulent dispersion (see Section 4.12). (d) Particle dispersion model for a large number of individual particles with random wind velocities (turbulence) generated by a stochastic (Markov) process.
From Lin (2012).
4.11.1 Models for Single Trajectories
The simplest Lagrangian models are trajectory models that track the time evolution of a single particle along its trajectory (or of an ensemble of neighboring particles to account for errors in airflow). Particles are assumed to retain their identity during their displacement. They transport physical quantities such as momentum, energy (i.e., potential temperature), water vapor, or chemicals. Trajectory models are simple and useful, for example, to track the fate of a pollution plume, or to determine the sources contributing to concentrations observed at a receptor site.
The Lagrangian trajectory of a particle located at r(t) with three spatial components (x, y, z) is derived by integrating
(4.268)
in which v(u, v, w) is the velocity vector along the particle trajectory. If the wind field v(r(t)) and the initial position of the particle are known, the trajectory is completely determined. The “zero acceleration solution,” which is computationally inexpensive but only first-order accurate, provides the position of the particle at time t + Δt as a function of the particle at time t
r(t + Δt) = r(t) + Δtv(r(t))
(4.269)
where Δt is the integration time step. A second-order accurate solution is obtained from the so-called “constant acceleration solution” or Pettersen’s scheme
(4.270)
This implicit expression is solved by an iterative method. Expressions with higher-order approximations can also be used. One const
ructs in this manner forward trajectories or, if the sign of the time interval is reversed, backward trajectories. Wind vectors used to construct the trajectory are typically interpolated from an archive of assimilated meteorological data available on a fixed grid and time interval. Errors on the calculated trajectories arise from the truncation in the finite difference scheme, the quality of the wind fields, the interpolation of the winds, poor knowledge of the unmeasured vertical wind component, and the starting position of the trajectory (Stohl, 1998). Back trajectories constructed to analyze the history of an air parcel are often useful only for a few days, beyond which the errors grow too large, especially if the backward-moving particles encounter convective situations with large unresolved vertical motions.
Local sources and sinks may affect the chemical variables transported along trajectories. The change in mixing ratio μi of a species i along a trajectory is obtained by integrating the Lagrangian form of the continuity equation:
(4.271)
Here, Si [expressed for example in ppm s–1] represents the net source rate of species i. As the particle is transported from departure point A at time tA to arrival point B at time tB, the change in mixing ratio is obtained by integration along the trajectory
(4.272)
Si may include a source term from emission or chemical production, available as input on the same grid and time interval as the winds; and a sink from first-order loss. The source term must be interpolated in the same way as the winds, and its value along the trajectory between points A and B is often estimated at the midpoint or as the average between values at A and B. The first-order loss is applied as an exponential decay along the trajectory.
4.11.2 Stochastic Models
Lagrangian stochastic (LS) models, also called Lagrangian particle dispersion models (LPDMs), simulate the transport and dispersion of chemicals by calculating the random trajectories of an ensemble of particles in the turbulent flow. In the zeroth-order formulation, the displacement of particles is treated as a Markov chain, where the probabilities of future states do not depend on the path by which the present state was achieved. The position r of each particle is determined by a sequence of random increments. The change in the three components (x, y, z) of the position r of a fluid element is derived from the stochastic differential equations
dx = αxdt + σxdξx dy = αydt + σydξy dz = αzdt + σzdξz
(4.273)
where the terms with αx, αy, and αz represent the deterministic motions and the terms with σx, σy, and σz represent the stochastic component. Components dξx, dξy, and dξz are uncorrelated, and chosen so that their mean values equal zero and their variances equal dt. One can show (Boughton et al., 1987; Sportisse, 2010) that the mean concentration obtained from the stochastic equations satisfies the Eulerian advection–diffusion equation if
(4.274)
and
(4.275)
where is the 3-D field of the mean (resolved) wind velocity and (Kx, Ky, Kz) are the diffusion coefficients for the three directions. Under these conditions, the particle position evolves over time step Δt as
(4.276)
The three components Δξx, Δξy, Δξz, known as the Wiener–Lévy process in the theory of Brownian motion, are stochastic components with mean values equal to zero and variances equal to Δt.
In the first-order method, the particle path is derived by a sequence of random increments applied to the turbulent velocity v′ (rather than the position) of the particle. The turbulent motions are represented again by a Markov process and the variation in the wind velocity by a stochastic differential equation established by Langevin (1908) to describe Brownian motion (Thomson, 1987):
(4.277)
The deterministic part of the acceleration is defined by vector a = (ax, ay, az)T and the random component is defined by the 3 × 3 matrix B = (bij). In most cases, only the diagonal components of B (bxx, byy, bzz) are retained. In these relations, the terms dξx, dξy and dξz denote again the spatial components of a Gaussian white noise with an average of zero and a variance of dt. Expressions for these terms are available for different types of turbulence (see Wilson and Sawford, 1996).
For a stationary and horizontally homogeneous turbulent flow with a constant mean wind directed in the x-direction and a constant air density, Luhar (2012) writes the following Langevin equations for the three components of the velocity
(4.278)
where denote the variances of the three spatial components of the wind velocity. The Lagrangian timescales in the three spatial directions are
(4.279)
Here, ε denotes the turbulent kinetic energy dissipation rate. For homogeneous turbulence in the vertical direction, . In atmospheric boundary layer problems, the mean horizontal wind velocity is usually large relative to turbulent fluctuations, and the Langevin equation is written only for the vertical wind component (Stohl, 2005).
Once the velocity components have been calculated, the position (x(t), y(t), z(t)) of the particle is obtained by integration of
dx = u dt dy = v dt dz = w dt
(4.280)
Lagrangian models can be run forward or backward in time. In the forward mode, particles released at one or more source locations are transported by the fluid motions through the model domain. In the backward mode, particles released at a receptor point are used to determine the upwind influences at that point (Figure 4.26). A frequent application of the backward mode is to determine the surface flux footprint contributing to the atmospheric concentration at a given location.
Figure 4.26 Lagrangian receptor model to determine the sources (“footprint”) of chemical concentrations measured at a receptor point at a given time. A large number of Lagrangian particles are released at the receptor point and observation time, and are then tracked back in time. They diverge as a result of random turbulence, wind shear, and convection. A surface footprint for the observations can be determined from the statistics of the particles reaching the surface at all backward times.
Adapted from Lin et al. (2003).
4.11.3 Global and Regional Three-Dimensional Lagrangian Models
The Lagrangian formalism can be applied for calculating the spatial distribution and temporal evolution of chemical composition in global or regional atmospheric domains. In these models, one follows the displacement and dispersion of a large number of particles (typically 106 or more). The particles maintain their integrity as the model integration proceeds, so that particles of different origins in the spatial domain never mix and the chemical species that they carry are not allowed to react. Accounting for chemical reactions is often done on an Eulerian grid interpolated from the Lagrangian particle information.
Lagrangian models have several advantages over their Eulerian counterparts. There is no numerical diffusion so that sharp gradients are preserved. The effective resolution in regions of particular interest can be readily enhanced by increasing the density of particles. The numerical integration of the trajectory equations is stable and the time step is therefore not limited by the Courant–Friedrichs–Lewy (CFL) criterion (see Section 7.3). The independence of calculations along individual trajectories is particularly attractive for massively parallel computing architectures.
Lagrangian formulations also have several disadvantages. The lack of mixing between air parcels generates spurious small-scale features that would be smoothed in the real atmosphere by turbulent mixing. Not accounting for mixing is problematic for simulating nonlinear chemistry and aerosol processes where mixing can greatly affect rates. Lagrangian models also generally cannot provide uniformly dense coverage of a given 3-D domain; depending on shear in the flow, particles may cluster in some regions of the domain while leaving other regions unsampled.
4.11.4 Semi-Lagrangian Models
A semi-Lagrangian approach is often adopted in global models. Here, the concentrations are calculated at fixed points on an Eulerian grid. The transport over a time step Δt = tn+1 – tn is calculated by initiating back traj
ectories over time Δt from the location of each gridpoint at time tn+1. The location of the departure point, i.e., the location of the corresponding particle at the previous time tn, is derived by integrating the trajectory backward in time. The departure point does not in general correspond to a gridpoint, and thus the mixing ratio at that departure point is derived by interpolating the values from neighboring grid points at time tn. A great advantage of semi-Lagrangian methods is that they are computationally stable for any time step. Chemical reactions are computed on the model grid, allowing for mixing and nonlinear processes. However, interpolation of the concentration field compromises mass conservation and this requires correction at every time step. Semi-Lagrangian algorithms have been combined with finite volume Eulerian methods (Lin and Rood, 1996) to overcome Courant number limitations. More details on semi-Lagrangian methods are provided in Section 7.8.
4.12 Atmospheric Plume Models
The atmospheric dispersion and chemical evolution of a plume originating at a source point r0 can be represented by Gaussian plume models. Such models have the computational economy of trajectory models while accounting for small-scale eddy motions (turbulent diffusion) and allowing for nonlinear chemistry.
Modeling of Atmospheric Chemistry Page 22