(4.208)
This equation does not suffer any approximation, so the accuracy of the solution is determined by the algorithm adopted to evaluate the integral of function F over a time step Δt. Different approaches are illustrated in Figure 4.13. In the simplest of them, the function F is approximated by a constant value over the time interval, equal to its value at time tn or at time tn+1. These approximations correspond to the first-order accurate explicit Euler forward and the implicit Euler backward methods, respectively. If function F is approximated by a linear function coincident with the value of F at times tn and tn+1, the accuracy of the method is improved and corresponds to the second-order trapezoidal rule. Finally, if the integration is performed over an interval 2 Δt from time tn–1 to tn+1 and if function F is approximated over this interval by the value of F at time tn, the method is also second-order accurate and corresponds to the leapfrog algorithm.
Figure 4.13 Geometric interpretation of numerical integration algorithms for ordinary differential equations: (a) Euler forward; (b) Euler backward; (c) trapezoidal and (d) leapfrog methods.
Personal communication of Paul Ullrich, University of California Davis.
More accurate quadrature algorithms can be implemented to calculate the integral of function F(Ψ, t). We can write, for example,
(4.209)
where cj < 1 denotes nodes within the interval Δt, and bj weights used for the quadrature. The values of Ψ(tn + cjΔt) are not known, but they can be estimated numerically by a series of s preliminary calculations called stages. In the case of the s-stage explicit scheme, one computes sequentially the following quantities
(4.210)
for j = 2, 3, …, s and with Φ1 = Ψn. Parameters s and aj,k are specific to a particular scheme. The solution at time tn+1 is then approximated by
(4.211)
When considering implicit s-stage schemes, the upper bound (j – 1) of the summation in expression (4.210) is replaced by a value equal or larger than j, for example by parameter s. Runge–Kutta and Rosenbrock methods are examples of multi-stage approaches and are discussed in Chapter 6.
Partial Differential Equations (PDEs)
We now illustrate the use of finite difference methods in the case of partial differential equations. We consider here the 1-D linear advection equation
(4.212)
where the scalar field Ψ(x, t) represents, for example, the mixing ratio of a chemical species and c a constant velocity of the flow in the x direction. The initial condition is provided by a specified function of space x: Ψ(x, 0) = Ψ0(x). As will be shown in Section 7.3, the explicit forward in time, centered in space (FTCS) Euler solution
(4.213)
is unstable for any time step. Here, represents the value of function Ψ(x, t) at discrete time tn and geometric location xi. Other numerical algorithms that are at least conditionally stable (stable for a sufficiently small value of the time step) are presented in Chapter 7.
In the case of the combined 1-D advection–diffusion equation (with K representing a diffusion coefficient):
(4.214)
the solution provided by the FTCS Euler forward scheme
(4.215)
is stable for sufficiently small time steps through the addition of the diffusion term. The Lax–Wendroff numerical scheme makes use of this property when applied to the pure advection equation: A numerical diffusion term proportional to the chosen time step is added to the FTCS scheme (see Section 7.3).
If an implicit time stepping scheme is adopted for the advection–diffusion equation
(4.216)
the solution is unconditionally stable (stable for all values of the velocity c and the diffusivity K). The system, however, is more difficult to solve. By rearranging the terms, (4.216) can be written as
(4.217)
where
(4.218)
(4.219)
(4.220)
The system of I equations, whose matrix is tridiagonal, can be solved for each time tn+1 by a matrix decomposition and back-substitution method (Box 4.4). The approach described here can be generalized to more than one spatial dimension. Again, linear algebraic equations are derived and can be solved either directly in the explicit case or through more complex methods when an implicit algorithm is used. A more detailed description of the methods used to solve the transport equations is given in Chapter 7.
Box 4.4 Tridiagonal Matrix Algorithm
The Thomas algorithm, a simple form of Gaussian elimination, can be used to solve the tridiagonal system of equations that arises from the discretization of the 1-D diffusion equation
αiΨi − 1 + βiΨi + γiΨi + 1 = δi (i = 1, 2, …., n)
where α1 = 0 and γn = 0. If we first define the modified coefficients
the solution is obtained by back substitution
The fluid dynamics equations, specifically the momentum equation, include nonlinear terms of the form u ∂u/∂x. The corresponding equation, expressed here along spatial dimension x,
(4.221)
produces interactions between atmospheric waves. As a result, new wave modes are generated, including waves shorter than can be explicitly represented by the model grid. The energy of these short waves is then folded back into longer waves, so that wave energy tends to spuriously accumulate in the spectral region near the cut-off wavelength. This process, called aliasing, can be a source of nonlinear numerical instability in dynamical models. Short waves can be eliminated by a filtering process (see Section 4.15).
4.8.2 Finite Volume (Grid Cell) Methods
The finite volume method provides another approach for solving PDEs. It is particularly well suited for applications in which mass and energy conservation is a critical consideration. In this method, the prognostic quantities are not defined on discrete nodes i (grid points), but are instead expressed as averages across specified finite control volumes. These control volumes correspond to the specified grid cells of the model. In the simple 1-D formulation, where the volume of cell i is replaced by the size of the spatial mesh Δxi = xi+½ – xi–½, the average ⟨Ψi (t)⟩ of a variable Ψ over cell i is given by
(4.222)
Consider, for example, a 1-D version of the continuity equation (similar to (4.1) but for a single geometrical dimension x, and for no chemical source):
(4.223)
where Ψ represents the mass or number density of a chemical species and F its flux in the x direction. The value of Ψ at time tn+1 = tn + Δt is obtained from the value at time tn by
(4.224)
The equation for the average of Ψ over grid cell i (denoted ⟨Ψi⟩) is
(4.225)
or, after integration over Δxi
(4.226)
The time derivative of this last expression is
(4.227)
where the value of flux F is estimated at both edges of grid cell i. Note that no approximation has been made to derive expression (4.227).
If (4.223) represents the continuity equation for density Ψ and the flux F = Ψ u is expressed as the product of this density by the velocity u in the x direction, (4.227) can be rewritten as
(4.228)
If the velocities are prescribed at the edges of the grid cells, and the corresponding values of Ψ are constructed by interpolation of the mean values in the neighboring cells, then (4.228) can be expressed as
(4.229)
This first-order differential equation can easily be integrated numerically (see example in Section 4.8.1).
An important advantage of this method (compared to the finite difference method) is that it does not require a structured geometry, but can easily be applied to different types of grids (see Section 4.8.3). Consider the more general case of the hyperbolic problem
(4.230)
where the flux vector F has components in all spatial dimensions. The spatial domain under consideration can be subdivided into finite volumes Vi, and the equation can be integrated over the volume of a cell i
(4.231)
&n
bsp; Noting that the first term is equal to Vi d⟨Ψi⟩/dt and applying the divergence (Stokes) theorem to the second term, we find
(4.232)
Here, the integral is taken along the walls Si of cell i and represents the material flowing across the boundaries of the cell; n is the unit vector perpendicular to wall element dS. Again, the values of F at the edge of the cell can be constructed by interpolation of average values in neighboring cells. As can be deduced from the equations, the finite volume method conserves the transported variables easily even for coarse grids: What is lost in one cell is gained by a neighboring cell. This approach is therefore well adapted to treat the advective transport of tracers on complex grids (see Chapter 7).
4.8.3 Model Grids
Different types of grids can be adopted to solve the finite difference or finite volume approximations to the model equations. The choice of grid must balance considerations regarding the nature of the problem to be solved, the need for accuracy and stability, and the computational resources at hand. Another important consideration in chemical transport models using offline meteorological fields is to ensure consistency with the grid of the parent meteorological model and thus to avoid interpolation errors.
Vertical Discretization
We first consider the vertical discretization of the continuity equation for scalar quantity Ψ (here again a mass or number density), and write, for example, for the vertical component
(4.233)
where z is the vertical coordinate adopted in the model (not necessarily the geometric altitude). If quantity Ψ and the vertical wind component w are defined on the same discrete levels (j – 1, j, j + 1), (4.233) can be approximated as
(4.234)
If, on the other hand, the vertical layers are staggered (Box 4.5) as shown in Figure 4.14, and if quantity Ψ is defined at levels j – 1, j, and j + 1, while the vertical wind component w is specified at levels, j – ½ and j + ½, the tendency resulting from the vertical component of the flux convergence is approximated by
(4.235)
or, if function Ψ is linearly interpolated,
(4.236)
Box 4.5 Grid Staggering
When applying a numerical algorithm to solve a system of differential equations, all unknown variables are not necessarily defined at the same grid points. In a staggered grid model, the different dependent variables are provided on different grid points, usually offset by half the grid size (see, for example, Figure 4.14 in the case of vertical grid staggering). Staggered grids are often adopted because the accuracy of the calculated spatial derivatives is improved. Further, the short wavelength components of the solution are often more accurately represented (Durran, 2010). In an unstaggered grid, the derivative Gi at location xi of a function Ψ is often calculated across a 2Δx interval
while, in a staggered grid, it can be derived across a single Δx interval
In the latter case, the derivative G is calculated on one grid (i – 1, i, i + 1) while function Ψ is defined on the other grid (i – 1/2, i + 1/2). Staggered grids have the advantage of effectively halving the grid size, so that the truncation errors are reduced. A standard configuration is to define wind components on one grid and other dependent variables such as pressure, temperature, concentrations, etc. on the other grid. The utility of this can directly be seen in the finite-volume approximation (4.229) of the advection equation. Since the effective grid size is smaller than the nominal size Δx, the time step Δt adopted for integrating the equations must be correspondingly reduced to maintain the stability of the solution.
Figure 4.14 Example of a staggered vertical grid: temperature T, horizontal wind components (u, v), and function Ψ (e.g., trace species density ρi) are derived at the center of the cell (midpoint level) while the vertical wind component w and the eddy diffusion coefficient K are derived at the edge of the cell (interface level).
In general circulation models of the atmosphere, if the temperature and the vertical component of the wind are defined on intermediate levels between the geopotential and the horizontal wind components, the vertical staggering refers to the so-called Charney–Philips grid. If the geopotential, the temperature, and the horizontal wind components are provided on the same levels, while only the vertical wind component is defined on intermediate levels, the staggering is referred to as the Lorenz grid.
Horizontal Discretization
As the PDEs of the model are transformed into their finite difference or finite volume analogs, the variables at each level of the model have to be defined on a finite number of grid points or grid cells. Here again, different approaches can be used. The simplest of them is to define all variables on the same grid points and implement the discretization procedure accordingly. Such a grid is called A-grid by Arakawa and Lamb (1977). Despite its simplicity, it is little used because it often produces noise in the numerical solution of the equations. Alternate staggered configurations are shown in Figure 4.15. The C-grid is often adopted in atmospheric models because it is readily adapted to the finite-volume solution (4.229) of the continuity equation. A more detailed discussion of the different staggered grids is provided by Arakawa and Lamb (1977) and by Haltiner and Williams (1980).
Figure 4.15 Unstaggered grid (A) and different staggered grids (B, C, and D) in two dimensions, using the Arakawa and Lamb (1977) classification. ψ is a scalar simulated by the model, for example chemical concentration; u and v are the zonal and meridional components of the wind velocity.
Grid Geometries
The longitude–latitude grid is frequently adopted to represent the horizontal distribution of atmospheric variables (Figure 4.16). It has the advantage of being orthogonal so that the finite difference or finite volume approximations of the PDEs can be easily derived. A disadvantage is that the geometric spacing in longitude Δx = a Δλ cos φ decreases gradually from the Equator to the pole. (Here λ represents the longitude, φ the latitude, and a the Earth’s radius). As discussed in Section 4.7.1, the time step Δt to be adopted in the model to ensure stability of the integration method is often bounded by the smallest geometric grid spacing. In the case of the advection equation, the time step applied in Eulerian schemes must be chosen in most cases to be smaller than the so-called Courant number, which is proportional to the grid size Δx (see Chapter 7 for definition and for examples). Thus, the existence of singular points at the poles introduces severe limitations on the adopted time step and leads to low computational performance.
Figure 4.16 Latitude–longitude grid with representation of model levels and grid cells.
Different approaches have been proposed to avoid the pole problem. The most straightforward is to apply either spatial or Fourier filters to eliminate the smaller-scale features (noise) appearing in the solutions near the poles. A second approach is to reduce the number of points in the longitudinal direction and hence to increase the longitudinal grid spacing in the vicinity of the poles. In Figure 4.17, a reduced grid is compared to a regular longitude–latitude grid. A third solution is to apply in the polar regions some backup algorithm that is not constrained by the Courant criteria (e.g., semi-Lagrangian methods, see Section 7.8).
Figure 4.17 (a): Latitude–longitude grid. The convergence of the longitudinal lines toward the poles reduces gradually the geographical distances in the longitudinal direction, which in turn limits the time step required to solve the partial differential equations. (b): Reduced grid in which the number of longitudes per circle of latitude decreases as one approaches the poles.
From Washington et al. (2009).
Other types of grids have been proposed to address the pole problem of the longitude–latitude grids. An increasingly frequent approach adopted in conjunction with finite volume methods (see Section 4.8.2) and referred to as “tiling” is to cover the surface of the sphere by some geometric shapes with no overlaps and no gaps. In the Voronoi tessellation approach, the surface of the Earth is partitioned into closed regions around pre-defined points, called seeds or generators of
the grid. By definition, a Voronoi cell for a given seed includes all geometric points that are closer to that particular point than to any other seed points in the model domain. The seed points can be specified so that the cells cover rather homogeneous areas of the surface (such as uniform ecosystems). They can also be chosen to have regular polygonal shapes (Figures 4.18 and 4.19). The icosahedral grid, for example, is characterized by a relatively uniform spatial resolution without any singularity. The grid elements can be triangular or hexagonal (with in this particular case 12 pentagonal cells). The cubed sphere grid also allows a relatively uniform resolution. In this approach, the existence of eight special “corner points” and plane boundaries requires special treatment. The Yin–Yang grid is the combination of two distinct longitude–latitude grid systems with mutually orthogonal axes and partial overlap (contact region). The advantage of homogeneous and highly isotropic model grids is that the spatial density of model results is relatively uniform, so that the data can be efficiently archived, analyzed, and remapped to other grid structures.
Modeling of Atmospheric Chemistry Page 20