Non-interpolating semi-Lagrangian schemes have also been developed (see, e.g., Ritchie, 1986). In this case, the vector that defines the back-trajectory is decomposed in the sum of a vector that reaches the grid point G closest to the departure point D and a residual vector pointing from this grid point to the departure point D. To determine the value of the transported function at grid point G, no interpolation is needed. The transport along the second vector is performed using a classic Eulerian method, for which the Courant number is always smaller than 1. The overall advection for a given time step is thus always stable, but the advection for the second substep has the dispersive/diffusive properties of the Eulerian scheme that is adopted.
As in the case of the pure Lagrangian methods, SLT algorithms are stable for relatively large time steps. The stability condition is provided by the Lipschitz criterion (trajectories may not cross each other), which is considerably less severe than the CFL condition. To illustrate the performance of the semi-Lagrangian method, Figure 7.27 shows the evolution of a “slotted” cylinder after six revolutions of solid-body rotation at uniform angular velocity about the center of the domain (Staniforth and Côté, 1991). The adopted Courant number of 4.2 is considerably larger than required for the stability of Eulerian schemes. For the same accuracy, the SLT method is more computationally efficient than Eulerian methods.
Figure 7.27 Advection (solid-body rotation) of a “slotted” cylinder using a semi-Lagrangian method with a constant uniform angular velocity about the center of the domain. The adopted Courant number is 4.2. A cubic-spline interpolator is adopted: (a) shows the initial condition and (b) the shape of the cylinder after six revolutions.
From Staniforth and Côté (1991). Copyright © American Meteorological Society, used with permission.
A major disadvantage of SLT grid point schemes is that they do not conserve mass. Numerical adjustment is necessary and different methods can be used for this purpose (Rasch and Williamson, 1990b; see also Section 7.9). The simplest is to compare the total mass Mi(tn+1) of tracer i at time tn+1 over the model domain to the total mass Mi(tn) at time tn, and apply a uniform multiplicative correction Mi(tn)/Mi(tn+1). to the mixing ratios at tn+1.
7.8.2 Finite Volume Based SLT Schemes
To avoid the artificial mass correction process required by SLT grid point schemes, conservative finite volume SLT schemes have been developed (Nair and Machenhauer, 2002; Nair et al. 2003; Zerroukat et al., 2002, 2007; Lauritzen et al., 2006). Rather than transporting the tracer mixing ratio at specific points in the model domain, these schemes advect variable finite volume elements that contain a specified mass of the tracer. By this process, the total mass (or the averaged mass density) of the tracer in the “departure” volume is equal to its mass in the “arrival” volume, ensuring mass conservation. In a 2-D formulation, the conservation of the total mass inside a Lagrangian grid cell of area A that moves and is distorted with the fluid motion is expressed as:
(7.210)
After discretization in time, we have
(7.211)
where and are the surface areas of the cell j at the departure and arrival time levels, respectively, and
(7.212)
is the mean density in grid cell j. The surface at time tn+1 coincides with a Eulerian grid cell of the model. The value of the mean density in the departure area is derived by interpolation of the solution at time tn. The solution at time tn+1 is obtained by a remapping process. The method is illustrated in Figure 7.28 in a 1-D configuration. In this simple case, one defines the mean density of a tracer in a grid cell of size Δx:
(7.213)
The solution at the arrival time tn+1 is
(7.214)
If the wind velocity u varies as a function of space and time, the size of the gridbox (noted Δx at the arrival location) varies during the back-trajectory step and becomes ΔxD ≠ Δx at the departure location.
Figure 7.28 Schematic representation of the 1-D finite volume mass conserving Lagrangian method of Lin and Rood (1996). The “volume” Δx at the arrival time tn+1 is represented by the grid cell AB, and its value ΔxD at departure time tn is represented by the distance A′B′. If the wind velocity u is not constant in space, the displacements at the cell interfaces Dj+1/2 and Dj–1/2 for points A′ and B′ are not equal and Δx is different from ΔxD (Lin and Rood, 1996).
Copyright © American Meteorological Society, used with permission.
Flux-form finite volume SLT schemes are conservative for tracer transport if the winds originate from a general circulation model (GCM) where the exact same scheme was applied to solve the Navier–Stokes equation for momentum. This requirement can be achieved in online GCM simulations of chemical tracers, but is generally not achievable in offline chemical transport models (CTMs) (Jöckel et al., 2001; Horowitz et al., 2003). This problem in offline CTMs often results from inconsistencies between the advection scheme and the surface pressure tendency provided by the dynamical model. A mass fixer is often applied to alleviate this problem (Horowitz et al., 2003), but may violate monotonicity requirements and generate non-physical transport.
7.9 Spectral, Finite Element, and Spectral Element Methods
The spectral method is widely used to solve the dynamical equations in GCMs. In this approach, the fields such as the temperature or the winds are represented by a series of continuous basis functions such as spherical harmonics in the horizontal plane (see Chapter 4) and by a finite difference formulation in the vertical direction. Spectral methods have also been used in these models to represent the advection of moisture. To illustrate the methodology, we consider again the 1-D case in the x-direction and we approximate the solution of the advection equation by the expansion
(7.215)
in which Φk(x) represents a set of orthogonal functions (e.g., elementary trigonometric functions) and ak are unknown coefficients that depend on time t. Expression (7.215) is introduced in the advection equation (7.17), leading to a system of K differential equations for the coefficients ak(t) that can be solved by standard methods. When spherical geometry is used, the solution is instead expressed as a function of spherical harmonics (see equation (4.252)).
With the proper choice of parameters, the spectral method provides accurate and stable results and can be mass-conservative. However, it is not shape-preserving (monotonic, positive definite); overshoots and undershoots can be produced (see Box 4.7). Negative values of the transported fields can be eliminated by using appropriate filtering and filling schemes, but these corrupt the correlations between tracers which is critical for nonlinear chemistry. In addition, the algorithm does not satisfy the criterion of locality. Because of these limitations, spectral methods are generally not well suited for chemical applications.
In the finite element method, the solution ψ(x, t) is also provided by expansion (7.215), but with the basis functions Φk(x) defined only on a small region of space [A, B] called finite elements (see Section 4.10). In the Galerkin approach, the coefficients ak are derived by requiring that the error arising from representing function ψ(x, t) by expansion (7.215),
(7.216)
be orthogonal to the basis functions. This condition is expressed by the integral over the domain [A, B]:
(7.217)
for all values of i ∈ [0, K]. If in (7.217) one replaces the error eK by expression (7.216), one obtains the system of K-coupled ODEs:
(7.218)
which is solved to obtain the coefficients ak(t). For this purpose, the time derivatives are usually approximated by finite differences.
The spectral element method (Patera, 1984) is a finite element technique in which a high-degree spectral method is applied within each element. As discussed by Nair et al. (2011), the spectral element method combines the geometric flexibility of the traditional finite element methods with the high accuracy, rapid convergence, and weak numerical dispersion and dissipation of the classical spectral methods. It is not inherently conservative, but can be engineered to ensure a u
ser-required level of mass conservation.
Although finite element and spectral element methods have so far mainly been used for engineering applications, they are now emerging as promising methods for atmospheric problems. Their local domain decomposition property makes them particularly well suited for massively parallel computer architectures, and they can be easily applied when the model domain is geometrically complex or when grid refinement is needed in specified atmospheric regions.
7.10 Numerical Fixers and Filters
An important requirement for the numerical method applied to solve the advection equation is that the solution be positive and that mass be conserved. This is not always the case, and if so one can apply a-posteriori corrections to restore these desired properties. Jablonowski and Williamson (2011) provide an extensive review of the use of fixers and filters in atmospheric models.
Fixers. Negative values in the transported variables (e.g., tracer mixing ratio) can be eliminated by introducing “numerical fixers” that borrow mass from surrounding (or downstream) grid points. Rasch and Williamson (1990a) describe a possible implementation of local and global fixers. When a negative value is encountered during a point-by-point scan of the calculated quantities on each horizontal surface, the immediate neighboring points are examined, and if sufficient mass is available to “fill the hole,” the negative value is set to zero and the values at the neighboring points are reduced proportionally. If there is not sufficient mass available, no action is taken. The application of this local filling method does conserve mass, but may not eliminate all negative values. In a second step, a global filter is applied in which the remaining negative values are set to zero; this violates mass conservation but a renormalization can be applied by scaling to the global masses in the domain, as described in Section 7.8.1 to enforce mass conservation in SLT schemes. This filling process produces diffusion and does not ensure monotonicity of the corrected fields; it can also be computationally expensive.
Filters. In several of the algorithms described previously, short waves may grow excessively in the solution of the advection equation, producing undesired noise and even catastrophic instability. These waves can be eliminated by appropriate smoothing or filtering. One option is to add a small diffusive term in the advection equation, which will smooth the solution and suppress high wavenumbers. In spectral models, numerical noise can easily be suppressed by omitting wavenumbers larger than a specified value and highlighting only the scales of interest. Spectral filtering can also be used in grid point models by applying a Fourier transform to the solution, which eliminates the high frequencies in the signal, and applying an inverse Fourier transform. Such filtering is often applied in polar regions where the meridians converge and a longitudinal grid point spacing becomes so small that, without filtering, the solution would become unstable. Finally, numerical filters such as the linear Shapiro filter presented in Section 4.15.2 are often applied to eliminate two-grid interval waves completely while having little effect on longer waves (Shapiro, 1971). Other high-frequency filters have been developed by Asselin (1972) and Forester (1977).
7.11 Concluding Remarks
In this chapter we have reviewed different numerical algorithms used to approximate the solution of the linear advection equation. No existing method fully addresses modelers’ requirements. The examination of simple numerical schemes reveals that high-order algorithms are generally not monotonic and occasionally produce undesired negative values. Low-order algorithms such as the upstream method preserve the sign of the solution, but are excessively diffusive. Thus, practical applications must adopt more elaborate schemes that address some of the drawbacks that characterize the simple methods. Modern schemes are often upstream-based Eulerian finite volume methods that are mass conservative, positive definite, and possess good phase-error characteristics. They may use adaptive time steps to meet CFL stability requirements or semi-Lagrangian options to get around these requirements. Specific, often complex nonlinear algorithms are developed to reduce the numerical diffusion that characterizes upstream methods. These complex schemes can yield significant improvement in accuracy, but often with enhanced computational costs.
Finite volume Eulerian methods, in which a subgrid distribution of the transported quantity is specified, provide highly accurate solutions and are popular in global CTMs. In many respects, they are superior to classic grid point methods. Computational cost depends on the user tolerance for numerical diffusion. The Prather scheme is regarded as a reference among Eulerian models. It produces accurate solutions with little diffusion. However, it has large computing time and storage requirements. A van Leer scheme may enable higher grid resolution, compensating for the lower accuracy.
Lagrangian methods are popular for source-oriented and receptor-oriented transport problems in which one is concerned with transport from a point source or transport contributing to concentrations at a receptor point. However, they do not provide the regular full-domain solution achievable by Eulerian methods and cannot properly represent nonlinear chemistry or aerosol microphysics. Semi-Lagrangian methods are very popular in global CTMs because their numerical stability is not as severely constrained by choice of time step as in the case of Eulerian schemes. They are sometimes used as a back-up scheme in cases where the regular Eulerian solver violates the CFL criterion.
In summary, there is no single advection scheme that is universally best. The choice of scheme depends on the type of problem being solved, the tolerance for different kinds of error, and the computational demands. For any scheme, it is important to verify that basic criteria of stability and mass conservation are met. The material in this chapter should enable readers to understand the issues associated with different advection schemes and to make informed choices in selecting appropriate schemes for their applications.
References
Allen D. J., Douglass A. R., Rood R. B., and Guthrie P. D. (1991) Application of a monotonic upstream-biased transport scheme to three-dimensional constituent transport calculations, Mon. Wea. Rev., 119, 2456–2464.
Asselin R. (1972) Frequency filter for time integrations, Mon. Wea. Rev., 100, 487–490.
Bermejo R. (1990) Notes and correspondence on the equivalence of semi-Lagrangian schemes and particle-in-cell finite element methods, Mon. Wea. Rev., 118, 979–987.
Boris J. P. and Book D. L. (1973) Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works, J. Comput. Phys., 11, 38–69.
Boris J. P. and Book D. L. (1976) Flux-corrected transport. III. Minimal-error FCT algorithms, J. Comput. Phys., 20, 397–431.
Bott A. (1989a) A positive definite advection scheme obtained by non-linear renormalization of the advective fluxes, Mon. Wea. Rev., 117, 1006–1015.
Bott A. (1989b) Reply, Mon. Wea. Rev., 117, 2633–2636.
Carpenter R. L., Droegemeier K. K., Woodward P. R., and Hane C. E. (1990) Application of the piecewise parabolic method (PPM) to meteorological modeling. Mon. Wea. Rev., 118, 586–612.
Chen Y. and Falconer R. A. (1992) Advection–diffusion modelling using the modified QUICK scheme, Int. J. Numer. Meth. Fluids, 15, 1171–1196.
Chlond A. (1994) Locally modified version of Bott’s advection scheme, Mon. Wea. Rev., 122, 111–125.
Colella P. and Sekora M. D. (2008) A limiter for PPM that preserves accuracy at smooth extrema, J. Comput. Phys., 227 (15), 7069–7076.
Colella P. and Woodward P. R. (1984) The piecewise parabolic method (PPM) for gasdynamical simulations, J. Comput. Phys., 54, 174–201.
Courant R., Friedrichs K., and Lewy H. (1928) Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann., 100, 32–74.
Courant R., Isaacson E., and Rees M. (1952) On the solution of nonlinear hyperbolic differential equations by finite difference, Commun. Pure Appl. Math., 5, 243–255.
Crowley W. P. (1968) Numerical advection experiments, Mon. Wea. Rev., 96, 1–11.
Davies H. C. (1983) Limitations of some common lateral boundary
schemes used in regional NWP models, Mon. Wea. Rev., 111, 1002–1012.
Dukowicz J. K. and Baumgardner J. R. (2000) Incremental remapping as a transport/advection algorithm, J. Comput. Phys., 160, 318–335.
Dullemond C. P. (2009) Numerical Fluid Dynamics, Lecture Notes, Zentrum für Astronomie, Ruprecht-Karls Universität, Heidelberg.
Farrow D. E. and Stevens D. P. (1994) A new tracer advection scheme for Bryan Cox type ocean general circulation models, J. Phys. Oceanogr., 25, 1731–1741.
Forester C. K. (1977) Higher order monotonic convective difference schemes, J. Comput. Phys., 23, 1–22.
Godunov S. K. (1959) A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics, Mat. Sb., 47, 357–393.
Gross E. S., Koseff J. R., and Monismith S. G. (1999) Evaluation of advective schemes for estuarine salinity simulations, J. Hydr. Engrg., ASCE, 125 (1) 32–46.
Horowitz L. W., Walters S., Mauzerall D. L., et al. (2003) A global simulation of tropospheric ozone and related tracers: Description and evaluation of MOZART, version 2, J. Geophys. Res., 108 (D24), 4784, doi: 10.1029/2002JD002853.
Iserles A. (1986) Generalised leapfrog methods, IMA J. Numer. Anal., 6, 381–392.
Jablonowski C. and Williamson D. L. (2011) The pros and cons of diffusion, filters and fixers in atmospheric general circulation models, In Numerical Techniques for Global Atmospheric Models, Lecture Notes in Computational Science and Engineering (Lauritzen P. H., Jablonowski C., Taylor M. A., and Nair R. D., eds.), Springer, Berlin.
Modeling of Atmospheric Chemistry Page 41