In finite-difference representations of the conservation equations, the flux form of the advection terms is often preferred to the advective form because of the immediate conservation of advected quantity. The scheme can be designed to further conserve higher-order moments, for example, the kinetic energy, which is important to the suppression of nonlinear instability. It is pointed out here that in most cases an advective form that is numerically equivalent to the flux form can be found, for schemes based on centered difference. This is also true for higher-order schemes and is not restricted to a particular grid type. An advection scheme that is fourth-order accurate in space for uniform advective flows is proposed that conserves both first and second moments of the advected variable in a nonhydrostatic framework. The role of the elastic correction term in addition to the pure flux term in compressible nonhydrostatic models is also discussed.
In finite-difference representations of the conservation equations (e.g., for momentum, mass, heat, and other substances), the flux form of advection terms is more commonly used in primitive equation models because of the immediate conservation of the advected quantities (e.g., momentum, mass, and water substances). The scheme can be designed to further conserve higher-order moments (here the moments are multiplied by the density in a nonhydrostatic system and by water depth in the case of shallow water equations), for example, the kinetic energy, which is important to the control of nonlinear instability. Second-order discretization schemes that conserve both the first (e.g., momentum) and the second moment (e.g., kinetic energy) of the advected quantities can be easily formulated in flux form (e.g., Lilly 1964). In general, numerical schemes that conserve more quantities that are conserved by the continuous counterpart perform better, especially in maintaining correct energy cascade in the discrete system. Arakawa (1966) and Arakawa and Lamb (1977) showed for hydrostatic systems that the conservation of both kinetic energy and enstrophy (vorticity squared) by advection schemes prevent systematic and unrealistic energy cascade toward high wavenumbers, a cause of nonlinear instability (Phillips 1959). Their pioneering work has since been extended and generalized by many researchers, examples include Sadourny (1972), Arakawa and Lamb (1981), Mesinger (1981), and Janjić (1984). These schemes are all presented in hydrostatic frameworks in which the flows are quasi-two-dimensional, however. Numerical schemes that conserve enstrophy are much harder to design in three-dimensional nonhydrostatic systems. Still, the conservation of kinetic energy and the second moments of scalar variables is also highly desirable in nonhydrostatic systems.
For the above reasons, most contemporary nonhydrostatic numerical models use conservative discretization schemes based on the flux-form advection formulation (e.g., Clark 1977; Anthes and Warner 1978; Xue and Thorpe 1991). The advective form is nonetheless also used in some models (e.g., Klemp and Wilhelmson 1978, referred to as KW78 hereafter) partly due to its simplicity when implementing higher-order (e.g., fourth order) schemes. It should be pointed out that it is certainly possible to design conservative discretization schemes based on the advective form. Early examples can be found in Grammeltvedt (1969), Janjić (1984), and Rančić (1988), among others. Most of these schemes were designed for quasi-two-dimensional hydrostatic flows.
Without special care, discretization schemes based on the advective form of advection terms usually do not conserve even the first moment, for example, the momentum. This is the primary reason why many modelers abandon the use of advective form (e.g., Wilhelmson and Chen 1982, WC82 hereafter; Tripoli and Cotton 1982, TC82 hereafter). Elastic (compressible) meteorological models usually do not predict the full air density using an exact mass conservation equation (KW78;TC82). To cast the advection term in an elastic system in flux form, a correction term due to compressibility (to be referred to as the elastic correction term hereafter) is present (or neglected by some modelers) in addition to the pure flux term (e.g., WC82; Dudhia 1993). It is pointed out in this note that there exists, in most cases, a numerically equivalent advective form to a flux form formulation, and neglecting the elastic correction term does not necessarily improve the conservation of the second moment. One of such equivalent advective forms is used by Skamarock and Klemp (1993), but its exact equivalence to a flux form is not well recognized, at least in the nonhydrostatic modeling community. Such equivalence appears to be implicitly recognized by Janjić (1984) and Rančić (1988), for shallow water equations formulated on Arakawa E grid.
2. The numerical equivalence
Without loss of generality, we consider only 2D equations for u and w and an equation for the scalar potential temperature θ in a nonhydrostatic elastic system, and concentrate on the advective part of the equations, which are written as
where the subscripts x and z denote partial differentiations. Here, fu, fw, and fθ represent the right-hand side (rhs) of the u, w, and θ equations, respectively, and their detailed form varies with applications. The other variables have their usual meanings.
The advection terms in Eqs. (1) are in advective form. In models that do not explicitly predict full air density, the above equations are often multiplied by a base-state density ρ(z) and recast into a flux form:
where u* ≡ ρu and w* ≡ ρw. Systems similar to Eqs. (2) are widely used in elastic cloud models (e.g., TC82). Note that in an anelastic system, the continuity equation is
and the terms enclosed in braces (curly brackets) in Eqs. (2) vanish. In the rest of this paper, we should refer to these terms as the elastic correction terms.
Adopting the following notations,
a second-order centered-difference form of Eqs. (2) on the Arakawa C grid can be written as follows:
It is well known that for an anelastic flow the advective process in the above system conserves momentum ρu, ρw, the kinetic energy ρ(u2 + w2)/2, and ρθ as well as ρθ2 (Lilly 1964). In a compressible system, the density (ρ) weighted divergence Div2 ≡ (δxu* + δzw*) is, though small in general, nonzero, and exact conservations no longer exist. WC82 and Ikawa (1988) choose to neglect the elastic correction terms in an attempt to preserve conservations. Doing so does make the first moment conservative. However, without Eq. (3) being satisfied, the second moment is not conservative either (see the appendix). In fact, neglecting them makes the model integration less stable in some cases (Ikawa 1988; Durran and Klemp 1983). Sound wave damping of a certain type alleviates its impact (Ikawa 1988; Skamarock and Klemp 1993). In a recent effort to compare the results of the Advanced Regional Prediction System (ARPS) model (Xue et al. 1995) against that of the KW78 model with WC82 formulation (which uses flux-form advection but neglects the elastic correction terms), Richardson (1999) found that a good agreement in the time evolution of a multicell convective storm simulation could not be obtained until the neglected elastic correction term is added back to the equation of water vapor mixing ratio in the KW78 model. For the comparison purpose, the ARPS code was modified to the maximum extent possible to match the WC92 version of the KW78 model, except for the treatment of advection. Certain simple bubble tests also show that the solutions including the elastic correction term are better behaved, in terms of producing more coherent structure of the rising bubble. Ikawa (1988) found that his mountain flow solution becomes unstable when the horizontally explicit and vertically implicit mode-split time integration scheme (as used in ARPS and WC82) is used together with the flux-form advection without elastic correction, unless extra effort is made to damp deep horizontal acoustic modes. These results suggest strongly that the elastic correction terms should not be neglected.
Finite differencing can be based on the advective form of Eqs. (2), which are then discretized as follows:
Similar finite differencing is done in Skamarock and Klemp (1993). However, it is not generally recognized that Eqs. (6) are exactly equivalent to Eqs. (5). It takes only some algebraic manipulations to prove for the u equation, for instance, that
Equations (5a) and (6a) are therefore exactly the same. Hence, there is no need to split the advection term into two parts to take advantage of the known conservation properties of flux form. From a practical point of view, formulating the term in advective form cuts the computations by half, and in addition, roundoff errors are also reduced simply because there are fewer floating-point operations. In an anelastic system where Eq. (3) is exactly satisfied, the equivalence between the two systems also holds. The same is true for other systems whose flux formulation of advection terms does not involve a correction term.
Numerical experiments were performed to verify the above observations. We have formulated a 3D compressible model using the split-explicit time integration scheme of KW78. Two runs with formulations (5) and (6) are found, without any explicit diffusion, to be stable for long time integrations, and the results are essentially identical (with the only difference being caused by machine roundoff errors). In contrast, when the advection terms were formulated, for example, for the u equation, in a straightforward fashion on a C grid as
explicit numerical diffusion must be added to control nonlinear instability. Without diffusion, the model run blows up after a certain time.
3. Extension to higher-order and other grid structures
The equivalent advective form can also be found for fourth-order differencing schemes as well as for schemes based on other types of grid structure. For example, for the u equation, a fourth-order accurate difference scheme on the Arakawa C grid can be formulated by a straightforward extension of the conservative second-order scheme:
and for scalar θ,
where the fourth-order (base state) density weighted divergence Div4 is defined at the scalar point as
Terms on the left-hand side of Eqs. (9) and (10) are in flux form whereas those on the rhs are in the equivalent advective form. In an anelastic system where Div4 in (11) is identically zero, the fourth-order advections in both forms conserve first and second moments (see the appendix for proof) of the advected quantities. Again neglecting the Div4 term on the left-hand side to yield a pure flux term, as is done in WC82, does not guarantee the conservation of second moments in a compressible model.
It should also be noted that the formulations in Eqs. (9) and (10) are only quasi-fourth-order accurate, mainly because the linear spatial average is introducing second-order truncation errors that do not cancel each other. High-order spatial averaging must be used to construct a truly fourth-order scheme, and more grid points will be involved in such a scheme. The fourth-order scheme used by WC82 is not truly fourth order for nonconstant advective flow either. Both schemes reduce to the same truly fourth-order scheme when the flow is uniform in space. The scheme in flux form as given in Eqs. (9) and (10) has been implemented in the model of Xue and Thorpe (1991), and the Advanced Regional Prediction System (Xue et al. 1995), and the experiments with a cold downburst/density current flow showed superior solutions to those obtained using the second-order counterparts as given in Eqs. (5) and (6).
Equivalent advection and flux forms based on centered difference can also be found for other types of grid structure, such as the Arakawa B grid used in the nonhydrostatic version of the Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (Dudhia 1993). The key step is to perform proper spatial averages on the advection terms as exemplified in Eqs. (6), (9), and (10), instead of directly evaluating them at the point of the predicted variable [as in Eq. (8)]. Such spatial averaging has the desired dealiasing effect through the conservation of first and second moments. Since it is not the purpose of this note to exhaust all possibilities, we will not go into more detail here.
In this note, we have pointed out that the flux form of the advection terms in conservation equations can, for most schemes based on centered finite difference, have a numerically equivalent formulation in advective form. This equivalent advective form has the same conservation properties as the corresponding flux form. Such equivalence can also be found for higher-order (≥4) schemes and on grid structures other than the Arakawa C type. These schemes are (nearly) conservative up to the second moments in the absence of (with weak) compressibility. The advective form is generally simpler than the flux-form counterpart. In certain compressible atmospheric models, in which an elastic correction term is needed in addition to the pure flux term that includes a base-state density, using the equivalent advective form can save up to half of the computational time. It is also pointed out that neglecting the elastic correction term for the sake of momentum conservation, as is done in at least one nonhydrostatic atmospheric model, is undesirable.
This work was supported by NSF through Grant ATM-8809862 to the Center for Analysis and Prediction of Storms (CAPS), University of Oklahoma, and by Grant NSF ATM-9909007 to the first author. The authors wish to thank Dr. A. Shapiro for reading the original manuscript. Most of the work described in this paper was conducted and a draft manuscript completed in the early 1990s when both authors were working at CAPS. Two anonymous reviewers are thanked for their constructive comments.
Formulating the advection terms in the fourth-order flux form as given in Eq. (9), the u equation becomes
Multiplying (A1) by u, rearranging the terms, and summing the equation over the entire model domain yields
In (A2), the terms inside the first summation on the rhs of Eq. (A2) are of the form (Ai+1/2 − Ai−1/2) or (Ai+1 − Ai−1); thus, all terms except those on the boundary cancel, generating no net energy in the interior domain. The second summation term on the rhs originates from the pure flux term in Eq. (9) and the third summation term comes from the divergence term in the same equation. Both terms represent energy sources when Div4 is not zero. As pointed out previously in the text, even if the Div4 term in Eqs. (9) or (10) [hence the third summation on the rhs of Eq. (A2)] is neglected, the conservation of second moments is not guaranteed. On the other hand when Div4 is zero (or very small), the above formulation possesses an exact (or a good) conservation (of the second moment). Finally, the conservation properties of w and scalar equations can be considered in a similar fashion. When the second moments of all three velocity components are conserved, so is the kinetic energy. The conservation of the first moment by Eq. (A1) is evident when Div4 = 0.
Corresponding author address: Dr. Ming Xue, School of Meteorology, University of Oklahoma, 100 East Boyd, Norman, OK 73019.