(6.72)
The specific expressions used in the Gear’s algorithm for different orders (1 to 6) are the following:
This algorithm is very robust, accurate, and, in most cases, stable up to order 6 (see Section 4.8.1). It requires, however, as in the case of the backward Euler algorithm, the resolution of a nonlinear algebraic system. In the Livermore Solver for Ordinary Differential Equations (LSODE) (Hindmarsh, 1977), the solution to (6.72) or equivalently to
is found by applying a Newton-Raphson iterative procedure
where r represents an iteration index and . This leads to the linear algebraic system
(6.73)
that is solved by using, for example, an LU decomposition technique (see Box 6.2). The Jacobian matrix J = ∂s/∂ψ, which appears in the predictor matrix (I−∆t γ J) is usually sparse, so that computationally efficient sparse-matrix techniques can be applied to solve the system. The Jacobian should in principle be re-evaluated at each step of the iteration. In most practical applications, however, it is calculated only at the start of the iteration, or occasionally re-evaluated as the iteration proceeds. Gear?s method offers a strategy to keep the solution error below a user-specified tolerance, by varying the order of the backward differentiation scheme and, when appropriate, by reducing the time step according to the intermediate results of the computation.
References
Byrne G. D. and Hindmarsh A. C. (1975) A polyalgorithm for the numerical solution of ordinary differential equations, ACM Trans. Math. Softw., 1, 71–96.
Dabdub D. and Seinfeld J. H. (1995) Extrapolation techniques used in the solution of stiff ODEs associated with chemical kinetics of air quality models, Atmos. Environ., 29, 403–410.
Gear C. W. (1971) Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, NJ.
Gong W. and Cho H.-R. (1993) A numerical scheme for the integration of the gas phase chemical rate equations in a three-dimensional atmospheric models, Atmos. Environ., 27A (14), 2147–2160.
Hairer E. and Wanner G. (1996) Solving Ordinary Differential Equations II. Stiff and Differential Algebraic Problems, Springer, New York.
Hairer E., Norsett S. P. and Wanner G. (1993) Solving Ordinary Differential Equations I. Nonstiff Problems, 2nd edition, Springer, New York.
Hairer E., Lubich C., and Wanner G. (2002) Geometric Numerical Integration, Structure Preserving Methods for Ordinary Differential Equations, Springer Verlag, Berlin.
Hertel O., Berkowicz R., and Christensen J. (1993) Test of two numerical schemes for use in atmospheric transport-chemistry models, Atmos. Environ., 27A, 16, 2591–2611.
Hesstvedt E., Hov O., and Isacksen I. (1978) Quasi-steady-state-approximation in air pollution modelling: Comparison of two numerical schemes for oxidant prediction, Int. J. Chem. Kinet., 10, 971–994.
Hindmarsh A. C. (1977), GEARB: Solution of Ordinary Differential Equations Having Banded Jacobian, LLNL, Report UCID 30059, Rev. 2.
Jay L. O., Sandu A., Potra F. A., and Carmichael G. R. (1995) Improved QSSA methods for atmospheric chemistry integration. Reports on Computational Mathematics No 67/1995, Dept. Mathematics, University of Iowa.
Press W. H., Teukolsky S. A., Vetterling W. T., and Flannery B. P. (2007) Numerical Recipes: The Art of Scientific Computing, 3rd edition, Cambridge University Press, Cambridge.
Ramaroson R., Pirre M., and Cariolle D. (1992) A box model for on-line computations of diurnal variations in multidimensional models: Application to the one-dimensional case, Ann. Geophys., 10, 416–428.
Rosenbrock H. H. (1963) Some general implicit processes for the numerical solution of differential equations, J. Comput., 5, 329–330.
Sandilands J. and McConnell J. (1997), Evaluation of a reduced Jacobian chemical solver, J. Geophys. Res., 102 (D15), 19073–19087
Sandu A., Verwer J. G., Blom J. G., et al. (1997a) Benchmarking stiff ODE solvers for atmospheric chemistry problems II: Rosenbrock solvers, Atmos. Environ., 31 (20), 3459–3472.
Sandu A., Verwer J. G., Van Loon M., et al. (1997b) Benchmarking stiff ODE solvers for atmospheric chemistry problems I: Implicit vs explicit, Atmos. Environ., 31 (19), 3151–3166.
Santillana M., Le Sager P., Jacob D. J., and Brenner M. P. (2010) An adaptive reduction algorithm for efficient chemical calculations in global atmospheric chemistry models, Atmos. Environ., 44, 35, 4426–4431.
Saylor R. D. and Ford G. D. (1995) On the comparison of numerical methods for the integration of kinetic equations in atmospheric chemistry and transport models, Atmos. Environ., 29 (19), 2585–2593.
Shimazaki T. (1985) Minor Constituents in the Middle Atmosphere, D. Reidel, Norwell, MA.
Verwer J. G. (1994) Explicit methods for stiff ODEs from atmospheric chemistry. In Numerical Mathematics Conference NUMDIFF-7, Halle, Germany.
Verwer J. G., Spee E. J., Blom J. G., and Hundsdorfer W. (1999) A second-order Rosenbrock method applied to photochemical dispersion problems, SIAM J. Sci. Comput., 20 (4), 1456–1480.
Young T. R. and Boris J. P. (1977) A numerical technique for solving stiff ordinary differential equations associated with the chemical kinetics of reactive flow problems, J. Phys. Chem., 81, 2424–2427.
Zhang H., Linford J. C., Sandu A., and Sander R. (2011) Chemical mechanism solvers in air quality models, Atmosphere, 2, 510–532, doi: 10.3390/atmos2030510.
7
Numerical Methods for Advection
7.1 Introduction
The distribution of chemical species in the atmosphere is affected by air motions ranging from the global circulation down to the millimeter scale, at which point molecular diffusion takes over to dissipate kinetic energy. Air motions conserve the mixing ratios of the transported species since the air molecules are transported in the same way as the species. A plume of an inert chemical species transported in the atmosphere may stretch and filament, but it retains its initial mixing ratio until the filaments have become thin enough for molecular diffusion to dissipate gradients.
Representing this conservative transport in an atmospheric model is a major challenge because models cannot resolve the full range of spatial and temporal scales involved. Even if they could, chaotic behavior in solving the equation of motion would prevent a deterministic representation of the flow. Assimilation of meteorological observations can force model winds to approximate the real atmosphere, but only on the large scales of the observational network and at the cost of small-scale numerical noise introduced by the assimilation process.
From a model perspective, it is useful to distinguish between transport by the resolved large-scale winds, which can be simulated deterministically; and transport by the unresolved small-scale winds, which must be represented stochastically. The distinction between large-scale and small-scale is defined by the resolution of the model. Transport by resolved winds is commonly called advection, while transport by unresolved winds is called eddy flow, turbulence, or (in the vertical) convection. This chapter focuses on the numerical schemes used to solve the advection problem. Schemes for smaller-scale unresolved transport are presented in Chapter 8.
Numerical methods should preserve the properties of the continuous partial differential equations (PDEs) that they attempt to approximate. Desirable properties of numerical methods for advection are listed by Rasch and Williamson (1990a), Williamson (1992), and Lauritzen et al. (2011). They include:
1. Accuracy. The solution must be close to the true state.
2. Stability. The solution must not diverge away from the true state.
3. Monotonicity. The solution should not generate spurious maxima or minima in mixing ratios. Since initial conditions for mixing ratios are positive, monotonicity implies positivity of the solution.
4. Conservation. In the absence of sources and sinks, total mass must be conserved during advection. The algorithm should also conserve the second moment (variance) of the advected quantity.
5. Transportivity. Transport should be downwind only.
6
. Locality. The solution at a given point must not be controlled by the concentrations far away from that point.
7. Correlativity. Relationships between species in the flow must be preserved.
8. Flexibility. An advection scheme is most useful if it can be implemented on different grids and at different resolutions; this makes it in particular applicable for adaptive grids.
9. Efficiency. A more computationally efficient algorithm facilitates simulations with higher resolution, over longer periods, and/or involving a larger number of transported species.
In Chapters 1 and 4 we drew a distinction between Eulerian and Lagrangian models for atmospheric transport (Figure 7.1). A Eulerian model solves the advection equation on a fixed reference grid, while a Lagrangian model tracks particles as they are transported in the atmosphere. Both have advantages and disadvantages. A Eulerian model provides a complete solution over the atmospheric domain with regular spatial resolution, but is subject to numerical noise and to stability constraints. A Lagrangian model has no limitations from numerical diffusion or stability, but it has uneven spatial resolution and cannot easily handle nonlinear chemistry. Eulerian and Lagrangian approaches are sometimes combined to benefit from the advantages of each, as in semi-Lagrangian advection schemes.
Figure 7.1 Eulerian and Lagrangian perspectives. In the Eulerian representation (a), the observer is located at fixed points (model grid points) and tracks the change in the calculated state variable ψ (e.g., mixing ratio C) as air parcels move by. In the Lagrangian representation (b), the observer tracks the change in the variable ψ in individual air parcels as they move with the flow.
Reproduced from Lin (2012).
We focus this chapter on the basic approaches to solve the advection equation, including description of some classic schemes. The schemes used in current models of atmospheric composition are based on these classic schemes, but often include refinements that we cannot present in detail.
Section 7.2 presents different forms of the advection equation. Section 7.3 reviews finite difference methods used to solve the advection equation in a Eulerian framework, while Section 7.4 focuses on finite volume methods. Flux-corrected methods introduced to preserve the monotonicity of the solution are discussed in Section 7.5. Selected advanced Eulerian numerical methods are presented in Section 7.6. Sections 7.7 and 7.8 describe the methods used in Lagrangian and semi-Lagrangian models.
7.2 The Advection Equation
An atmospheric species transported in a model is commonly called a tracer. The atmospheric advection of a tracer i is determined by its local mass density ρi and by the wind field v. The corresponding mass flux Fi is the product
Fi = ρiv
(7.1)
As shown in Chapter 4, the local rate of change in the density due to advective transport is the divergence ∇•Fi of this flux. This defines the continuity equation for an inert tracer (no local sources or sinks):
(7.2)
Equation (7.2) is the Eulerian flux form of the advection equation, previously discussed in Chapter 4. Transport in a Eulerian model is often expressed in terms of mass fluxes across grid cell interfaces. This can be expressed by integration over a finite volume Vc (usually a model grid cell) of the advection equation (7.2):
(7.3)
The first term in this integral equation represents the time evolution of the mean density ⟨ρi⟩ of tracer i in the finite volume Vc. The second term can be transformed by applying the divergence (or Gauss–Ostrogradsky) theorem (see Appendix E), which states the equivalence between the volume integral over Vc, and the surface integral over the closed boundary of volume Vc. Equation (7.3) then becomes:
(7.4)
where Fi represents the flux vector of tracer i across surface Sc of the boundary of the cell and n is a unit outward vector normal to the cell boundary. This expression states that the change of the mean density inside the finite cell is determining by the net flux of material in and out of the cell. Equations (7.2) and (7.4) are often called conservative forms of the advection equation.
When applied to a 2-D (e.g., horizontal) model, (7.4) becomes
(7.5)
where the integral is now calculated along the boundary line Lc of the 2-D cell with area Ac. Figure 7.2 illustrates the different elements needed to solve (7.5) in the case of a 2-D hexagonal grid cell. Finally, in the case of a 1-D problem (along direction x), the expression becomes
(7.6)
where Δx represents the size of the 1-D grid cell. and are the tracer fluxes at the right and left edges of the cell (positive rightward). As stated in Chapter 4, these expressions constitute the basis for finite volume methods; they are particularly suitable when solving the equations on complex or irregular grids.
Figure 7.2 Formulation of tracer advection using (7.5) in the case of a 2-D hexagonal grid cell with surface Ac and boundaries Lc. The flux Fi across a particular side of the hexagon is schematically represented. Unit vector n is normal (outward) to the cell boundary and ⟨ρi⟩ is the average density of tracer i in the cell.
When expressed in terms of mass mixing ratio μi = ρi/ρa where ρa is the density of air, the continuity equation becomes (advective form):
(7.7)
or
(7.8)
where the total derivative operator (derivative along the flow) is expressed by:
(7.9)
These forms were previously derived in Chapter 4. Equation (7.8) specifies the invariance of the mixing ratio along flow trajectories. By contrast, tracer densities ρi along flow trajectories may change because air is a compressible fluid (see example below). Equation (7.7) is the Eulerian advective form of the advection equation, while (7.8) is the Lagrangian form.
Solution of the advection equation in 3-D models involves operator splitting along individual dimensions. Thus, the 1-D advection equation is solved successively in the three dimensions to obtain the 3-D solution. We will focus our discussion in the following sections on the 1-D problem as it is most relevant for model applications. In this case, the flux-form advection equation (7.2) becomes:
(7.10)
and the advective-form equation (7.7) becomes:
(7.11)
where u(x, t) denotes the 1-D wind velocity. Rewriting the flux-form equation (7.10) as
(7.12)
shows that it is identical to the advective form (7.11), but with an additional term ρi ∂u/∂x that describes the compressibility of the flow. This term acts as a source of ρi when u(x, t) decreases with x (compression) and as a sink when u(x, t) increases with x (expansion). If the wind velocity is constant in space and time (u(x, t) = c), the flux form of the advection equation is simply
(7.13)
Thus the advection equation applies identically to density and mixing ratio in an incompressible flow.
Figure 7.3 shows the solution of the 1-D advection equation for the two cases of a constant and a spatially variable wind velocity. For a constant velocity, the advection of both the density and the mixing ratio is represented by a simple translation (without deformation) of the initial function in the direction of the velocity. If the velocity decreases with space, the initial distribution of both quantities is distorted as the material is advected. The value of the maximum mixing ratio is unchanged, but the maximum value of the density is enhanced. Advection can thus modify extrema of tracer densities in a diverging flow.
Figure 7.3 One-dimensional advection in the x-direction of a scalar function, noted here q(x, t), for a constant velocity u (left) and a spatially varying velocity u(x) (right). The initial distribution of the function is represented by the solid line. On the left panel, this initial function (density or mixing ratio) is translated in the direction of the constant velocity u. On the right panel, the dotted curve (labeled pure advection) corresponds to a case where the term –q ∂u/∂x is omitted and is therefore representative of the advection of a tracer mixing ratio. In this case, the spatial distribution of the function is modified, but the maximum value of th
e function is unchanged. The dashed curve is obtained by including the term –q ∂u/∂x and is therefore representative of the advection of a tracer density. In this case, the area under the curve is maintained during the advection process, but the maximum value of the function is not preserved.
From C. P. Dullemond with permission.
In many applications, tracers are not only advected but are also diluted by turbulent diffusion. The 1-D advection–diffusion equation is given by:
(7.14)
where K [m2 s–1] is a diffusion coefficient. The relative importance of advection versus diffusion is measured by the dimensionless Péclet number Pe:
(7.15)
where L is a characteristic length for the advection. The Péclet number can be regarded as a measure of the ratio of the diffusive timescale to the advective timescale. For conditions typical of horizontal flow with a wind velocity c ~10 m s–1, a characteristic length L for long-range transport ~1000 km and an eddy diffusion coefficient K of 105 m2 s–1, the Péclet number is ~100 and the transport is thus dominated by advection. For vertical transport in the boundary layer with a typical vertical velocity c ~0.01 m s–1, a characteristic length L of 100 m, and an eddy diffusion coefficient of 100 m2 s–1, the Péclet number is equal to 10–2 and diffusion becomes dominant. The solution of the transport equation for a boxcar function (shock front) subject to simultaneous advection and diffusion under a Péclet number of approximately 1 is shown schematically in Figure 7.4. Note the gradual deformation of the shape of the function under the influence of diffusion. The area under the curve, however, is conserved.
Modeling of Atmospheric Chemistry Page 35