Figure 4.18 Examples of grids: icosahedral triangular (a), icosahedral hexagonal (b) and cubed sphere (c).
Source: David Donofrio, Lawrence Berkeley National Laboratory (a and b).
Figure 4.19 The Yin–Yang grid system.
Reproduced From Qaddouri (2008).
The choice of grid resolution in a model must be adapted to the scales of the processes that one wishes to resolve. Because of computational limitations, the use of a very high-resolution grid system covering the entire model domain may not be feasible; using a grid with fine resolution in areas that are presumed to be of interest, with coarser resolution elsewhere, may be an appropriate approach to address a given problem. Mesh refinement provides high resolution in selected areas of the model domain by embedding a higher resolution grid that resolves small-scale processes into a coarser grid that captures the large-scale features (Hubbard, 2002). A succession of nested grids (zoom capability) can be implemented to derive increasingly high-resolution features as one approaches a given point of the domain. The nesting can be one-way, with the coarser grid providing dynamic boundary conditions to the finer grid but no reverse effect of the finer grid on the coarser grid. It can be two-way, with full exchange of information between the two grids. Two-way nesting is far more difficult to implement because of numerical noise at the interface between grids. Figure 4.20 (a and b) shows two different approaches used in the development of multi-resolution models. Unstructured grids (Figure 4.20, c) in which the grid meshes (often triangles or tetrahedra) are distributed as irregular patterns, can also be adopted to enhance the spatial resolution in specified areas. They are often used in oceanic models. The location of the region of interest may also change with time, as in the case of the long-range transport of a pollution plume. A solution is to use a dynamic adaptive grid, in which the location of the nodes is constantly modified during the model integration to achieve high spatial resolution in areas where the gradients in the calculated fields are large (Odman et al., 1997; Srivastava et al., 2000). Adaptative grids used in atmospheric dynamical modeling (e.g., Dietachmayer and Droegemeier, 1992; Skamarock and Klemp, 1993) can also be adopted in chemical transport models to better represent the effects of multi-scale sources (Tomlin et al., 1997; Box 4.6).
Figure 4.20 Model grids using (a) a cubed sphere with successive zooming capabilities for the Western USA; (b) a icosahedral system with a stretched grid and gradual zooming capability, and (c) an unstructured grid used here for tsunami simulations.
From the National Center for Atmospheric Research (NCAR).
Box 4.6 Dynamically Adaptive Grids in Chemical Transport Models (Srivastava et al., 2000; Odman et al., 2002; Garcia-Menendez et al., 2010)
Large-scale models with limited spatial resolution may need to resolve small-scale processes in regions of high concentrations and strong concentration gradients. Static nested grids can provide high resolution over specific regions of interest (such as an urban center) but do not adjust to changing locations of these regions (as for long-range transport of a pollution plume). In the dynamic grid adaptation method, the grid resolution changes automatically at each time step to capture and follow small-scale features of interest such as plumes or frontal boundaries. The overall structure of the grid and the total number of grid points are fixed, but the locations of the grid points evolve through the simulation according to a user-defined weight function that determines in which areas grid nodes must be clustered. At a given grid point j, the weight function wj can be expressed, for example, by
where k is an index for the chemical species, nk is the concentration, and ∇2 is the Laplacian operator which measures the curvature in the concentration field. Resolution requirements of different species can be accounted through the choice of the different coefficients αk (Odman et al., 2002). The new spatial position of grid point i in a horizontal plane is calculated from
where vector Pj (j = 1 to 4) represents the original position (before movement) of the centroids of the four grid cells that share grid point i (see Box 4.6 Figure 1, left panel). After the repositioning of the grid points, the concentration must be interpolated on the displaced grid nodes, and other parameters such as the emissions and meteorological variables must be redistributed on the adapted grid.
Box 4.6 Figure 1 Development of a dynamically adaptive grid to represent the evolution of a plume.
Sources: M. T. Odman, private communication; Steyn and Rao (2010).
The two panels in Box 4.6 Figure 1 on the right show the spatial distribution of concentrations originating from a point source as simulated by a standard air quality model (CMAQ) at 1.33 km horizontal resolution (middle panel) and by the same model, but with a dynamically adapted grid (right panel). Sharper gradients with higher concentration peaks are obtained when the model grid is dynamically adapted.
Consistency Between Vertical and Horizontal Resolutions
To adequately resolve sloping features in atmospheric models including fronts and slantwise convective systems, and correctly represent the 3-D transport of chemical species associated with such dynamical systems, the grid increments in the vertical and horizontal directions should not be specified independently. Rather, the ratio between vertical and horizontal grid point spacings should be of the order of typically 0.005–0.02 m m–1 for an appropriate representation of mesoscale features (Warner, 2011). Thus, for a horizontal resolution of 100 km, typical of global models, the spacing between model layers should be of the order of 0.5–2 km. For a horizontal resolution of only 10 km, as often adopted in regional (limited-area) models, the vertical spacing should be reduced to 50–200 m. The lack of consistency between horizontal and vertical resolution in dynamical models can generate spurious waves during the simulation and thus undesired noisy fields.
4.9 Spectral Methods
An alternative approach to the finite difference and finite volume algorithms for global models is the spectral method developed by Orszag (1970) and Eliasen et al. (1970), and implemented for the first time by Bourke (1974). In this method, the horizontal distribution of the atmospheric variables is represented by a finite expansion of periodic functions (waves), called basis functions. The orthogonality of these functions (see Box 4.10) allows the derivation of coupled ODEs for the expansion coefficients, which vary with time and height. This method has been important for climate modeling but is rarely used in chemical transport models because chemical concentrations are strongly affected by local forcing processes that cannot be easily described by waves. In addition, spectral transport algorithms can produce negative concentrations.
To introduce the spectral method, we first consider a function Ψ(x) defined on the spatial interval [–π, +π]. By applying a Fourier transform and assuming that the function repeats itself with a period 2π, we can write
(4.237)
where k represents the wavenumber. Expansion coefficients ak and bk are given by the following integrals over the [–π, +π] interval
(4.238)
(4.239)
(4.240)
Function Ψ(x) may also be expressed as the sum of sine or cosine functions of different amplitudes dk, wavenumbers k, and phases φk (Box 4.7):
(4.241)
The component with the lowest frequency is called the fundamental, and any higher-frequency parts are referred to as harmonic components.
Box 4.7 An Example of Fourier Decomposition
We present here a simple illustration of the spectral decomposition method. A function y(t) equal to 1 on interval [0, 0.5] and to –1 on interval [0.5, 1], as shown in Box 4.7 Figure 1, can be approximated by
The Fourier coefficients with an even numbered index disappear in this particular case.
Box 4.7 Figure 1 Fourier decomposition of a square function ys. The fundamental y1 = 3/π sin 2πt, and the two harmonics y3 = 4/(3π) sin 6πt and y5 = 4/(5π)sin (10πt) are shown together with the sum y′ = y1 + y2 + y3 (dotted line).
Reproduced with permission From Cruse (2006).
&
nbsp; An equivalent, but often more convenient, form of the Fourier series is given by
(4.242)
with
(4.243)
This coefficient ck is related to ak and bk in (4.237) by
(4.244)
(4.245)
with c0 = a0. In practical applications, the number of terms retained in the Fourier series is necessarily limited. A truncation is therefore applied when index k reaches a specified value K, and function Ψ(x) is now approximated as
(4.246)
or
(4.247)
Accuracy increases for higher values of K but so do computing costs.
The calculation of an approximation for the spatial derivative is straightforward:
(4.248)
or
(4.249)
with coefficients
αk = k bk βk = − k ak γk = ik ck
(4.250)
The calculation of the derivatives of a function is thus easy to perform once the Fourier expansion coefficients of the function have been derived. Fast Fourier transform (FFT) algorithms are routinely used in general circulation models to transfer information between the grid space and the spectral space (see Appendix E).
To solve PDEs such as the advection–diffusion equation (Box 4.8), the solution is approximated by an expansion such as (in one dimension)
(4.251)
in which the coefficients ak are only a function of time t. When expression (4.251) is introduced in a PDE, one obtains for each value of k an ODE for ak that is solved by a classic finite difference method. The algorithm requires that K differential equations be solved consecutively, which is usually computationally less expensive than solving the equation by a finite difference method at J grid points. Indeed, in most applications where function Ψ(x, t) is relatively smooth, the value of K can be considerably smaller than the value of J. When function Ψ(x) contains discontinuities or if too few terms are included in the Fourier series (too small value of K), the solution provided by the spectral decomposition method is characterized by large overshoots and undershoots (Box 4.9), and can produce undesired negative values. This is why the spectral approach does not provide satisfying results when applied to the transport of atmospheric species with strong gradients.
Box 4.8 Fourier Transform of the Advection–Diffusion Equation
In the case of the 1-D advection–diffusion equation
with the solution approximated by (4.251), we obtain the differential equations
which are solved by an appropriate finite difference algorithm.
Box 4.9 Gibbs Phenomenon
The Gibbs phenomenon, named after J. Willard Gibbs, describes the peculiar manner in which the Fourier series of a piecewise continuously differentiable periodic function behaves at a jump discontinuity. The use of truncated series instead of infinite Fourier series leads to overshooting and undershooting of the true function. This effect can be illustrated by considering a step function whose value Ψ(x) is –1 for –π < x < 0 and 1 for 0 < x < π. This function can be expressed by the infinite series of periodic functions
If the function Ψ(x) is approximated by only the first term in the series, the maximum and minimum values of Ψ are ±4/π (or ±1.27), instead of 1 and –1, respectively. This misrepresentation of the field with overshoot and undershoot is called the Gibbs phenomenon. This effect decreases when more terms are included in the truncated series, and is less acute when the field contains no discontinuities.
Box 4.10 Orthogonal Functions
Two functions Ψn and Ψm are said to be orthogonal if the integral of their product is zero. Thus, if δnm is the Kronecker delta (equal to zero if m ≠ n and one if n = m) and A is a normalization constant
For example, when representing function Ψ(x) by a series of elementary trigonometric functions as in (4.248),
the orthogonal relationships are
Application of the spectral method to the sphere can be done by expanding the functional forms Ψ(λ, μ, t) of the different variables as a function of longitude λ [0, 2π] and sine of latitude μ [–1, 1]) using normalized spherical harmonics (see Figure 4.21):
(4.252)
where are the spectral coefficients, which are the unknowns to be determined as a function of time t. The choice of parameters M and N(m) define the truncation of the expansion, and are discussed below. The spherical harmonics are the eigenfunctions of the Laplace operator on the sphere that verify the relation
(4.253)
They are expressed as a combination of sines and cosines (or equivalently by complex exponentials) to represent the periodic variations in the zonal direction, and by real associated Legendre functions of the first kind (see Box 4.11) to account for the variations in the meridional direction. Thus,
(4.254)
Here, index m represents the zonal wavenumber; its highest value M specifies the number of waves retained in the zonal direction. Index n – |m| is called the meridional nodal number.
Figure 4.21 Representation of the characteristics of three spherical harmonics with total wavenumber n = 6. (a): zonal wavenumber m = 0, (b): m = 3 and (c): m = 6.
From Williamson and Laprise (1998).
Box. 4.11 Associated Legendre Polynomials of the First Kind
Associated Legendre polynomials of degree n and order m are real functions of x and solution of the Legendre equation
They can be expressed analytically by
The first polynomials of degree n(0 to 5) for m = 0 are (Box 4.11 Figure 1)
To satisfy the condition of orthonormality, the Legendre functions must obey
and the normalized Legendre function is then expressed by
Box 4.11 Figure 1 First Legendre polynomials of degrees 1 to 5 displayed as a function of x.
Copyright 2010 J. Maddock, P. A. Bristow, H. Holin, X. Zhang, B. Lalande, J. Råde, G. Sewani and T. van den Berg.
The type of truncation to be adopted for expression (4.252) is determined by the relation between the number of waves allowed in the zonal and the meridional directions. If N is chosen to be equal to M, the truncation is said to be triangular. If it is such that N = M + |m|, it is called rhomboidal (Figure 4.22). Triangular truncation is the universal choice for high-resolution models, while rhomboidal truncation is often adopted in the case of low-resolution atmospheric models.
Figure 4.22 Rhomboidal (a) and triangular (b) truncations.
We have seen that chemical concentrations are better described in the physical grid space than in the spectral space. Pseudo-spectral models allow certain processes to be treated on a physical grid while others are treated by the spectral method. In this approach, the governing equations are solved in the physical space by applying, for example, the finite difference method. The spatial derivatives, however, are calculated analytically after converting the physical quantities from the physical space to the spectral space. Rapid forward and inverse transforms of different variables between the physical and spectral space are thus required at each time step, and can be performed very efficiently by the FFT technique (see Appendix E). The aliasing problem arising from nonlinear terms in the transform process is avoided if the number of grid points in the physical grid is equal to 3M + 1 in the zonal direction. The number of points in the meridional direction must be (3M +1)/2 for the triangular truncation and 5N/2 for the rhomboidal truncation. A transformed grid that is often adopted is the Gaussian grid with grid points equally spaced along the longitudes, but not along the latitudes; their location in the meridional direction is defined by the roots μm (m = 1, M) of
(4.255)
The Gaussian grid has no grid point at the pole. In a reduced Gaussian grid, the number of grid points in the zonal direction decreases toward the poles, which keeps the zonal distance between grid points approximately constant across the sphere.
One often wishes to characterize the spatial resolution of a spectral global model in terms of grid spacing L rather than by the highest wavenumber M. The definition of grid spacing that i
s equivalent to a given spectral resolution is not straightforward. Laprise (1992) suggests four possible approaches, whose corresponding grid spacing estimates differ by about a factor of 2. One simple approach is to calculate the average distance between grid points on the Gaussian grid (or equivalently the spacing L1 [km] between longitudinal grid points at the Equator). With the triangular truncation, the equivalent grid spacing is
(4.256)
where a (6378 km) represents the Earth’s radius and M is again the largest wavenumber in the zonal direction. A second measure L2 is half of the shortest resolved zonal wave at the Equator:
(4.257)
A third approach assumes that an equal area of the Earth’s surface is assigned to every piece of information contained in the spherical series, with (M+1)2 real coefficients in the case of a triangular truncation. This yields an area resolution 4πa2/(M + 1)2, corresponding to a length L3
(4.258)
The fourth definition is based on (4.253) that expresses the spherical harmonics as the eigenfunctions of the Laplace operator on the sphere. By equating the eigenvalue of the highest resolved mode with the corresponding eigenvalue of Fourier modes in Cartesian geometry, Laprise (1992) deduces a fourth possible expression L4 for the spatial resolution:
Modeling of Atmospheric Chemistry Page 21