Figure 7.11 Stencils describing several numerical algorithms for the approximate solution of the 1-D advection equation.
Figure 7.12 Two-dimensional advection of a square wave of width 20Δx and initial concentration = 1 advected over 50 grid cells in each coordinate direction for a Courant number of 0.25. (a) Exact solution; (b) first-order upstream; (c) leapfrog; (d) QUICKEST; (e) Lax–Wendroff with flux limiters, (f) MPDATA.
Reproduced from Gross et al. (1999) with permission from the American Society of Civil Engineers (ASCE).
In summary, first-order methods such as the upstream algorithm are characterized by numerical diffusion in the solution and, as a result, tend to reduce the amplitude of peaks and to smooth spatial gradients that are present in the initial tracer distributions. High wavenumber components are also eliminated. Dispersion, which is common to the simple high-order methods described above, tends to distort the solution since it causes all waves, and specifically the small waves, to lag the true solution. Numerical diffusion is generally viewed as a lesser evil because it merely causes loss of information, while dispersion generates spurious information.
7.3.5 Generalization to Variable Wind Speed and Grid Size
The previous discussion has highlighted some fundamental properties of different Eulerian algorithms. In practical applications, the wind in the x-direction may not be uniform (c is replaced by u(x, t)), and the discretization interval Δxj may vary along the spatial dimension x. In that case, the space derivative in the conservative 1-D flux-form equation:
(7.101)
is replaced by its second-order centered finite-difference approximation:
(7.102)
where
and Δxj = xj − xj − 1, Δxj + 1 = xj + 1 − xj. This center-difference scheme can be applied in the case of explicit, implicit, or Crank–Nicholson algorithms.
7.3.6 Mass Conservation
If we integrate the 1-D advection equation (7.17) over the spatial interval [A, B] and between time levels tn and tn+1, we find the conservation expression
(7.103)
where we have again assumed c to be fixed. Condition (7.103) must be met for tracer mass to be conserved with ΨA and ΨB representing boundary conditions. In the simple case of the FCTS scheme, it is easy to evaluate the left-hand side of this integral relation:
(7.104)
or
(7.105)
if
Thus, under the conditions adopted here, the accumulation of the conservative tracer Ψ in the domain [A, B] is proportional to the net flux (c Ψ) at the boundaries A and B. The finite-difference analog has preserved the integral expressed by the continuum equation (7.103). If the fluxes across the external boundaries are zero or if the domain is periodic with ΨA = ΨB, mass is fully conserved in the domain. If the constant velocity c is replaced by a velocity u(x, t) that varies with space and time, tracer conservation will be generally obtained if one considers the finite difference analog of the flux-form equation (7.2) but not its advective form (7.7).
7.3.7 Multidimensional Cases
The 1-D advection problem can be generalized to multiple spatial dimensions. In a 2-D Cartesian space (x, y), the flux-conservative form of the advection equation is expressed by
(7.106)
where u and v are the velocities of the wind components in the x and y directions, respectively. If we assume constant grid spacing, the discretization leads to the following expression
(7.107)
where Δx and Δy are the grid spacings in the x and y directions, and the indices i and j refer to these two directions respectively. If the differences in space are estimated at time tn (explicit method), the algorithm is unconditionally unstable, but it can be made stable by adding numerical diffusion as in the 1-D case. If estimated at time tn+1 (implicit method), the algorithm is unconditionally stable. Expression (7.107) can be extended to three dimensions and solved for variable grid spacing. A Crank–Nicholson form of (7.107) can also be easily derived. The matrix of the system becomes rapidly very large and is not banded as in the 1-D case. Three-dimensional advection in practical applications is usually performed by operator splitting, with successive numerical solution of the 1-D advection equation over each dimension for individual time steps.
7.3.8 Boundary Conditions
The resolution of hyperbolic equations such as the advection equation applied to a limited spatial domain requires that a condition be imposed at the boundary through which material flows into the domain. In the 1-D case with a domain [a, b], the condition must be specified at point x = a if the velocity is positive (c > 0) and at x = b in the opposite situation (c < 0). In some cases, the numerical algorithm requires that an additional condition be provided at the outflow boundary. By applying such a condition without precaution, the problem becomes ill-posed, and the algorithm may provide unstable solutions. This is the case when centered differences are used to represent space derivatives. Consider a 1-D flow on a domain [a, b] with a positive constant velocity c and an inflow boundary condition Ψ(a, t) = H(t) at location x = a. We approximate the spatial derivative by the leapfrog (CTCS) scheme (j = 1, J):
(7.108)
At grid point j = 1, the solution is easily calculated:
(7.109)
since is specified. The calculation of the solution at grid point j = J,
(7.110)
is not straightforward because, in a well-posed problem, no value should be imposed at point J + 1. If a value is nevertheless imposed at this outflow boundary, e.g., or , the scheme will produce unrealistic (unphysical) reflections that propagate upstream. An illustration is provided by Figure 7.13 that shows the advection of a bell-shaped function by a leapfrog scheme with Ψ = 0 at both limits of the domain. As the signal reaches the downwind boundary, spurious wave reflections (saw-toothed artifacts) are produced and propagate upstream. This is generalizable to multidimensional problems: Spurious reflections often occur at the lateral boundaries of a limited-domain nested model driven by boundary conditions from a larger-domain model. Even if the imposed boundary conditions verify the analytic solution of the advection equation, some reflections are to be expected since the numerical solution is slightly different.
Figure 7.13 Advection of a bell-shaped function (arbitrary units) by a leapfrog scheme over a domain of 1000 km with a Courant number of 0.1 (c = 20 m s–1, Δt = 50 s, Δx = 10 km) and boundary conditions of zero at the edges of the domain. Initial condition (a), solution after six hours (b), after eight hours (c) and after ten hours (d).
From P. Termonia, Royal Meteorological Institute of Belgium.
In practical atmospheric applications, the sign of the wind velocity and hence the direction of the flow at the boundary of the domain frequently change as the model simulation proceeds. It is therefore difficult to identify the boundary at which a condition should be specified. This problem is usually addressed by imposing some conditions along the entire boundary of the model domain (ill-posed condition), while adding in the equations a relaxation term that damps the high-frequency signal produced at the downwind boundaries. In this case, the original advection equation is modified as (Davies, 1983):
(7.111)
where the relaxation coefficient λ(x) is different from zero only in the boundary zones (a few grid cells near the inflow and outflow boundaries, called buffer zones) and is an externally specified field chosen to be close to the expected solution. The value of λ(x) and the width of the relaxation zone (typically 2Δx to 6Δx) need to be optimized to avoid the reflection of waves while minimizing perturbation to the solution. To ensure stability, the relaxation term should be estimated at time tn+1.
An alternative damping scheme is to add a diffusion term to the advection equation:
(7.112)
where the diffusion coefficient K(x) is non-zero only near the boundary zones. Davies (1983) discusses the stability conditions for this approach.
7.4 Elementary Finite Volume Methods
In finite volume approach
es, rather than considering the values of function Ψ at specified points of a model grid, one calculates the average of this function over defined grid cells. The grid points are now viewed as the centers of grid cells, often called gridboxes. The cell boundaries are called cell edges, walls, or interfaces.
7.4.1 One-Dimensional Formulation
In the 1-D problem (x-direction), the location of the cell center is noted xj, while the locations of the cell interfaces are noted xj–1/2 (left side) and xj+1/2 (right side). For each grid cell (j) (whose size is assumed here to be constant and equal to Δx):
(7.113)
except at the left (j = 1) and right (j = N) boundaries of the model, where we adopt
(7.114)
The average value Ψj of the variable distribution ψ(x, t) inside the cell is
(7.115)
If ψ(x, t) is a tracer concentration, then Ψj is the mean tracer concentration and Ψj Δx the tracer mass in grid cell (j).
The exact (analytic) solution of the 1-D advection equation with fixed wind velocity c, when integrated over a time period Δt, is
ψ(x, t + Δt) = ψ(x − cΔt, t)
(7.116)
or, when the integral form is considered instead,
(7.117)
Recognizing that the first term in (7.117) is equal to , and defining x ' = x − cΔt, one can write:
(7.118)
Splitting this integral into different contributing parts, one writes equivalently:
(7.119)
or
(7.120)
The mean tracer concentration in grid cell j at time tn+1 is thus obtained by adding the mean value of the tracer concentration that enters grid cell (j) to the existing mean value in that cell at time tn and removing the mean value that is transported downstream from cell (j) to cell (j + 1). In this expression, the mass leaving the upwind donor cell equals the mass entering the neighboring downwind receptor cell. The finite volume method is therefore perfectly mass-conserving, which is its main advantage.
If and represent the mean fluxes through the left and right interfaces of grid cell j, respectively, averaged over time Δt = tn + 1 − tn, we write equivalently to (7.120):
(7.121)
For the period Δt during which the subgrid function ψ(x) is assumed to remain unchanged, donor cell (j – 1) transfers a mass to receptor cell (j):
(7.122)
where the interval [xj–1/2 – cΔt, xj–1/2] corresponds to the shaded area in cell (j – 1) (Figure 7.14). Similarly, the mass transferred from donor grid cell j to receptor grid cell (j + 1) is
(7.123)
Figure 7.14 Representation of the donor-cell scheme in one dimension x. During time step Δt, the shaded area in cell j is transported in the x-direction to cell j + 1.
Different implementations of the finite volume method (i.e., different assumptions for the subgrid distribution of ψ(x) inside each cell) lead to different estimates of the fluxes at the interfaces of the grid cells. The simplest assumption is that ψ(x) is uniform inside each cell, so that the state Ψj at the grid center is identical to the state everywhere in the grid. This is the donor cell method. If the subgrid function ψ(x) is assumed to vary linearly with position x inside each grid cell, the algorithm is called piecewise linear. If ψ(x) is a second-order polynomial the algorithm is called quadratic or piecewise parabolic; see examples in Figure 7.15.
Figure 7.15 Representation of the spatial distribution of a tracer within four grid cells: zeroth-, first-, and second-order polynomials.
Donor-cell algorithm
In the simple donor-cell scheme with c > 0, the subgrid function ψ(x) is uniform inside each cell, and the mass Fj–1/2 Δt advected from cell (j – 1) to cell (j) over a time period Δt is equal to Ψj − 1 c Δt = Ψj − 1 α Δx, Simultaneously, a mass equal to Fj + 1/2 Δt = Ψj c Δt = Ψj α Δx is displaced from cell (j) to cell (j + 1). From (7.121), we find:
(7.124)
Similarly, if c < 0, we have
(7.125)
As shown by Figure 7.14 displayed for c > 0, term accounts for the mass transferred from the donor grid cell during the time step Δt, while term Ψj (1 – α) Δx represents the mass that remains in cell j during this time period. The addition of these two terms represents the resulting mass in gridbox (j) at time tn+1. This expression, derived in the simple case where ψ(x) is assumed uniform, is identical to the upstream formula (7.74). The method is first-order accurate and is therefore characterized by large numerical diffusion.
Piecewise linear algorithm
In the case of the piecewise linear approach, we write for xj − 1/2 < x < xj + 1/2:
ψ(x) = Ψj + bj(x − xj)
(7.126)
where Ψj denotes the value of linear function ψ(x) at the center of the cell (also the mean value of ψ(x) in the cell), and bj is the slope of the function inside the cell. For a fixed velocity c > 0 and an equally spaced grid, the flux at the left cell interface is:
(7.127)
or
(7.128)
Similarly, for c < 0, one finds
(7.129)
From expression (7.119), one can easily deduce that the average value of subgrid function ψ(x) in cell (j) at time tn+1 is given by:
(7.130)
if c > 0, and
(7.131)
if c < 0. This algorithm can be viewed as an extension of the donor-cell scheme with a correction term that disappears if the slope of function ψ(x) is equal to zero.
The value of bj is expressed as a function of the values of function Ψ at the center of adjacent cells. Different options are: (1) centered slope or Fromm method, (2) upwind slope or Beam–Warming method, and (3) downwind slope (equivalent to the Lax–Wendroff algorithm). The values of the bj coefficients are respectively
The resulting algorithm for c > 0 is in the case of the Fromm scheme (which is upwind biased)
(7.132)
In the case of the Beam–Warming scheme (also upwind biased), it is (see also expression 7.82)
(7.133)
and, in the case of the Lax–Wendroff scheme, we find the three-point stencil (spatially centered) expression that is identical to (7.52)
(7.134)
Figure 7.19 in Section 7.5 shows the numerical solution for advection of a step function obtained with the second-order accurate Fromm and Beam–Warming methods. The solution from the Lax–Wendroff algorithm was previously shown in Figure 7.10. In all three cases, the solution is not monotonic.
Other implementations of the finite volume method with specific subgrid distributions of function ψ(x) (e.g., the algorithms of Russell and Lerner, 1981; Colella and Woodward, 1984; Prather, 1986) are discussed in Section 7.5.
7.4.2 Two-Dimensional Formulation
When extended to two dimensions (see Figure 7.16), the finite volume algorithm is expressed as:
(7.135)
where F and G represent the mean fluxes in the x- and y-direction respectively. The challenge in defining accurate algorithms is to properly formulate the fluxes at the interfaces as a function of the dependent variables in the neighboring cells.
Figure 7.16 Representation of orthogonal flux components F (in the x-direction) and G (in the y-direction) across cell interfaces in two dimensions. The flux form adopted for the algorithm ensures mass conservation.
Another approach is to solve the discretized form (7.6) for the mean density ⟨ρ⟩j inside grid cell j (Dukowicz and Baumgardner, 2000; Lipscomb and Ringler, 2005; Miura, 2007; Skamarock and Menchaca, 2010):
(7.136)
Here, Aj is the area of the cell, Fj the mass flux across the interfaces of the cell, n a unit vector perpendicular to the cell boundaries, and Lc the length of each cell edge. The sum applies to all cell edges. The flux can be estimated through a remapping algorithm, as depicted in Figure 7.17 for a hexagonal cell where the fluid velocity v at one point of each cell boundary (e.g., center of the cell edge) is projected backward to define the upstream flux-area Am (shaded parallelogram). The me
an density in area Am is derived by a polynomial fit using the mean densities in the neighboring cells at time level tn. The flux Fj is then derived from the mass contained in area Am that is displaced across the cell edge over a time interval Δt with velocity v. The mass originating from all neighborhood cells is remapped onto cell j and provides the mean density in this cell at time level tn+1. The accuracy of the scheme depends on the order of the polynomial that is adopted. In the incremental remapping of Dukowicz and Baumgardner (2000) and of Lipscomb and Ringler (2005), all endpoints on the grid are tracked back, so that the upstream flux area is a polygon. The method bears some similarities with semi-Lagrangian schemes discussed in Section 7.8.
Figure 7.17 Schematic representation of the 2-D remapping algorithm of Miura (2007) in the case of a hexagonal cell. The shaded region represents the mass advected through the cell boundary over a time step Δt.
Redrawn from Skamarock and Menchaca (2010).
7.5 Preserving Monotonicity: Flux-Corrected Transport
As stated in Section 7.4, the solutions provided by high-order accurate algorithms are not monotonic. Preserving monotonicity in the solution of the advection equation is an important requirement for chemical transport models. The generation of new extrema or “ripples” in the vicinity of steep gradients (including shocks and discontinuities of the solution) is unacceptable in most applications. Correction techniques have therefore been proposed to eliminate these unphysical maxima or minima caused by numerical dispersion in high-order algorithms.
Modeling of Atmospheric Chemistry Page 38