where the trigonometric functions have been approximated using Taylor expansions. The resulting lag per unit time in the deplacement of the waves increases with the square of their wavenumber k since
(7.57)
The propagation of the waves is therefore fastest for the shortest wavelengths (and thus for wavelengths that approach the grid size Δx), which generates ripples in the advected signal. This type of behavior, shown here in the case of the Lax–Wendroff algorithm, is common to all second-order schemes and, as stated by Godunov’s theorem, the monotone behavior of a numerical solution cannot be assured for linear finite difference methods with more than first-order accuracy. This theorem introduces a major limitation in the development of numerical schemes that treat advection.
Implicit Schemes
The implicit or Euler backward-in-time, centered-in-space (BTCS)
(7.58)
with the recursive expression
(7.59)
is first-order accurate in time and second-order in space. The amplification coefficient derived from the von Neumann analysis is:
(7.60)
The resulting amplification factor
(7.61)
is smaller than unity for any value of the Courant number. The method is therefore unconditionally stable, allowing for the adoption of any arbitrary time step Δt, which is a great advantage. The solution, however, cannot be retrieved as easily as in the case of explicit schemes. In the implicit case, a system of J algebraic equations (if J is the number of grid points that are not associated with boundary conditions) must be solved, which is computationally expensive. In the 1-D case, the system of equations is tridiagonal and can be solved efficiently with the Thomas algorithm (see Box 4.4). Another limitation of the method is that it has limited accuracy, with the shortest wavelengths being more rapidly attenuated that the longer wavelengths.
Accuracy of the solution can be improved by combining the FTCS and BTCS approaches. The Crank–Nicholson algorithm, written as
(7.62)
is second-order accurate in time and space. It is implicit because it includes terms evaluated at time tn+1 on the right-hand side. The amplification coefficient is given by
(7.63)
and the amplification factor is equal to 1 for all wavenumbers and all Courant numbers:
|g(k)| = 1
(7.64)
The algorithm is thus unconditionally stable.
Matsuno Scheme
The Matsuno scheme is a two-step explicit–implicit algorithm that is first-order accurate in time and second-order accurate in space. The first step is to predict an intermediate value of the transported quantity at time level n + 1 by using a simple FTCS (Euler forward) approach:
(7.65)
The predicted values are then substituted into the space derivative term, and a correction step is applied:
(7.66)
By eliminating the intermediate terms Ψ*, one derives after some manipulations:
(7.67)
This explicit expression approximates an advection equation with an additional diffusion term that approaches zero for very small Δt:
(7.68)
The amplification coefficient is:
g(k) = 1 − i α sin (kΔx) − α2sin2(kΔx)
(7.69)
with corresponding amplification factor:
(7.70)
Even though this scheme bears some resemblance to implicit schemes, it is stable for the usual Courant condition (α ≤ 1) of the explicit method rather than the condition that applies to implicit methods (unconditional stability).
If, in the Matsuno scheme, the correction step is replaced by
(7.71)
in which the space derivative term is calculated as the average between the intermediate estimates and the estimates at time tn, we obtain the Heun scheme (see Table 7.1), which is second-order accurate in space and time like the Crank–Nicholson algorithm described earlier. The method is unconditionally unstable unless a small diffusion term is artificially added to the advection equation. In this case, the scheme becomes conditionally stable. If, in the Heun scheme, the predictor step is a leapfrog algorithm (see Section 7.3.3), we obtain the method proposed by Kurihara (1965), which is second-order accurate in space and time, stable for the Courant condition, and free of numerical diffusion. Unlike the leapfrog method, it is not subject to drift, but it does not provide the exact solution for α = 1.
Table 7.1 Elementary algorithms for solving the 1-D advection equation
Method Algorithm Stability Accuracy Remarks
Euler Forward Unconditionally unstable Δt, Δx2
Lax Stable for α < 1 Δt, Δx2 Diffusive
Leapfrog Stable for α < 1 Δt2, Δx2 Dispersive
Lax–Wendroff Stable for α < 1 Δt2, Δx2
Implicit Unconditionally stable Δt, Δx2
Crank–Nicholson Unconditionally stable Δt2, Δx2
Matsuno Stable for α < 1 Δt, Δx2 Diffusive
Heun Unconditionally unstable Δt2, Δx2
Kurihara Stable for α < 1 Δt2, Δx2 Not diffusive
Fourth-order (implicit) Unconditionally stable Δt, Δx4
Upstream (α > 0) Stable for α < 1 Δt, Δx Monotonic diffusive
Upstream (α < 0) Stable for α < 1 Δt, Δx Monotonic
diffusive
Figures 7.7 and 7.8 show the amplification factors |g| as a function of parameters k Δx and c Δt/Δx for several of the algorithms described above.
Figure 7.7 Amplification factor for different numerical methods as a function of parameter k Δx for a Courant number α = 0.5.
Figure 7.8 Amplification factor for different numerical methods as a function of the Courant number α for parameter kΔx = π/2 (corresponding to wavelength L = 4Δx).
Fourth-Order in Space Method
The algorithms discussed so far use low-order explicit or implicit forms of the finite difference equations. These algorithms can be extended to higher-order formulations. For example, the forward-in-time implicit form of the fourth-order approximation
(7.72)
is unconditionally stable. The matrix corresponding to this system is a banded matrix with five terms that can be inverted with a fast method (Press et al., 2007).
7.3.2 Methods Using Space-Uncentered Differences
In the algorithms discussed previously, the space derivative ∂Ψ/∂x is approximated by a second-order accurate centered difference. An alternative approach is to adopt a first-order accurate backward-in-space finite difference,
(7.73)
When introduced in the 1-D advection equation (7.17) together with a forward in time derivative, one obtains the upstream (or upwind differencing) method (Courant et al., 1952; Godunov, 1959). Consistent with physical considerations, this algorithm provides a solution that depends on the behavior of Ψ in the direction from which the flow emanates, and not from the function downstream. Thus, for c > 0, we write a forward-in-time, backward-in-space (FTBS) expression:
(7.74)
or equivalently:
(7.75)
For c < 0, the advection equation is approximated by a forward-in-time, forward-in-space (FTFS) expression:
(7.76)
or
(7.77)
The amplification coefficient (for c > 0) is:
g(k) = 1 − α[1 − cos (kΔx)] − iα sin (kΔx)
(7.78)
with amplitude
(7.79)
The amplitude remains below unity as long as the Courant condition (α ≤ 1) is verified. The phase Φ(k) is given by:
(7.80)
At the stability limit, when α = 1, the amplitude |g(k)| = 1 and the phase Φ = −k Δx. In this case, the solution provided by the upstream scheme is exact. In the general case with α < 1, the solution is dampened (numerical diffusion) with the highest wavenumbers (or smallest wavelengths) more rapidly attenuated than the lower wavenumbers. This explains why the sharp corners of the square waves in Figures 7.10 and 7.12 are rounded by the upstream metho
d.
The upstream method is monotonic and sign-preserving, but it is only first-order accurate in space and time and suffers therefore from numerical diffusion. This point can be intuitively understood by noting that the algorithm expressions (7.74) and (7.76) approximate to second-order in Δx and Δt the advection–diffusion equation:
(7.81)
with diffusion coefficient K = 0.5(c Δx − c2Δt) = 0.5 c Δx(1 − α).
Uncentered methods other than the upstream scheme have been proposed to reduce excessive numerical diffusion. For example, the approximation proposed by Warming and Beam (1976)
(7.82)
is second-order accurate in time and space and is stable for 0 ≤ α ≤ 2. It is equivalent to a Lax–Wendroff scheme in which the centered space differences are replaced by backward differences.
The Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme of Leonard (1979) employs four points to approximate the first-order space derivative. For a constant wind velocity c ≥ 0 and grid spacing Δx, the advection equation (7.17) is first discretized as a centered-in-space, time-forward explicit scheme:
(7.83)
where the values of the advected quantity at the left (j – 1/2) and right (j + 1/2) edges of cell j are determined by a quadratic interpolation. One derives, for example:
(7.84)
so that
(7.85)
The scheme is second-order accurate in space, but it is unstable unless some dissipation is added to the advection equation. Other formulations of the QUICK scheme (i.e., explicit, implicit, or semi-implicit approaches) are available (Chen and Falconer, 1992). A more elaborate algorithm, called QUICKEST, also proposed by Leonard (1979), is third-order accurate in time and space, and is stable for pure advection if α ≤ 1. In that scheme, the value of the function at the right edge is given by:
(7.86)
The QUICK and QUICKEST schemes often generate overshoots and undershoots. They can therefore produce negative solutions. This problem is addressable by imposing flux-limiters in the integration scheme (see Section 7.5).
Finally, the algorithm proposed by Farrow and Stevens (1994), which can be regarded as an adaptation of the QUICK scheme, is expressed as a predictor-corrector integration scheme
(7.87)
(7.88)
It is third-order accurate in space and second-order in time. A von Neumann stability analysis indicates that it is stable for Courant numbers smaller than approximately 0.6.
7.3.3 Multilevel Algorithms
In the numerical schemes discussed in previous sections, the time derivatives are approximated by a two-level forward difference. We now consider methods in which information from several earlier time levels are used to calculate the value of function Ψ at time tn+1.
The regular leapfrog method (Courant et al., 1928), which is second-order accurate in time, is based on a centered-in-time and centered-in-space (CTCS) approximation of the advection equation:
(7.89)
or
(7.90)
In this three-level algorithm, the solution “leapfrogs” from time level (n – 1) to time level (n + 1) over the time level (n) at which the space derivative term is computed.
The von Neumann stability analysis provides a quadratic equation for the amplification coefficient, whose two roots are:
g(k) = ± [1 − α2sin2(kΔx)]1/2 + i α sin (kΔx)
(7.91)
If |α sin (k Δx)| > 1, the square root term is completely imaginary, and the modulus |g(k)| for one of the two roots is larger than 1, indicating instability. If |α sin (k Δx) ≤ 1|, which is verified for all wavenumbers when |α| ≤ 1 (CFL condition), the modulus is unity:
|g(k)| = {[1 − α2sin2(kΔx)] + [α sin (kΔx)]2}1/2 = 1
(7.92)
and the phase shifts for the two roots (±) are respectively
Φ+ = − sin−1(α sin(kΔx)) and Φ− = π + sin−1(α sin(kΔx))
(7.93)
For Courant stable conditions, the amplitude of all waves is preserved, not dissipated. This represents the major advantage of the method. When α = 1, the method provides the exact solution (correct amplitude and phase). If α < 1, computational dispersion occurs as phase errors, particularly for short waves, and leads to some spurious numerical oscillations.
Expression (7.91) with a ± sign shows that the leapfrog algorithm generates two solutions with different amplification functions. One of them, called the physical mode, represents the meaningful solution, while the second one, referred to as the computational mode, is a mathematical artifact without any physical reality. This second solution propagates in the direction opposite to the flow and changes sign for every time step; it generates therefore noise that needs to be filtered with an appropriate method (see Section 4.15.4). The effect of the computational mode is visible in Figure 7.9, which shows the advection of a cosine-shaped function. Undesired oscillations with negative values of the function are produced upwind from large spatial gradients. The computational mode is most strongly excited when the initial conditions are characterized by sharp gradients.
Figure 7.9 Advection by the leapfrog scheme of a cosine-shaped function with a half-width resolution of 12Δx. The uniform grid is composed of 500 cells. The Courant number adopted in this example is 0.5. The solution is shown after 1600 time steps Δt.
From Smolarkiewicz (2006).
The leapfrog algorithm tends to decouple odd and even grid points. Although, in principle, the solutions at these two types of grid point should not diverge, in practice they often do so as time progresses, causing checkerboarding of the solution. The problem can be addressed by adding a small dissipative term, by discarding the solutions at one of the two types of grid points, or by switching occasionally to an alternate advection scheme for just one time step.
Different improved leapfrog schemes have been proposed (Kim, 2003). The upwind leapfrog scheme introduced by Iserles (1986),
(7.94)
or
(7.95)
is characterized by a considerably lower phase error than the regular leapfrog scheme. Accuracy can be increased by adopting a fourth-order accurate spatial discretization:
(7.96)
for which the von Neumann stability analysis leads to:
(7.97)
One derives easily that the scheme is stable if
(7.98)
For α < 0.73, the scheme is stable for all harmonics. The use of higher orders for the calculation of the time derivative also improves the accuracy of the solution. For example, the four-level algorithm:
(7.99)
is very accurate and leads to exact solutions when α = 1/2 or α = 1/3. However, it is unstable for α > 1/2.
Higher-order multi-stage methods are more accurate, but have the disadvantage of generating a larger number of computational modes. An interesting case is the third-order Adams–Bashforth scheme (see 4.197)):
(7.100)
because the two undesired computational modes that are produced in this case are strongly damped if |α| < 0.72. No filtering is required in most applications and this makes the algorithm particularly attractive, even though the solution is not positive definite. For higher values of |α|, the amplitude of one of the computational modes becomes larger than 1, and the scheme becomes unstable.
7.3.4 Performance of Elementary Finite Difference Algorithms
Table 7.1 summarizes the properties of the different algorithms presented previously. Figure 7.10 shows a comparison for the 1-D advection of an initial square function with constant wind speed. The instability of the Euler forward algorithm is manifested in the large oscillations. The upstream algorithm preserves sign and is free of oscillations (negligible phase lag), but it is very diffusive. The Lax and Matsuno algorithms are also very diffusive. The leapfrog algorithm is not diffusive and conserves the concentration variance but it produces oscillations with undesirable negative values. The Lax–Wendroff method is slightly diffusive and produces small unwanted oscillations
with negative values. Filters are generally applied to avoid unphysical negative values, but such filters may destroy the conservation properties of the numerical algorithms. The stencils associated with some of the algorithms are shown in Figure 7.11. Finally, Figure 7.12 shows the performance of different algorithms in the case of the diagonal advection of a square wave in a 2-D domain. Again, one notes the strong numerical diffusion associated with the upwind scheme and the presence of large oscillations in the case of the leapfrog scheme. The multidimensional definite advection transport algorithm (MPDATA) (Smolarkiewicz, 1984; see Section 7.6) provides positive definite solutions, but with large overshoots. The Lax–Wendroff scheme with flux limiters (see Section 7.5) performs best and the QUICKEST algorithms exhibit significant oscillations (Gross et al., 1999).
Figure 7.10 Comparison between exact analytic (red) and numerical (black) solutions of the 1-D advection equation for a square function. The velocity c is constant. The different numerical algorithms are labeled in the panels. The original square function is centered at x = 20 and the adopted Courant number is equal to 0.5. The periodic boundary condition for the advected field is zero at x = 0 and x = 100. The results are shown after 40 time steps. The Euler forward algorithm is unstable (note the different scale used for the y-axis).
Modeling of Atmospheric Chemistry Page 37