Figure 7.4 Solution of the 1-D advection–diffusion equation for a Péclet number of the order of 1. The initial function at time t = t0 is represented by a boxcar function. The gradual displacement of the function is due to the advection term and deformation (spread) of the shape of the function is caused by the diffusion term.
Adapted from Slingerland and Kump (2011).
When considering a discretized form of the advection–diffusion equation, one often introduces the numerical Péclet number as:
(7.16)
where Δx represents the grid size of the model. The numerical Péclet number measures the relative importance of advection and turbulent diffusion at the smallest spatial scale resolved by the model. The advection–diffusion equation is either parabolic (diffusion-dominated) or hyperbolic (advection-dominated), depending on the Péclet number. Numerical algorithms treating atmospheric diffusion are discussed in Chapter 8.
7.3 Elementary Finite Difference Methods
We examine here the fundamental properties of several advection algorithms applied to the simple 1-D advection equation (along direction x), with a constant velocity c (taken to be positive) and fixed grid spacing Δx. Generalization to variable wind velocity and grid spacing is presented in Section 7.3.5. As in the previous chapters dealing with numerical algorithms, we represent the field of the transported quantity by the generic mathematical symbol Ψ. The advection of a non-negative scalar function Ψ is described by the first-order hyperbolic PDE:
(7.17)
To solve this equation, initial and boundary conditions must be specified. The initial condition can be expressed as
Ψ(x, 0) = G(x)
(7.18)
where G(x) represents the spatial distribution of the tracer distribution at t = 0. The resulting analytic solution is simply the translation without any change in shape of function G(x) in the x-direction at a velocity c. Thus,
Ψ(x, t) = G(x − ct)
(7.19)
When the spatial domain, rather than being infinite, extends from x = a to x = b, a condition must be specified at one boundary of the domain. An advection problem (hyperbolic equation) is well-posed if a boundary condition on the value of Ψ (Dirichlet condition) is imposed at the inflow boundary. If c > 0, the condition must be expressed at x = a:
Ψ(a, t) = Ha(t)
(7.20)
while if c < 0 the condition must be applied at x = b:
Ψ(b, t) = Hb(t)
(7.21)
A periodic boundary condition, such as on a sphere, can be imposed as
Ψ(b, t) = Ψ(a, t)
(7.22)
In this case, the mass that flows out of the domain at boundary x = b flows back into the domain at boundary x = a.
Function Ψ(x, t) can be represented in an unbounded or periodic domain as a Fourier series with components (harmonics) characterized by their wavenumbers k [m–1] corresponding to wavelengths L = 2π/k:
(7.23)
In order to analyze the fundamental properties of different elementary algorithms, we assume that the problem is periodic in space and consider that a single harmonic (k) of function Ψ(x, t) at grid point xj = j Δx and time t takes the value
Ψk(xj, t) = Ak(t)eik jΔx
(7.24)
Here, Δx represents the grid spacing, assumed to be uniform, and j is the grid index. The advection equation becomes:
(7.25)
In this idealized case, the advection velocity c can be interpreted as the phase speed of the wave defined by Ψk. Because c is constant, the phase speed is independent of wavenumber k. Therefore all waves propagate at the same speed. Accurate numerical algorithms have to preserve this property as much as possible. In this particular case, the group speed cg, which is indicative of the rate at which the energy propagates, is equal to c and is independent of the wavenumber k. This analytic form of the advection equation will later be compared to the forms provided by numerical analogs.
7.3.1 Methods Using Centered Space Differences
When adopting a uniform grid spacing Δx, the space derivative ∂Ψ/∂x can be approximated by a centered space difference to yield the second-order accurate expression:
(7.26)
where O(Δx2) is the truncation error. For a wave with wavenumber k, the finite difference expression can be written as
(7.27)
Under these conditions, the approximate form of the advection equation
(7.28)
can be compared with the exact analytic equation (7.25). When centered finite differences are used to approximate space derivatives, the phase velocity c* associated with the numerical solution
(7.29)
varies with wavenumber k, while the true phase velocity c is independent of k (see Figure 7.5). Thus, even though all wavenumber components that characterize function Ψ should move at exactly the same speed, the shorter wavelengths are trailing the longer waves. As a result, the different Fourier components of the advected function are displaced along axis x at different velocities and the numerical solution is distorted. This property that arises from the space differencing is named numerical dispersion and leads to phase errors. In space-centered approximations, the shortest wavelength that can be resolved by the model (L = 2Δx or kΔx = π) does not move at all since its phase speed is zero (see Figure 7.5). For long waves (small values of kΔx), the phase speed c* provided by the numerical scheme approaches the true value of c. Figure 7.5 also shows the variation of the group velocity ,
(7.30)
for the centered difference scheme. For wavelengths between 4Δx and 2Δx, the group velocity is negative. This means that energy can propagate upstream, which is an undesirable property of the scheme.
Figure 7.5 Ratios between the phase velocity c* associated with the numerical solution and the true phase velocity c as a function of kΔx for the second-order and fourth-order space derivatives. Here k refers to the wavenumber of the different Fourier components of the signal and Δx to the grid spacing. The corresponding wavelengths L for three particular waves are indicated. The graph highlights the lag in the advection of the waves relative to the advective motion. This effect is more pronounced for the shortest wavelengths. The ratio between the group velocity cg resulting from the algorithm and the true group velocity c is also shown in the case of the second-order derivative.
When the space derivative is approximated by a fourth-order scheme over a grid with constant spacing,
(7.31)
the phase velocity c* becomes:
(7.32)
This improves over the second-order scheme, as shown in Figure 7.5, but still fails for wavelengths close to L = 2Δx.
Euler Forward Scheme (FTCS)
A simple numerical method to solve (7.17) is the Euler forward scheme, which approximates the time derivative by a forward difference and the space derivative by a centered difference:
(7.33)
where Δt is the time step and Δx the grid spacing (both assumed to be constant), and n and j are the indices referring to time and space, respectively, with
tn = n Δt and xj = j Δx
This algorithm, which is also referred to as the FTCS method (forward-in-time, centered-in-space), is first-order accurate in time and second-order accurate in space. Its solution,
(7.34)
where
(7.35)
is the Courant number, can be easily computed because the algorithm is explicit (the solution at time tn+1 for each point xj is derived directly from quantities that are already known at time tn). The method is also one-step because only one calculation is required to advance the integration from time level tn to the new time level tn+1. It is a two-level scheme since only two time levels (tn and tn+1) are involved in the calculation.
An important consideration when assessing an algorithm is its stability. To address this, we apply the von Neumann’s stability analysis presented in Box 7.1. For the FTCS algorithm, the amplification coefficient g(k) as a function of wavenumber k is:
&nbs
p; g(k) = 1 − iα sin (kΔx)
(7.36)
The amplification factor,
(7.37)
is greater than one for all values of the wavenumber k. As a result, any numerical error produced by the algorithm grows exponentially with time. As highlighted in Box 7.1, the Euler FTCS algorithm is thus unconditionally unstable and must be rejected.
Box 7.1 The von Neumann Stability Analysis
The von Neumann analysis provides a methodology for assessing the stability of numerical algorithms. It applies to linear PDEs with periodic boundary conditions. We first consider the analytic solution of the advection equation by noting that any function Ψ can be expressed as the superposition of an infinite number of waves. We perform therefore a discrete Fourier transform of function Ψ, which is advected in the x-direction with a positive and constant velocity c. Consider the advection of a single wave harmonic with wavenumber k. At time t, this harmonic is expressed by:
Ψk(x, t) = a eik(x − ct)
After a time interval Δt corresponding to a displacement of the wave over a distance cΔt, function Ψ takes the form:
Ψk(x, t + Δt) = a eik[x − c(t + Δt)] = Ψk(x, t)eikcΔt
The amplification coefficient or gain g(k) for harmonic k is the complex function defined as the ratio between function Ψk after and before the advection step. Thus,
with a modulus (also called amplitude and here amplification factor)
|g(k)| = g(k) ⋅ g*(k) = 1
and a phase
φ(k) = kcΔt
Here g* is the complex conjugate of g. The amplification factor |g(k)| represents the relative change in the amplitude of the harmonic of wavenumber k after one computational time step.
These relations highlight two properties that should be reproduced as closely as possible when numerical approximations to the exact solution are sought: (1) the amplitude of all wave harmonics is unchanged during an advection process, and (2) the phase of harmonics k varies according to φ = α kΔx, where α is the Courant number and Δx the grid spacing of the model. If, for example, the value of g(k) resulting from the use of a numerical approximation is less than 1, the amplitude of the wave is reduced by advection. Conversely, if it is larger than 1, the amplitude of the wave grows and the method becomes rapidly unstable. Thus, the numerical stability of an algorithm for advection requires that |g(k)| ≤ 1 for all waves resolved by the model. Similarly, errors on the phase cause waves of a spectrum to lag the displacement of other waves, leading to numerical dispersion.
To apply the von Neumann analysis to a numerical algorithm, let us assume that the advection equation (7.17) is approximated by the FTCS algorithm (centered difference scheme for the space derivative):
where α is the Courant number. We apply a discrete Fourier transform and consider a single harmonic with wavenumber k:
It results from the above FTCS formulation that:
The transfer function based on the FTCS algorithm is therefore:
Its amplitude is given by
and the phase Φ(k) associated with this particular algorithm is derived from:
A comparison of these expressions with the values derived from the analytic solution shows that the FTCS scheme is unconditionally unstable since the modulus of the transfer function is larger than one for all values of wavenumber k, even for very small time steps. The phase error E(k) is given by:
E(k) = Φ(k) − φ(k) = Φ(k) − α kΔx
The von Neumann analysis is adopted in this chapter to assess the stability conditions of different numerical schemes. It is a necessary and sufficient condition for stability in the case of linear finite difference equations with constant coefficients. It is a necessary but not sufficient stability condition for nonlinear equations.
However, the FTCS method can become conditionally stable if a numerical diffusion term is added to the right-hand side of the advective equation, with diffusion coefficient K:
(7.38)
In this case, the discretized equation becomes:
(7.39)
Here α is again the Courant number, and
(7.40)
is the so-called Fourier number. The amplification coefficient for wavenumber k derived from the von Neumann’s analysis becomes
g(k) = 1 + 2β[cos(kΔx) − 1] − iα sin (kΔx)
(7.41)
with the corresponding amplitude (Figure 7.6)
(7.42)
and phase
(7.43)
Necessary and sufficient conditions for the stability of this scheme are
(7.44)
Thus, the FTCS method can be used if a small diffusion term is added to the advection equation and the time step is sufficiently short.
Figure 7.6 Amplification factor |g(k)| as a function of kΔx when the advection equation is approximated by the FTCS algorithm for a constant velocity c and the Courant number equal to 0.5. Parameter β = kΔt/Δx2 represents the effect of added diffusion to the advection equation. The case with β = 0 corresponds to pure advection and is unconditionally unstable (amplification factor > 1). When diffusion is added, the solution is stable (|g(k)| < 1) for β < 0.5, but becomes unstable for larger values of β.
Lax Scheme
In the Lax method (Lax, 1954), the term used in the FTCS scheme is replaced by ½ , and the approximate form of the advection equation becomes:
(7.45)
When applying the von Neumann stability analysis, we obtain the amplification coefficient:
g(k) = cos (kΔx) − i α sin (kΔx)
(7.46)
whose amplitude
(7.47)
remains less than or equal to unity for Courant numbers α ≤ 1.
This stability condition α ≤ 1, which applies to many Eulerian schemes, is called the Courant–Friedrichs–Lewy (CFL) criterion. Over a time step Δt, the displacement of tracer should never exceed a distance larger than the grid spacing Δx. When a longitude–latitude grid is adopted in the model, this condition imposes severe limitations in the vicinity of the pole where the grid spacing in the longitudinal direction becomes very small. This is often circumvented by the application of numerical filters (see Section 7.10) or by the use of a reduced grid (see Section 4.7.3). The CFL condition also requires that increases in the spatial resolution of a model be accompanied by a proportional decrease in the value of the time step.
The stabilization of the solution in the Lax scheme can be understood by rearranging (7.45) as
(7.48)
which is the FTCS form of equation
(7.49)
Thus, the Lax scheme stabilizes the FTCS solution by adding a numerical diffusion term to the advection equation with diffusion coefficient Δx2/2Δt. Because of this numerical diffusion, a disturbance at grid point j propagates not only to downwind grid point (j + 1), but also to upwind point (j – 1). There results a transportivity error.
Lax–Wendroff Scheme
The Lax–Wendroff scheme (Lax and Wendroff, 1960, 1964), also called the Leith (1965) or Crowley (1968) scheme, is designed to provide the minimum amount of added numerical diffusion required to provide stability to the FTCS solution. The value of advected function Ψ at time level n + 1 and at grid point j is derived from:
(7.50)
where the value of Ψ at half time level (n + ½) and at half grid point (j + ½) is estimated using the Lax scheme:
(7.51)
The resulting approximation to the advection equation:
(7.52)
is second-order accurate in time even though it involves only two time levels. Here again, the third term in the right-hand side of the equation can be viewed as a diffusion term added to the FTCS scheme. In contrast to the Lax method, the effective diffusion coefficient c2Δt/2 is proportional to the time step Δt, so that numerical diffusion can be reduced by adopting smaller time steps.
The amplification coefficient is:
g(k) = 1 − i α sin (kΔx) − α2[1 − cos (kΔx)]
(7.53)
with amplification
factor
|g(k)| = {1 − 4α2[1 − α2]sin4(kΔx/2)}1/2
(7.54)
The stability criterion is again satisfied for α ≤ 1. Some amplitude dampening occurs in the solution for all α < 1, but it is weak for wavelengths that are large compared to the grid spacing Δx. For the smallest wave that can be resolved by the grid (L = 2Δx or equivalently kΔx = π), the amplification coefficient is as low as zero for α2 = 0.5 so the wave disappears. However, for a wave with L = 4 Δx, the dissipation is already considerably smaller; in this case the minimum amplification coefficient is 0.8 for α2 = 0.5.
The phase Φ(k) associated with wavenumber k (see Box 7.2) is deduced from
(7.55)
An important consideration is the nature of the phase errors produced by second-order algorithms such as the Lax–Wendroff scheme. If one applies a Taylor expansion to the trigonometric functions appearing in (7.55), one can show that the phase error for wavenumber k is given by
(7.56)
Modeling of Atmospheric Chemistry Page 36