(7.173)
where
are area integrals in which
z = (y − xj)/Δx
Again, function ψj (z, t) can be approximated by an area-preserving polynomial (7.172). Substitution of this polynomial into the above expressions yields the following expressions for integrals Ij + 1/2+ and Ij + 1/2−:
(7.174)
and
(7.175)
This polynomial fitting method is not exempt from localized unphysical oscillations near sharp spatial gradients. To address the problem, Bott (1989a, 1989b) introduces nonlinear flux limiters and imposes that the total amount of outflow from gridbox j during a time step Δt be limited to Δx/Δt. In addition, the flux Fj+1/2 is set to zero if it does not have the same sign as the velocity uj+1/2. These two conditions are fulfilled if, rather than using (7.173), the flux Fj+1/2 is expressed as
(7.176)
where
and
Here ε is a small value added to avoid numerical unstable situations if .
Spatial functional distributions other than polynomials have been used to interpolate the dependent variables between grid points. Spalding (1972), for example, uses an exponential fitting technique, which prevents the spurious oscillations associated with the Bott scheme near sharp gradients, and ensures positivity of the solution. The method, however, is diffusive and computationally expensive. Chlond (1994) uses a hybrid scheme in which the polynomial and exponential interpolation methods are combined. A switch is used so that the polynomial scheme is applied in regions where the distribution of the transported quantity is smooth, and the exponential fitting technique is applied near sharp gradients.
The Prather Scheme
The algorithm presented by Prather (1986) is an extension of the Russell and Lerner scheme in which the advected function inside a grid cell is represented by a second-order polynomial. In each grid cell, we represent the 3-D (x, y, z) distribution of the tracer mixing ratio ψ by
ψ(x, y, z) = a0 + axx + axxx2 + ayy + ayyy2 + azz + azzz2 + axyxy + ayzxy + axzxz
(7.177)
within a rectangular gridbox of volume V = Δx Δy Δz with 0 ≤ x ≤ Δx, 0 ≤ y ≤ Δy, and 0 ≤ z ≤ Δz. The same function can also be expressed by a linear combination of orthogonal second-order polynomials Kk:
ψ(x, y, z) = m0K0 + mxKx + myKy + myyKyy + mzKz + mzzKzz + mxyKxy + myzKyz + mxzKxz
(7.178)
where mk are moment coefficients. By definition, the orthogonal functions satisfy the conditions:
(7.179)
where dV = dx dy dz. Prather (1986) provides ten orthogonal polynomials that apply to the algorithm in three dimensions:
K0 = 0
with appropriate normalization factors. He also provides linear expressions that relate coefficients ak and mk. An upstream method is used to transport simultaneously the zeroth-order (mass), first-order (slope), and second-order (curvature) moments of the tracer distribution in each grid cell. The moments Si are defined by
(7.180)
(7.181)
(7.182)
(7.183)
Parallel expressions are derived for Sy, Syy, Sxy, Sz, Szz, and Sxz.
To illustrate the method, we consider a 2-D problem and assume that the velocity c > 0 of the flow is uniform and directed in the x-direction. The distribution of the tracer mixing ratio inside a 2-D gridbox is expressed by
ψ(x, y) = m0K0 + mxKx(x) + mxxKxx(x) + myKy(y) + myyKyy(y) + mxyKxy(x, y)
(7.184)
Within a cell j of volume V, we define the sub-volume VR of the fluid that will be removed from this cell and added to the neighboring cell j + 1 over a time step Δt (Figure 7.23):
VR = cΔt Δy Δz
(7.185)
The volume VL of the fluid remaining in the original cell j is
VL = (Δx − cΔt)Δy Δz
(7.186)
The method involves two consecutive steps:
First step: Decomposition of the moments Sk for each grid cell into sub-moments and associated with the fraction of the tracer that is advected to the downwind grid cell and the fraction of the tracer that remains in the grid cell during time step Δt. We have
where the Courant number α = cΔt/Δx = VR/V is assumed to be smaller than 1. The fraction of the tracer remaining in the cell j during time step Δt is located in the sub-box of volume VL extending from bounds xj–1/2 to xj+1/2 – cΔt. The corresponding sub-moments are:
Second step: Advection step and reconstruction of the moments for the complete grid cell. The advection is performed by transporting moments from grid cell j to adjacent cell j + 1, while maintaining moments in their original box. For time tn+1, the moments for the entire gridbox j can be reconstructed by calculating the updated moments:
From these new moments derived on the full grid cell, one derives the coefficients mk using (7.180)–(7.183), and from there the spatial distribution of the tracer mixing ratio inside the grid cell.
Figure 7.23 (a) Advection along direction x of moments during a time step Δt from grid cell j to j + 1 after decomposition of moment S into sub-moments SL, which remains in gridbox j, and SR, which is transferred from gridbox j to adjacent downwind gridbox j + 1. Coefficient α = VR/(VR + VL) = cΔt/Δx is the Courant number. (b) Comparison of the distribution of a tracer advected across 200 gridboxes by the basic upstream scheme (0), the first-order moments method (1) and the second-order moments method (2). Note the presence of overshoots and undershoots in the second-order scheme. Positivity is ensured by placing limits on the high-order moments (curve 2 + L).
From Prather (1986).
The Prather method is less diffusive than the slope scheme (Figure 7.23b), but it adds to the computational and memory requirements because at each time step ten moments must be computed and stored for each grid cell in the 3-D case. As with other high-order methods, the scheme produces overshoots and undershoots that can lead to negative solutions. Placing a limit on the magnitude of the higher-order moments can ensure positivity of the solution (Figure 7.23b). The scheme is absolutely stable for Courant numbers ranging from 0.2764 to 0.7236 (Prather, 1986), but is marginally unstable over the rest of the domain [0, 1]. Phase errors are extremely small.
7.7 Lagrangian Methods
Lagrangian advection methods divide the atmosphere into a large number of air parcels and follow the displacement of their centroids as a function of time. Tracer mixing ratios are conserved in these displacements. The trajectory of the centroid is determined from the wind velocity at the centroid location. Because the wind velocities are generally provided at the discrete points where observations are performed or at the grid points of an Eulerian meteorological model, their values must be interpolated at the location of each centroid. Although the motion of the air parcels does not follow any grid, model results can still be provided at regularly spaced grid points by interpolation from the randomly located air parcels situated in the vicinity of these grid points.
In the general case of a multidimensional model with variable wind velocity v(r, t), where r is the parcel location, the parcel trajectories are calculated from:
(7.187)
A second-order accurate solution of this differential equation is obtained by solving, for example, the implicit expression (4.270):
(7.188)
where r(t) is the location of the air parcel at time t (departure point) and r (t + Δt) is the position of the parcel at time t + Δt (arrival point). This implicit equation can be solved by an iteration and interpolation procedure (Kida, 1983; Stohl, 1998).
The use of computationally expensive iterative implicit methods can be avoided if the position vector r(t) is expanded by a Taylor series in which terms of order higher than (Δt)2 are neglected:
(7.189)
Here, the wind velocity v = (dr/dt)t and the acceleration γ = (dv/dt)t = (d2r/dt2)t are calculated at the departure point (time level t). The acceleration γ is easily obtained from the spatial variation of the velocity field if the local rate of change
in the velocity Δv/Δt can be neglected over the time interval Δt. Thus:
(7.190)
In the above equations, a random wind velocity v′ is often added to the mean wind velocity v to account for small-scale turbulent motions. An example of a simulation produced over a period of three weeks by a Lagrangian particle dispersion model following a volcanic eruption is shown in Figure 7.24.
Figure 7.24 Simulation for the period May 10–21, 2010 of the vertically integrated concentration [mg m–2] of atmospheric ash resulting from the eruption of the Eyjafjalljokull volcano in southern Iceland. Transport was simulated using the Lagrangian particle dispersion model FLEXPART (Stohl et al., 1998), which traces the displacement of a large number of particles by the mean winds to which random motions representing turbulence and convection are superimposed. The model is driven by assimilated meteorological data on a 0.18°×0.18° grid with 91 levels in the vertical.
Reproduced from Papayannis et al. (2012).
In global Lagrangian models, the position of an air parcel is often defined by its longitude λ, latitude φ, and altitude z (or pressure p). Thus, after a time step Δt, an air parcel originally located at point [λ(t), φ(t), z(t)] is displaced to a new point whose position [λ(t + Δt), φ(t + Δt), z(t + Δt)] is derived from
(7.191)
(7.192)
z(t + Δt) = z(t) + Δz
(7.193)
where a is the Earth’s radius. The geometric displacements Δx, Δy, and Δz are expressed as a function of the components (u, v, w) of the wind velocity v and the components (γx, γy, γz) of the wind acceleration γ by
(7.194)
From (7.190), the three components of the acceleration (γx, γy, γz) are expressed by
(7.195)
(7.196)
(7.197)
Lagrangian methods have several advantages over Eulerian methods. First, since all tracers follow the same trajectory, a single calculation of the air parcel displacement can be used to immediately infer the transport of all tracers. Second, the stability of the algorithm is not constrained by the value of the Courant number as in the explicit Eulerian methods, so the adopted time step is limited by accuracy rather than by stability considerations. Other desirable requirements are met: during the parcel displacement, mass is conserved and the sign of the transported function is maintained. Thus, unless interpolation procedures are not carefully performed, the method also guarantees monotonicity, transportivity and locality of the solution.
The Lagrangian methods also have several disadvantages. First, inaccuracies in the interpolation of the wind velocities lead to errors in the calculation of the parcel trajectories, and these errors can accumulate as the time integration proceeds. Second, the initially uniform distribution of air parcels may become highly irregular over time as a result of errors in the wind interpolation. As a result, the tracer concentration may become under-determined in certain parts of the domain while being over-determined in others. Third, Lagrangian transport does not allow for mixing between air parcels even when they are closely located. As a result, contrary to the Eulerian algorithms that often produce excessive diffusion, Lagrangian methods require that some diffusive mixing be added to account for interactions between air parcels. This is critical in particular for the treatment of nonlinear chemistry and aerosol microphysics.
7.8 Semi-Lagrangian Methods
Semi-Lagrangian transport (SLT) methods combine important advantages of the Lagrangian and Eulerian methods. The upstream SLT method (Figure 7.25) consists of using Lagrangian back-trajectories to relate concentrations on a regular Eulerian grid at the end of a model time step to the concentrations at the beginning of the time step. We first consider a version of the SLT scheme in which points in the atmosphere with given tracer mixing ratios are displaced with the flow during a model time step to reach locations coincident with the model grid points. We then consider a finite volume version of the SLT scheme in which volumes of air are displaced with the flow (with no mass transfer through their boundaries) to reach at the end of the model time step a volume of air that is coincident with a model grid cell.
Figure 7.25 Schematic representation of the semi-Lagrangian method in two dimensions. A parcel located at the arrival point A at time level tn+1 was located at the departure point D at time level tn. The value of a conserved quantity such as the mixing ratio of a passive tracer is unchanged as the parcel is displaced from point D to point A during the time step. The value of the quantity at the departure point D is derived by interpolation from the neighboring grid points at time tn. In the finite volume version of the SLT scheme, one considers the displacement of a given volume of air (a surface in two dimensions) from its departure position (blue area) to its arrival position (gray area) coincident with a Eulerian grid cell of the model.
Reproduced with permission from Peter Hjort Lauritzen (personal communication), National Center for Atmospheric Research.
7.8.1 Grid Point Based SLT Schemes
During a time step Δt = tn+1 – tn, one determines the backward trajectory of the atmospheric points that reach the different Eulerian grid points of the model at time level tn+1. The trajectories are calculated using interpolated gridded velocities. Consider an arrival point A at time tn+1 on the Eulerian grid (Figure 7.25). The location rD of the upwind departure point D at time tn is determined by backward integration of (7.187). For example, we write the equation
(7.198)
which has to be solved iteratively since the velocity vD depends on the location of the departure point, which is not known a priori. In general, the departure points do not coincide with model grid points. The tracer mixing ratio at departure point D and time level tn is determined by the interpolation of the surrounding values of the mixing ratio at the closest regular Eulerian grid points. For an inert tracer, the mixing ratio μ remains constant along the trajectory between D and A. Thus:
μ(rA, tn + 1) = μ(rD, tn)
(7.199)
The starting point of the backward trajectory can be derived by replacing the 3-D advection problem with three successive 1-D problems (Seibert and Morariu, 1991). For the advection along the x-axis, the trajectory is determined by integrating
(7.200)
between the departure point xD and the arrival point (coincident with grid point xA = xj; see Figure 7.26). Thus,
(7.201)
Approximating the wind along the trajectory by
(7.202)
(7.201) can be solved analytically and the departure point xD is found to be
(7.203)
where the partial derivative ∂u/∂x is numerically calculated as
(7.204)
In the simple 1-D case in which the wind velocity u = c is constant and positive (see Figure 7.26), the departure point is given by
xD = xA − cΔt
(7.205)
and the mixing ratio at the departure point D is provided for example by a linear interpolation between the values at the closest grid points (indices m – 1 = j – (L + 1) and m = j – L in Figure 7.26)
(7.206)
When the Courant number α = cΔt/Δx is smaller than 1, the departure point D is located between grid point xj–1 and the arrival point xj (A), and it is straightforward to show that
(7.207)
because the mixing ratio remains constant during the displacement of the parcel. Adopting the more classic notation, we write
(7.208)
More generally, if the Courant number α has a ceiling (smallest larger integer) of L, such that α′ = L – α is positive and smaller than unity, then the departure point D is located in the grid cell [j – L, j – L + 1]. The interpolation formula (7.208) then becomes:
(7.209)
Numerical approximation (7.209) is equivalent to expression (7.75) that describes the notoriously dissipative Eulerian upstream method. Thus, if one uses a linear interpolation as implemented here, the semi-Lagrangian method is excessively diffusive. By using higher-order interpolation schemes, the i
ntensity of the diffusion can be considerably reduced. Cubic spline functions (Bermejo, 1990) or biquadratic polynomials (Lauritzen et al., 2010) are often adopted. Williamson and Rasch (1989) and Rasch and Williamson (1990b) examine several possible interpolators and assess their ability to preserve the shape of the advected fields. Accurate interpolation schemes add to the computational costs of the method.
Figure 7.26 Representation of the semi-Lagrangian scheme in the 1-D case for a variable wind speed c(x, t). The trajectory is represented by the curve DA (in red). The arrival point A at time tn+1 is coincident with grid point xj of the Eulerian model grid. The location of the departure point D is derived from a back-trajectory calculation and is not coincident with a model grid point. It is located between grid points xj–(L+1) and xj–L. (In this figure, L = 2, but it could be larger for longer time steps). By approximating the variable wind velocity by its value at the midpoint M (or by the average between the wind speeds at the departure and arrival points), the curve DA (in red) can be approximated by the straight line D′A (in blue), and the departure point D by point D′. The approximation can be improved by iteration.
Modeling of Atmospheric Chemistry Page 40