Bott’s forward-in-time, positive-definite, area-preserving flux-form advection algorithm is extended to higher orders. The new algorithm versions are tested with two-dimensional flow fields: a purely rotational flow and a purely deformational flow. The results suggest that a balance between required accuracy and reasonable computational cost is obtained for schemes with orders of 3–6, with distinction to the fifth-order algorithm. It is also shown that modified odd-order schemes successfully avoid significant phase-speed errors and furnish results that are very close to that obtained with the next-higher even-order scheme but with a lower computational cost.
Numerical modeling of atmospheric transport phenomena often requires the solution of the advection equation, and a large variety of numerical methods have been developed in past decades to solve it. However, many of these schemes produce undesirable results such as large numerical diffusion and damping or negative values for an expected positive-definite quantity.
The development of high-order schemes, such as the various higher-order versions of the upstream scheme (Tremback et al. 1987), allowed solutions with very small numerical diffusion, and the advent of positive-definite advective schemes (e.g., Smolarkiewicz 1983), based on the concept of flux-corrected transport (FCT) method (Boris and Book 1973), represented a significant improvement to the evaluation of the advective transport of positive-definite quantities.
At the end of the 1980s, Bott (1989a) proposed an advection scheme that is area preserving and positive definite, in which a weighting flux treatment was introduced. Subsequently, it was proved that the coefficients of the polynomial fitting of an advected variable Ψ can be directly obtained from the assumption of the area preservation (Bott 1989b). As pointed out by Bott (1989b), the area-preserving flux-form (APF) algorithm aims at getting better results than the previous flux-form algorithms by reducing both amplitude and phase errors. Further improvements have been made to the APF algorithm by Bott (1992), Easter (1993), and Chlond (1994).
This note proposes an extension of the APF scheme to higher orders and shows the results of tests performed with a rotational flow field and a purely deformational, periodically inverted flow field (Seibert and Morariu 1991).
2. Extension of the APF algorithm to higher orders
According to Tremback et al. (1987), the even-order versions of the conventional upstream scheme are more attractive than the odd-order ones and the sixth-order scheme presents the most appropriate balance between accuracy and computational cost.
It can be verified that, in the case of a constant grid spacing and not considering the flux-weighting procedure and positive-definite or monotone constraints, the l-order primary version of Bott’s APF scheme (Bott 1989a) can be reduced to the (l + 1)-order scheme presented in Tremback et al. (1987). Due to the similarity of the schemes, there is an expectation that Bott’s APF method behaves like Tremback’s upstream algorithm. It is natural to try to extend the APF scheme to its higher orders and to verify if the fifth-order Bott’s algorithm offers a better balance between accuracy and numerical effort than the fourth-order scheme, for instance.
In the APF scheme, the advected variable Ψ, within the grid box j, is described by
where x′ = (x − xj)/Δx and l gives the order of the scheme.
To extend it to orders from 5–8, we first solve the following linear system for the coefficients aj,k:
For the even-order cases, m ranges from −l/2 to l/2. For the odd cases, two sets of coefficients can be obtained for −(l − 1)/2 < m < (l + 1)/2 and for −(l + 1)/2 < m < (l − 1)/2. It can be seen that these versions of the odd-order schemes (respectively, “a” and “b”) differ from the choice of one additional point in opposite directions. As was pointed out by Easter (1993), they produce significant, but opposing, phase errors. The coefficients of “c” versions may be obtained by a simple average between “a” and “b” coefficients, and it can be verified that the phase errors produced in this manner are much less significant than those related to “a” and “b” versions.
Table 1 shows the coefficients aj,k obtained by solving the system for l = 5–8, with only the “c” versions of the odd-order schemes.
3. Numerical experiments
The higher-order versions of the APF scheme with positive-definite constraints are subjected to two-dimensional tests involving distribution functions with a zero background. Only the results with the “c” version of the odd schemes are presently analyzed because of their smaller phase errors when compared to “a” and “b” versions. An alternating time-splitting procedure was adopted, with a one-dimensional “x” advection step followed by “y,” then “y” followed by “x.”
a. Rotational flow field—Wide distribution
The purely rotational flow field test carried out by Smolarkiewicz (1982) and Bott (1989a,b) is reproduced. In this test, a cone with maximum height 3.87 and base radius 15Δx, centered at the point (50Δx, 75Δy) of the 100Δx × 100Δx domain, is subjected to a rigid-body rotation with constant angular velocity ω = 0.1. A time step Δt = 0.1 is adopted, which provides a complete revolution at each 628 steps.
Table 2 shows the results for this test, and Fig. 1 contains the analytical and numerical solution for six revolutions with l = 2, 3c, and 5c versions of the APF algorithm (just the region of the domain 25Δx ≤ x ≤ 75Δx, 50Δy ≤ y ≤ 100Δy is represented). Numerical diffusion is reduced when the order of the advection scheme increases. A significant improvement in the numerical solution is already obtained when l increases from 2 to 3c. For orders greater than 6, it is clear that the growth of the computational cost becomes more important than the reduction of the numerical diffusion.
As can be seen, good results are already achieved with the l = 3c scheme. However, more severe tests with narrower distributions demonstrate the importance of using l ≥ 5 schemes.
b. Rotational flow field—Narrow distribution
A rotational flow field test, quite similar to that proposed by Chock (1991), was performed. A conical hill with base radius 4Δx and a peak of 100 is centered at the point (7Δx, 17Δy) of a 32Δx × 32Δy domain. It is assumed that a time of 7200;t2π is required to complete one revolution. The maximum Courant number in this simulation is 0.42.
Table 3 and Fig. 2 show clearly that the use of higher-order versions of the APF scheme is strongly recommended in the case of narrow distributions, with a significant reduction of the effects of the numerical diffusion, even in the case of long-term simulations. In Fig. 3, we show the evolution of the peak value of the numerical solution for several orders of the APF scheme.
In this case, to achieve good results, the l ≥ 5 schemes should be used and odd-order schemes are recommended, because their performance is nearly equal to that of the next-higher even-order schemes. The comparison between the increase in peak value (Table 3) and the increase in CPU time (Table 2) shows that a significant improvement in the numerical solution is obtained when the order of the scheme passes from 3c to 5c. If the order of the scheme increases more (from 5c to 7c, for instance), the computational cost becomes more important, even if 20 revolutions are considered. In simulations with a very large number of time steps, the use of the l = 7c scheme is certainly profitable, but l = 5c provides a better balance between efficiency and accuracy in the present case.
c. Deformational flow field
A purely deformational flow field, similar to that proposed by Seibert and Morariu (1991), is adopted for the present tests. One advantage in using this flow field test instead of using the multiple vortices deformational flow field is that the comparison between numerical and analytical solutions become very simple.
The components of the velocity are given by u = αx; v = αy, with a value of α = 0.1.
A pyramid (Ψmax = 1.0; base radius 4Δx) is placed in the center of a 32Δx × 32Δy grid. When the total time of integration reaches t = ln2, the flow is inverted. New flow inversions are made at each 2 ln2 units of time.
The analytical solution for this test for t = 2n ln2 (n integer) is identical to the initial condition. A time step Δt = 0.069 is adopted; the maximum one-dimensional Courant number in this test is 0.111.
Table 4 and Fig. 4 show numerical results for t = 80 ln2 (800 iterations). As in the previous simulations, significant improvements with a reasonable computational cost can be achieved for l = 5c. It should be noted that the overall performance of the APF scheme became worse in the deformational case, compared with the rotational case.
The APF algorithm is extended to higher orders (5–8) and the explicit formulas for the coefficients are presented. Only the positive-definite version has been used, but all the present development can be also applied to the monotone version.
Numerical experiments suggest that a good balance between accuracy and computational cost can be achieved with some of the schemes from l = 3 to 5, depending on the case. In the rotating wide cone test, an accurate result is already obtained with l = 3. However, the rotating narrow cone test and the deforming pyramid test show clearly that the order of the scheme should be increased at least to l = 5. Higher-order schemes (seventh, eighth) exhibit some improvements to the numerical solution, but it seems that they do not compensate for the growth of the CPU time. Also, the increase in the number of neighboring points required to evaluate Ψ makes these schemes disadvantageous, except in cases with narrow distributions and a very large number of time iterations.
It was also clearly verified that the use of modified odd-order schemes (Easter 1993) is more advantageous. The odd-order schemes and the next-higher even-order schemes lead to nearly equal amplification factors but present different phase speed errors, which are bigger in the odd-order case, at least in its original form. In contrast with these original odd-order schemes, the modified odd-l-order algorithms, due to the reduction of phase errors, offer solutions quite similar to those obtained with the even-(l + 1)-order, with a considerably smaller computational cost.
Therefore, the (l = 5c) scheme seems to have the most appropriate combination of accuracy and efficiency even in the cases of narrow distributions or deformational fields. Our results agree with Tremback et al.’s (1987) statement about the superiority of their sixth-order version of the forward-in-time advection algorithm and lead us to conclude that the l = 5c APF scheme is strongly recommended for many atmospheric applications.
The authors wish to thank Dr. Murilo Pereira de Almeida and Dr. José Carlos Parente de Oliveira of the Cloud Physics and Mesoscale Laboratory (Federal University of Ceará) for their comments.
Corresponding author address: Dr. Alexandre A. Costa, Dept. of Atmospheric Science, Colorado State University, Fort Collins, CO 80523.