Ψtrue(tn + 1) = exp [(A1 + A2)Δt] Ψ(tn)
Adopting a splitting method, the equation is solved over the splitting time interval Δt by treating the two processes A1 and A2 sequentially:
We find
Ψsplitting(tn + 1) = exp [A2Δt] exp [A1Δt] Ψ(tn)
The local error e (error made during one single time step Δt) due to the splitting process is
e(Δt) = Ψtrue(tn + 1) – Ψsplitting(tn + 1) = {exp[(A1 + A2)Δt] – exp [A2Δt] exp [A1Δt]} Ψ(tn)
which is zero only if the two matrices commute (A1 A2 = A2 A1). If the exponentials are approximated to the second-order as
where I is the identity matrix, it is easy to show that the local splitting error is of order O(Δt2) if the two processes do not commute. By summing over the entire integration time, one finds that the global splitting error associated with the above algorithm is of the order O(Δt), and the method is therefore said to be a first-order method.
When one of the two operators (say, B) is stiff, the A–B scheme described above can be slightly modified to avoid transient disturbances (Sportisse, 2000). In this case, rather than imposing an initial condition to the second substep that accounts for the calculation of the first substep, one prefers to add a source term during the second substep and replace (4.320) by
(4.321)
This scheme represents an explicit integration with the non-stiff operator A and an implicit integration with the stiff operator B.
The algorithm described above is first-order accurate (on global error, see Box 4.12) with respect to the splitting time Δt. The splitting error can be reduced by adopting a symmetric splitting scheme. One option, which is computationally expensive, is to calculate at each time step the average between the solutions derived for different possible orderings of the operators. When only two operators A and B are considered, the solution is simply the average between the values calculated by the A–B and B–A schemes:
(4.322)
The local error associated with the solution is fourth-order.
A more efficient approach introduced by Strang (1968) and presented here in the simple case of two operators A and B is to integrate the equations as a symmetric sequence of operators A–B–A. The integration is performed first for operator A over the time interval [0, Δt/2], then for operator B over the full interval [0, Δt] and finally for operator A again over the interval [0, Δt/2]. Thus
(4.323)
(4.324)
(4.325)
This Strang approach can be generalized to a larger number of operators. For example, advection (A), diffusion (D), and chemical (C) in the sequence AΔt/2DΔt/2CΔtDΔt/2AΔt/2 where the subscript denotes the time interval over which the integration is performed.
Strang splitting leads to a second-order approximation. There is no operator splitting error if the operators commute, but this is not generally the case in atmospheric applications. For example, advection commutes with diffusion only if either the wind or diffusion fields do not vary in space. Diffusion commutes with chemistry only if the chemical source term is linear in concentration and independent of the spatial variable (Lanser and Verwer, 1998).
The computational cost of solving the 3-D transport equation on a grid of N points is generally proportional to N3, but decreases to about 3N if the problem is reduced to a set of three 1-D equations by operator splitting. In the case of pure advection, where the flux divergence is written as
(4.326)
the 3-D advection equation can be solved by three sequential 1-D operators applied to the 3-D concentration field represented by vector Ψ:
Ψ(1) = Ax(Ψn)
(4.327)
Ψ(2) = Ay(Ψ(1))
(4.328)
Ψn + 1 = Az(Ψ(2))
(4.329)
Here the advection operators Ax, Ay, Az solve the corresponding 1-D advection equations. For Ax, for example, the advection equation is
(4.330)
and Ax(Ψn) solves this equation over the time interval [tn, tm+1]. Numerical methods for solving the 1-D advection equation are presented in Chapter 7. A side advantage of splitting 3-D advection into 1-D operators is that it allows the use of different time steps in different directions to address CFL criterion limitations. With stronger winds in the x-direction than in the y-direction, and weak winds in the z direction, one can use for example the following arrangement of operators:
Ψn + 1 = Ax , Δt/4Ay , Δt/2Ax , Δt/4Az , ΔtAx , Δt/4Ay , Δt/2Ax , Δt/4(Ψn)
(4.331)
Diffusive transport over a model time step can also be decomposed into three different operators, one in each direction. Here, numerical stability requires that an implicit approach be adopted (Section 8.6). Each direction can be treated completely separately, as in the advection case. However, it is usually preferable to adopt an alternating direction implicit (ADI) method in which, during a third of the time step, diffusion is solved implicitly in one direction and explicitly in the other two directions. Consider the simplest case of a constant diffusion coefficient K and constant air density. The diffusion term takes the form
The ADI method updates the value of the 3-D gridded field Ψ from Ψn to Ψn+1 over a time step Δt = tn+1 – tn with the sequence
(4.332)
(4.333)
(4.334)
where
γ = KΔt
(4.335)
and Dx applied to grid point Ψi,j,k is
(4.336)
with equivalent forms for Dy and Dz. The ADI method is easily generalized to cases in which the diffusion coefficient is variable.
4.15 Filtering
Meteorological models tend to produce undesirable noise caused by dispersion errors in the numerical integration of the dynamical equations. This noise affects the atmospheric distributions of chemical species transported by the model. In order to keep model simulations numerically stable, some form of dissipation may need to be introduced (Jablonowski and Williamson, 2011). Numerical diffusion in solving the advection equation in Eulerian models (Chapter 7) often serves as an implicit filter to damp undesired noise. Filters are explicit if they are implemented by the addition of terms in the governing equations or if they are applied a posteriori as a correction to the calculated fields. Here we describe explicit filters that damp the noisy small-scale waves with wavelengths L of the order of (2–4) Δx without reducing substantially the amplitude of the better-resolved scales.
As we will see, numerical filters generally involve the application of a diffusion term to the solutions of the dynamical equations in order to damp small-scale features. Aside from ensuring numerical stability, an additional purpose of this damping is to mimic turbulence-related processes that are unresolved by the model grid. The introduction of a numerical filter as diffusion operator can thus be based on physical as well as numerical considerations.
4.15.1 Diffusive Filters
A diffusive filter involves the addition of a diffusion term to the dynamical equations as
(4.337)
where q is a positive integer, 2q the order of the diffusion, and K2q the adopted diffusion coefficient [with units m2q s–1]. Setting q = 1 corresponds to the second-order diffusion previously introduced in Section 4.2.3 to parameterize subgrid turbulent mixing. It is also often applied as an artificial sponge in the upper layers of dynamical models to avoid spurious reflection of waves at the top boundary. Second-order filters are not very scale-selective and may negatively impact the well-resolved waves produced by the model. More scale-selective hyper-diffusion schemes with higher values of q can be used. The fourth-order hyper-diffusion scheme with q = 2, called bi-harmonic diffusion or super-viscosity, is often adopted. The chosen value of the diffusion coefficient is somewhat arbitrary and is often regarded as a tuning parameter of the model; it has to be as small as possible to avoid dissipation of well-resolved physical waves, and large enough to ensure numerical stability of the computed solution.
4.15.2 Digital Spatial Filters
Digital spatia
l filters are local filters that take into account only neighboring grid points. A widely used digital filter is the linear Shapiro filter (Shapiro, 1970; 1975) that is based on constant-coefficient grid point operators of order m. The order of the filter determines the width of the numerical stencil (i.e., the number of neighboring grid points involved in the filtering operator). The Shapiro filter eliminates short waves from the calculated fields and thus functions equivalently to the addition of a diffusion term in the dynamical equations. In the 1-D case, we consider a function Ψ(x) defined over the interval (–∞ < x < +∞) with values Ψi at discrete points xi = iΔx. The smoothed value {Ψi} at point xi using a first-order Shapiro filter is the weighted average between the unsmoothed value Ψi at point xi with weight (1 – S) and the two adjacent unsmoothed values at points xi–1 and xi+1 with weight S/2. Thus,
(4.338)
or equivalently,
(4.339)
where S is the so-called smoothing element. We note three important properties of this filtering operator: (1) it is symmetric in space; (2) it involves only three values of x (the filter is local); and (3) over a large number of points, the averages of the smoothed values approach the averages of the unsmoothed values.
To derive the properties of the filter, we express at grid point xi = iΔx the Fourier component of functions Ψ and {Ψ} for wavenumber k [m–1] corresponding to wavelength L = 2π/k, as
Ψk(xj) = Ak exp [ikxi]
(4.340)
and
{Ψk(xi)} = {Ak} exp [ikxi]
(4.341)
The response function g(k) = {Ak}/Ak of the filter is
(4.342)
The value of S must be chosen so that the response function is positive and smaller than 1. This requires 0 ≤ S ≤ ½. In most applications, the value S = ½ is adopted, and the first-order filter operator is expressed by
(4.343)
with the corresponding response function
(4.344)
The filter can be applied more than once to achieve greater smoothing of the solution.
One can also use higher-order filters whose amplitude response are provided by
(4.345)
where the order of the filter m is an integer multiple of 2. Define the difference operator δ as
δΨi = Ψi + 1/2 − Ψi − 1/2
(4.346)
with, for example,
δ2Ψi = δ(δΨi) = δΨi + 1/2 − δΨi − 1/2 = Ψi + 1 − 2Ψi + Ψi − 1
(4.347)
The filter operators of orders 1, 2, 4, …, m are defined successively as
(4.348)
(4.349)
(4.350)
…….
(4.351)
For example, the filter operators of orders 1 and 2 are
(4.352)
(4.353)
Shapiro filters totally eliminate the shortest resolvable wave corresponding to two grid cells (L = 2Δx), and damp to a lesser extent the amplitude of the other small-scale resolvable waves (L = 3Δx and 4Δx).
Figure 4.29 shows the response function of the Shapiro filter for different orders m after 1 and 1000 applications, respectively. It highlights the cumulative character of the filtering operation, specifically in the case of the low-order filters, which strongly damp relatively large waves (up to the 12Δx wave in the case of a second-order filter with 1000 applications). The order-8 filter is often adopted since, after a large number of applications, it effectively eliminates all waves with wavelengths less than 4Δx, but preserves almost entirely the waves with wavelengths larger than 6Δx.
Figure 4.29 Response function for 1-D Shapiro filters of order 2, 3, 4, 8, and 16 with S = ½ after (a) one application and (b) 1000 applications.
Reproduced from Jablonowski and Williamson (2011).
If the first-order Shapiro filter is applied twice in sequence with smoothing factors S equal to +½ and –½, respectively, we find an intermediate value
(4.354)
and a final filtered value
(4.355)
The resulting response function of the two-stage operator is
(4.356)
This particular filter, referred to as the Shuman filter (Shuman, 1957), is very effective for attenuating short waves while preserving large waves: the L = 2Δx wave is totally eliminated and the amplitude of the response function is equal to 0.75 for L = 4Δx and 0.98 for L = 8Δx. The effect of the filter is very small for wavelengths L > 8Δx.
A 2-D function Ψ(x, y) can be filtered by applying the 1-D Shapiro operator successively in directions x and y (with the corresponding indices i and j). The resulting operator involves nine discrete points:
(4.357)
with a response function
(4.358)
Here kx and ky represent the wavenumbers in the x and y directions, respectively. An alternative is to define the 2-D smoothing function as
(4.359)
which results in a five-point operator
(4.360)
with a response function
(4.361)
In the preceding discussion, we have assumed that the spatial domain is infinite and we have therefore ignored the influence of the domain boundaries. This is not of concern if the domain is periodic, as in global model applications. For limited domains, values of the field Ψ imposed at the boundaries may strongly influence the filtered function inside the domain, especially if the filtering is repeated a large number of times. For a simple three-point operator in one direction, the influence of the boundary propagates to the interior of the domain by one grid point from both ends at each application of the filtering process. Shapiro (1970) discusses the boundary effects in detail and shows that, for example, a successful filtering procedure with periodic boundaries may result in the spurious growth of undesired waves inside the domain if fixed conditions are imposed at the boundaries. Some adjustments in the filtering procedure can be applied to limit the influence of the boundary conditions.
4.15.3 Spectral Filters
Spectral filters are commonly used in global longitude–latitude grid point models to damp noise in the vicinity of the pole where convergence of the meridians reduces the longitudinal spacing Δx between grid points. This reduced spacing can lead to violation of the CFL criterion and numerical instability (see Chapter 7). Short waves resulting from such numerical instability are damped or eliminated by applying a 1-D Fourier filter in the zonal direction. The grid data are first transformed into the spectral space via Fourier transform, and the resulting Fourier coefficients am corresponding to dimensionless wavenumber m are modified to become
{am} = F(m) am
(4.362)
with F(m) being the so-called response function. In practical applications, the Fourier filter is applied only poleward of a cut-off latitude φc, and the strength of the filter is gradually enhanced toward the pole. This can be accomplished by increasing the number of wavenumbers m affected by the filtering process and by choosing a response function whose value decreases with latitude. Figure 4.30 shows two examples of response functions; the first one expressed by
(4.363)
corresponds to a strong filter and the second one
(4.364)
to a weaker filter. Here, φ denotes latitude and Δλ the angular longitudinal resolution. The positive integer parameter q can be chosen to modify the strength of the filter.
Figure 4.30 Response function of a strong Fourier filter with q = 1 (solid line) and a weaker (dashed line) filter represented as a function of the dimensionless zonal wavenumber at latitudes of 45°, 60°, and 85°, respectively. The cut-off latitude is assumed to be 40°.
Reproduced from Jablonowski and Williamson (2011).
In the final step of the filtering process, the fields are converted back into the grid point space by an inverse Fourier transform. The advantage of the Fourier filter is that it can be made very scale-selective and dependent on latitude; the drawback is that all data along latitude rings are needed. The use of local spectral filters has be
en proposed for models that utilize local spectral methods like the discontinuous Galerkin approach or the spectral element method (Section 4.10). See Vandeven (1991) and Boyd (1998) for more details.
4.15.4 Time-Smoothing Filters
In certain cases, it is appropriate to apply a slight time smoothing to the solution of a differential equation. In the Robert–Asselin filter (Robert, 1969; Asselin, 1972), the solution Ψn at time tn is replaced by
{Ψn} = Ψn + ν({Ψn − 1} − 2Ψn + Ψn + 1)
(4.365)
and the high frequencies in the solution are damped. Note that the value of the function Ψ at time tn–1 is replaced by its already filtered value. Such a filter with a parameter ν of the order of 1% is often applied when the leapfrog advection scheme is used, since it removes the odd–even checkerboard pattern generated by that numerical algorithm.
4.16 Interpolation and Remapping
Models quantify the values of variables at discrete locations (grid points in Eulerian models, location of individual particles in Lagrangian models) and discrete time steps. One often needs model information intermediate between these discrete points. If the value of a function F(r) is known for a sequence of N + 1 selected values ri of the independent variable r, the process by which one estimates the value of this function for intermediate values is called interpolation. Interpolation methods may be used for comparing model results to observations. They may also be used to convert model variables to a different grid, a procedure known as regridding. Yet another application is for converting between a Lagrangian particle field and an Eulerian grid, as is done in semi-Lagrangian transport schemes (Chapter 7) or for treating nonlinear chemistry in an otherwise Lagrangian modeling framework.
Modeling of Atmospheric Chemistry Page 24