Abstract

Starting from three Eulerian second-order nonlinear advection schemes for semi-staggered Arakawa grids B/E, advection schemes of fourth order of formal accuracy were developed. All three second-order advection schemes control the nonlinear energy cascade in case of nondivergent flow by conserving quadratic quantities. Linearization of all three schemes leads to the same second-order linear advection scheme. The second-order term of the truncation error of the linear advection scheme has a special form so that it can be eliminated by modifying the advected quantity while still preserving consistency. Tests with linear advection of a cone confirm the advantage of the fourth-order scheme. However, if a localized, large amplitude and high wavenumber pattern is present in initial conditions, the clear advantage of the fourth-order scheme disappears.

The new nonlinear fourth-order schemes are quadratic conservative and reduce to the Arakawa Jacobian for advected quantities in case of nondivergent flow. In case of general flow the conservation properties of the new momentum advection schemes impose stricter constraint on the nonlinear cascade than the original second-order schemes. However, for nondivergent flow, the conservation properties of the fourth-order schemes cannot be proven in the same way as those of the original second-order schemes. Therefore, demanding long-term and low-resolution nonlinear tests were carried out in order to investigate how well the fourth-order schemes control the nonlinear energy cascade. All schemes were able to maintain meaningful solutions throughout the test.

Finally, the impact was examined of the fourth-order momentum advection on global medium-range forecasts. The 500-hPa anomaly correlation coefficient obtained using the best performing fourth-order scheme did not show an improvement compared to the tests using its second-order counterpart.

1. Introduction

The advent of massively parallel computers invigorated the interest of atmospheric modelers for explicit time differencing for adjustment terms and Eulerian finite-difference advection schemes. The explicit schemes are local in space and thus avoid global dependencies. On the other hand, in contrast to the semi-Lagrangian schemes, the Eulerian schemes require relatively narrow and constant width halos, which simplify and reduce communications.

Historically, perhaps the dominant two schools of thought concerning general criteria for designing Eulerian numerical schemes have been those insisting on the importance of formal accuracy and relying on implied or explicit artificial diffusion or filtering to handle computational problems, and those emphasizing conservation of important, typically quadratic properties of flow, and preservation of important features of differential operators. The former approach is exemplified by the early monograph by Kreiss and Oliger (1973) and papers by Crowley (1968), Boris and Book (1973), Zalesak (1979), Collela and Woodward (1984), and van Leer (1985). More on the further work along these lines can be found in Thuburn (1997) and Shchepetkin and McWilliams (1998). Being straightforward and intuitive, this approach attracted many modelers. The latter approach was pioneered by Arakawa (1966, 1972) and Arakawa and Lamb (1977, 1981), and adopted by several other modelers (e.g., Sadourny 1975a,b; Janjic 1977, 1984; Tripoli 1992; Jacobson 2001). The approach introduced by Arakawa found followers also in applied mathematics, where it is known as “mimetic” method, for mimicking features of continuous physical systems that are being modeled (e.g., Bochev and Hyman 2006; Hyman and Shashkov 1999; Hyman et al. 2002).

The Arakawa approach addressed the problem of excessive accumulation of energy at small scales due to the nonlinear energy cascade that could lead to nonlinear instability (Arakawa 1966, 1972). To control the nonlinear energy cascade, Arakawa suggested to conserve energy and enstrophy in case of nondivergent flow, and designed a numerical scheme that can do that, the renowned Arakawa Jacobian. Thus, the Arakawa approach is based on replicating fundamental properties of nonlinear flows.

On the other hand, the concept of formal accuracy is restricted to linear discrete equations, and does not take into account the nonlinear nature of geophysical flows. However, one may hope that with increasing formal accuracy, the discrete equations will converge to continuous equations, and therefore will approximate the nonlinear features of the flow in an acceptable way.

In practice it is not always easy to draw a clear line between the two approaches since specific designs for a variety of problems encountered in atmospheric modeling often combine features of both. For example, the first-order moments are often conserved in high-order formal accuracy schemes. Similarly, schemes of higher formal order of accuracy can be constructed that have advanced conservation properties (e.g., Abramopoulos 1991; Arakawa 1966; Rancic 1988; Rancic et al. 2008). Such schemes will be developed here as well but using a different approach, and their performance will be compared to that of second-order conservative schemes in idealized cases and in extended real-data runs. Even though the results obtained in few examples cannot be generalized, we may hope that they will give us some indication on relative importance of higher formal order of accuracy and conservation.

2. Derivation of the fourth-order schemes

We will restrict ourselves to the semistaggered Arakawa B/E grids. For simplicity, consider a regular, square Arakawa B grid in plain geometry, shown in the left panel of Fig. 1 together with coordinate systems (x, y) and (x′, y′) that are rotated by 45° with respect to each other. In the figure, symbol h denotes mass variables, while symbol v denotes horizontal velocity components. Following Janjic (1984), as on the E grid, on the B grid the horizontal velocity components can be defined in terms of velocity potential χ and streamfunction ψ using the simplest centered differencing and averaging operators:

 
formula

Here, s denotes an arbitrary coordinate axis and Δ is a finite difference of the values at two neighboring points in the direction of the coordinate axis s. With these definitions, using the distribution of points carrying the velocity potential χ and streamfunction ψ in the right-hand-side panel of Fig. 1, the velocity components along the coordinate axes on the B grid take the following form:

 
formula

With the definitions (1) and (2), the topology of the grid B′ becomes equivalent to that of the grid B.

Fig. 1.

(left) The Arakawa B grid and (right) its equivalent representation B′ using velocity potential χ and streamfunction ψ.

Fig. 1.

(left) The Arakawa B grid and (right) its equivalent representation B′ using velocity potential χ and streamfunction ψ.

Three nonlinear energy cascade controlling, second-order accuracy B-grid momentum advection schemes will be considered for upgrade. Two of them were originally formulated for the E grid (Janjic 1984). Since the E grid is obtained by rotating the B grid for 45°, the alternate coordinate system (x′, y′) shown in Fig. 1. is the “native” coordinate system for the E grid. Thus, applied to shallow-water equations, the Janjic (1984) schemes can be written in the following form:

 
formula

and

 
formula

Here, h is the height of the free surface, and the mass fluxes appearing in (3) and (4) are

 
formula

The wind tendencies and winds u and υ on the B grid can be obtained from the wind tendencies and winds along the axes x′, y′.

Schemes (3) and (4), among other things, conserve both energy and enstrophy (vorticity squared), but with enstrophy computed in different ways. Scheme (3) conserves enstrophy as computed using the most accurate five-point second-order approximation of Laplacian operating on the streamfunction points of the B′ grid:

 
formula

where repeated subscripts denote repeated application of the finite-difference operator they accompany. On the other hand, scheme (4) conserves enstrophy as computed from the B-grid velocity components. In this case the vorticity expressed in terms of the streamfunction takes the following form:

 
formula

As shown in Janjic (1984), and later in this paper, the small differences in the formulation and properties of the two schemes make a big difference in their nonlinear performance.

The third scheme considered,

 
formula

was proposed by Arakawa (1972). In case of nondivergent flow, this scheme can be obtained as an arithmetic mean of the energy and enstrophy conserving schemes (3) and (4) (Janjic 1984). It does not conserve enstrophy exactly, but conserves other quadratic quantities such as and through which it imposes control over the nonlinear cascade.

Upon linearization, all three second-order scalar advection schemes (3), (4), and (8) take the same form:

 
formula

Here, u0 and υ0 are the components of the constant advection speed v0 along the coordinate axes x and y, superscripts 2x and 2y associated with overbars denote averaging over two grid intervals in the direction of the coordinate axis indicated, and Z is the advected quantity.

Denoting the finite-difference operator on the right-hand side by A, the truncation error of the scheme (9) is

 
formula

Here, d = Δx = Δy is the grid size. Using operator A, (10) can be replaced by

 
formula

where is the simplest five-point approximation of the Laplacian defined by (6). After rearrangement, one obtains

 
formula

and finally,

 
formula

where

 
formula

Thus, because of the specific form of the leading term of the truncation error of the nonlinear cascade controlling second-order schemes (3), (4), and (8), a fourth-order accuracy scheme (13)(14) is obtained only by modifying the argument of the advection operator A. Note that the consistency of the scheme is not violated by this modification. This is very convenient because fourth-order schemes can be introduced without changing the existing code, and with only minimal additional computational effort.

Note that following an analogous procedure, a fourth-order scheme could be derived using in (11) any linear combination of the finite-difference Laplacians (6) and (7), or any other approximation of Laplacian with the order of formal accuracy not lower than second. We chose Laplacian (6) because among the second-order Laplacian approximations it has the smallest truncation error, and generally steepest ascent of eigenvalues with increasing wavenumber. As will be explained later, the latter is an advantage in the context of the generalized conservative nonlinear fourth-order schemes that reduce to the linear form (13)(14).

3. Relative phase speeds of the second- and fourth-order schemes

Let the constant advection speed be

 
formula

where

 
formula
 
formula

are the components of the phase speed along the x and y axes. The lowercase letters k and l denote the components of the wavenumber vector along these axes. Then, the relative phase speed for the second-order scheme has the following form:

 
formula

where

 
formula

In case of the fourth-order linear advection scheme, a correction factor appears due to the modification of the argument of the advection operator, so that the relative phase speed takes the following form:

 
formula

The expressions in (18) and (20) are plotted on the left and right panels of Fig. 2, respectively, in the first quadrant of the admissible wavenumber range. As can be seen from the figures, as expected, the fourth-order scheme significantly improves the relative phase speed.

Fig. 2.

Relative phase speed for the (left) second-order scheme and (right) fourth-order scheme.

Fig. 2.

Relative phase speed for the (left) second-order scheme and (right) fourth-order scheme.

4. Tests of the linear second- and fourth-order schemes

The linear second-order scheme in (9) and the fourth-order scheme in (13)(14) were tested in solid-body rotation experiments. The resolution was 8 km and the integration domain had 101 by 101 grid points. A full rotation of 2π is completed in 48 h. As shown in the top-left panel of Fig. 3, in the first test, the initial conditions for Z were specified in the form of a cone with the basis with the radius of 100 km, and the height of 100 arbitrary units. Initially, the cone was centered at the grid point with indices i = 67 and j = 34. To minimize the errors due to time stepping, the second-order Adams–Bashforth scheme with a very short time step of only 1 s was used for time integration.

Fig. 3.

Linearized schemes. (top) Initial conditions, and solutions after one full rotation using the (bottom left) second-order scheme and (bottom right) fourth-order scheme.

Fig. 3.

Linearized schemes. (top) Initial conditions, and solutions after one full rotation using the (bottom left) second-order scheme and (bottom right) fourth-order scheme.

The results of integrations with the second-order scheme in (9) and the fourth-order scheme in (13)(14) after one full rotation are shown in the bottom-left and -right panels of Fig. 3, respectively. As can be seen from the figure, the adverse impact of computational dispersion is clearly visible in the integration using the second-order scheme. The minimum values are negative and reach 10% of the initial height of the cone. In contrast, in the experiment with the fourth-order scheme, a significant improvement has been achieved, and the minimum values, although still negative, are reduced to a few percents of the initial cone height. The rms error after one full revolution was 2.33 for the second-order scheme, and only 0.005 59 for the fourth-order scheme. Thus, the experiments confirmed the advantages of higher formal order of accuracy in this example.

However, fitting high-order polynomials to noisy data may be inaccurate in between data points, which may result in inaccurate estimates of spatial derivatives in case of schemes of higher formal order of accuracy. More specifically, although the higher-order schemes improve the phase speeds in most part of the admissible wave range, they are not free of the computational dispersion. The relative phase speed of the shortest resolvable waves still remain zero, and therefore, the relative phase speeds of higher-order schemes generally decay faster with increasing wavenumbers in the small-scale range than those of the second-order schemes. This behavior is clearly visible in Fig. 2. However, the steeper decent of the relative phase speed with increasing wavenumber in the case of higher-order schemes means larger group speed errors in the small-scale range.

To visualize this problem, let a wavy pattern of the following form:

 
formula

be superposed on the conical initial condition of the preceding experiment. In (21), q0 is the amplitude of the noise equal to the height of the cone, xc and yc are the coordinates of the center of the base of the cone, and Rc is its radius. Note that the dominant noise wavelength is about six grid intervals. The new initial conditions are shown in the top-left panel of Fig. 4. In atmospheric models such small-scale noisy patterns can be generated (e.g., by forcing by mountains).

Fig. 4.

Linearized schemes. (top) Noisy initial conditions, and solutions after one full rotation using the (bottom left) second-order scheme and (bottom right) fourth-order scheme.

Fig. 4.

Linearized schemes. (top) Noisy initial conditions, and solutions after one full rotation using the (bottom left) second-order scheme and (bottom right) fourth-order scheme.

The solutions after one full rotation using the second-order scheme in (9) and the fourth order scheme in (13)(14) are shown in the bottom-left and -right panels of Fig. 4, respectively. As can be seen from the figure, with both schemes the small-scale part of the initial conditions disperses leaving the large-scale cone behind. The trailing noise in the fourth order scheme is stronger and stays closer to the large-scale cone presumably due to reduced phase speed errors. Both solutions are noisy and visual impression is no longer sufficient to determine which one is closer to the true solution. However, the rms errors for the two schemes are now comparable: 16.18 m for the second-order scheme and 15.52 m for the fourth-order scheme. Such a result could be expected considering that the superposed noisy pattern had spectral components with significant amplitudes in the spectral ranges with substantial dispersion errors with both schemes.

In atmospheric models the situation is further complicated by the nonlinearity of the advective process, and possibly numerical or physical diffusion on coordinate surfaces of terrain-following coordinate systems. So, the formal accuracy of a scheme may not be automatically transferred into the actual accuracy of computations in an atmospheric model.

5. Implementation of schemes of higher formal order of accuracy in atmospheric models

In models using schemes with formal order of accuracy in space higher than second, it has been the usual practice to apply second-order accuracy schemes to the adjustment terms, and to use higher-order accuracy schemes only for advection terms in the momentum and thermodynamic equations. Thus, the overall formal accuracy in space of the algorithms used in such models is still of the second order.

However, there is an inherent difficulty in distinguishing adjustment terms from advection terms in systems that are more complex than the shallow-water equations. Namely, as pointed out early in the history of development of atmospheric models (e.g., Kurihara 1968; Arakawa 1972), but also recently in a somewhat different context (Klemp et al. 2003), in order to conserve energy and provide consistent energy transformations between kinetic and potential energy in hydrostatic limit, it is important to discretize the omega-alpha term of the thermodynamic equation consistently with the discretization of the pressure gradient force in the equation of motion. The affected terms of the nonhydrostatic horizontal momentum advection equation in terrain-following coordinates and the thermodynamic equation are

 
formula
 
formula

In (22) and (23) (Janjic et al. 2001), ɛ is the ratio of vertical acceleration and gravity, and the other symbols used have their usual meaning. Thus, since the pressure gradient force is computed using a second-order scheme, the omega-alpha term on the right-hand side of (23) must be second-order accurate, too. Then, for consistency, the temperature advection term that is not shown in (23) should be second-order accurate as well. Therefore, a higher-order accuracy advection scheme can be used only for momentum advection. Of course, the preceding argument does not apply if higher-order schemes are used for all terms, including the adjustment ones.

6. Properties of nonlinear advection schemes and nonlinear tests

After conversion to the fourth-order accuracy, the schemes (3), (4), and (8), respectively, take the following form:

 
formula
 
formula

and

 
formula

Here, denotes the argument modifying operator in (14). For convenience, the original second-order schemes, (3), (4), and (8) will be called c2, e2, and a2, respectively, while their fourth-order extensions will be denoted c4, e4, and a4.

The resulting nonlinear fourth-order schemes for velocity components are quadratic conservative and reduce to the Arakawa Jacobian form in case of nondivergent flow. In case of general flow the conservation properties of the momentum advection schemes impose stricter constraint on the nonlinear energy cascade than the original energy conserving second-order schemes. As can be easily verified, if summation is performed over a closed domain, or a domain with double periodic boundary conditions, from (24) to (26):

 
formula

where the overbar denotes the summation over the entire domain. Let the velocity components appearing in (27) be expended in finite orthogonal series of the following form:

 
formula

where ξn are normalized orthogonal functions. Then, applying the operator (14) to (28):

 
formula

Here, λn are the eigenvalues of the finite-difference Laplacian . Finally, assuming that h in (27) is constant in order to avoid difficulties with triple products in (27), substituting the expressions for velocity components (28) and (29) in (27),

 
formula

and taking into account an analog of Parseval theorem for finite sums,

 
formula

Thus,

 
formula

As can be seen from (32), the presence of the term associated with the finite-difference Laplacian imposes stricter control on the amplitudes of the wavenumber components than the conservation of energy does. Recall that at the end of section 2 we praised Laplacian (6) for its generally steepest ascent of eigenvalues with increasing wavenumber among the second-order Laplacian approximations. From (32) it becomes clear why this feature is beneficial.

However, the conservational properties of the fourth-order schemes in the case of nondivergent flow cannot be proven in the same way as those of the original second-order schemes. Therefore, nonlinear tests were performed following the linear ones in order to check whether and how well the fourth-order schemes control the nonlinear energy cascade. A difficult nonlinear test introduced in Janjic (1984), and subsequently used also by Rancic (1988), Rancic et al. (2008), and Gavrilov (1997), was carried out. Few nonlinear schemes have passed this test producing meaningful results.

In the test, nonlinear shallow-water equations are solved in a rotating rectangular domain with mirror image boundary conditions and with dissipation reduced to a minimum (“Rotating Flat Square Earth”). The domain was covered with only 17 × 17 grid points. The tests were run for 116 days of simulated time. The leapfrog scheme with weak Asselin filter was used for time integration.

The performance of the schemes with respect to nonlinear energy cascade is assessed monitoring the evolution of the quantity (Janjic 1984):

 
formula

which is normalized by its initial value. The quantity in (33) is not conserved, either by the physical system or the numerical schemes considered. However, being proportional to the fourth power of the eigenvalues of the finite-difference Laplacian, a large increase of this parameter would signal a major downscale energy transfer.

The results of the tests of the second-order schemes c2 (blue), e2 (green), and a2 (red) are shown in the top panel of Fig. 5 (cf. Rancic et al. 2008). The results for the fourth-order schemes c4, e4, and a4 are shown in the bottom panel of Fig. 5 using the same colors. All six schemes performed reasonably well. Among the second-order schemes the best result was obtained with the scheme c2 that conserved enstrophy as computed by the most accurate second-order approximation of Laplacian operating on the streamfunction points. This scheme was closely followed by the Arakawa (1972) scheme a2, while the scheme conserving enstrophy computed from the B-grid velocity components e2 was distant third. The fourth-order schemes c4, e4, and a4 ranked in the same order as the second-order schemes from which they had been derived, and generally were competitive with the second-order schemes in preventing accumulation of energy at the small scales due to the nonlinear cascade throughout the experiments. The scheme e4 actually demonstrated a notable improvement over its second-order counterpart, which may indicate that its nonlinear behavior benefited from enhanced formal accuracy. However, the schemes c4 and a4 did not show any qualitative improvement in this test over the second-order schemes c2 and a2, which were already exceptionally good.

Fig. 5.

Noise measure as a function of time over 116 days in (top) second-order and (bottom) fourth-order very low-resolution nonlinear tests.

Fig. 5.

Noise measure as a function of time over 116 days in (top) second-order and (bottom) fourth-order very low-resolution nonlinear tests.

7. Medium-range forecasting tests

The National Centers for Environmental Prediction (NCEP) Nonhydrostatic Multiscale Model on the B grid (NMMB; Janjic 2005; Janjic and Black 2007) was used for medium-range global forecasting tests. The model was initialized and verified against operational analyses generated by NCEP’s Global Forecasting System (GFS).

The NMMB has been designed for a broad range of spatial and temporal scales and follows the general modeling philosophy of the NCEP’s Weather Research and Forecasting (WRF) NMM gridpoint dynamical core (Janjic et al. 2001; Janjic 2003). However, in contrast to the WRF NMM that uses the Arakawa E grid, the NMMB was reformulated for the Arakawa B grid. The NMMB uses the regular latitude–longitude grid for the global domain, and a rotated latitude–longitude grid in regional applications. The nonhydrostatic component of the model dynamics is introduced through an add-on module so that it can be turned on or off depending on resolution.

In the model, “isotropic” finite-volume horizontal differencing is employed that conserves a variety of basic and derived dynamical and quadratic quantities and preserves some important properties of differential operators. The second-order energy and enstrophy conserving scheme c2 is used for horizontal advection of momentum, but the c4 fourth-order scheme that turned out to be the best in the nonlinear tests of the preceding paragraph is available as an option, too. The conservation of energy and enstrophy by the c2 scheme, and the experimentally demonstrated control over the nonlinear energy cascade by the c4 scheme, improve the accuracy of the nonlinear dynamics of the model on all scales, and renders the model suitable for extended integrations without excessive dissipation.

“Across the pole” polar boundary conditions are specified in the global limit. The polar filter selectively reduces tendencies of the wave components of the basic dynamical variables that would otherwise propagate faster in the zonal direction than the fastest wave propagating in the meridional direction.

In the vertical, the hybrid pressure-sigma coordinate is used (Simmons and Burridge 1981; Eckermann 2009). The forward–backward scheme (e.g., Janjic 1979) is employed for horizontally propagating fast waves, and an implicit scheme is used for vertically propagating sound waves.

Slightly off-centered Adams–Bashforth scheme is applied for nonsplit horizontal advection of the basic dynamical variables and for the Coriolis force (Janjic et al. 2001; Janjic 2003). Instead of being slightly unstable, due to off-centering, the scheme becomes weakly dissipative. Even though the instability of the second-order Adams–Bashforth scheme is very weak, and can be tolerated in practice, we preferred the weakly damping off-centered scheme since strictly speaking unstable schemes do not converge. Note that even though the off-centered scheme is formally of the first order accuracy, the actual magnitude of its truncation error remains close to that of the second-order Adams–Bashforth scheme due to very small off-centering. On the B/E semistaggered grids the advection time step can be only about twice longer than the forward–backward adjustment step because longer time steps can be used for the adjustment terms than on the C grid. So the Adams–Bashforth scheme for advection with the same time step as that used for the adjustment terms costs the same as a two-step iterative scheme with a twice longer time step. On top of that, since there is no time splitting, there is no need for iterating the adjustment terms, and the short time steps reduce the time-stepping errors. So, the noniterative, nonsplit Adams–Bashforth scheme offers significant savings, and at the same time its short time step reduces numerical errors.

To eliminate computational stability problems due to thin vertical layers, the Crank–Nicholson scheme is used to compute the contributions of vertical advection. As a compromise between requirements for computational affordability and accuracy, a fast Eulerian conservative and positive-definite scheme was developed for model tracers (Janjic et al. 2009). Conservative monotonization is applied in order to control oversteepening within the conservative and positive-definite tracer advection scheme. A variety of physical parameterizations have been coupled to the model. The resolution used in the tests was 0.33° in the meridional and 0.47° in the zonal direction. The model had 64 vertical levels. This resolution is comparable to that used operationally for medium-range forecasting at NCEP.

The results of the tests performed with the second-order scheme c2 and the fourth-order scheme c4 are summarized in Fig. 6. In the figure, the anomaly correlation coefficients at 500 hPa calculated over the globe are shown as functions of forecast time for the second-order scheme c2 (blue diamonds) and the fourth order scheme c4 (purple squares). The results shown represent an average of 111 cases run twice daily between 17 November 2009 and 16 January 2010. A few cases are missing because of random computer maintenance outages.

Fig. 6.

Anomaly correlation coefficients calculated over the globe and averaged over 111 cases for the (blue diamonds) second-order scheme c2 and the (purple squares) fourth-order scheme c4.

Fig. 6.

Anomaly correlation coefficients calculated over the globe and averaged over 111 cases for the (blue diamonds) second-order scheme c2 and the (purple squares) fourth-order scheme c4.

As can be seen from the figure, the average scores for the two schemes are very similar, and do not confirm the common expectation that higher-order schemes produce more accurate forecasts. To illustrate how close to each other the scores are, zoomed in difference between the scores for the fourth-order scheme c4 and the second-order scheme c2 is shown in Fig. 7. The difference increases somewhat as the forecasts approach day 9, but by then both scores drop well below 0.6, which is often used as the threshold for usefulness of forecasts.

Fig. 7.

Difference between anomaly correlation coefficients calculated over the globe and averaged over 111 cases for the second-order scheme c2 and the fourth-order scheme c4.

Fig. 7.

Difference between anomaly correlation coefficients calculated over the globe and averaged over 111 cases for the second-order scheme c2 and the fourth-order scheme c4.

A reason for such result may be the presence of small-scale forcing and small-scale features produced by steep and high orography in real-data runs. As demonstrated in section 4, the rms errors of linear second- and fourth-order advection schemes can be similar if spectral components with significant amplitudes are present in spectral ranges with large computational dispersion errors. Also, one should have in mind that the second-order scheme c2 used in the tests is very specific because of its conservation properties and the properties of differential operators it preserves. To compute advection tendencies, scheme c2 uses information from nine grid points. Similarly, straightforward fourth-order schemes also use information from nine points, but from a different gridpoint stencil. Considering that the real flows are nonlinear, favorable features of scheme c2 with respect to the nonlinear energy cascade may outweigh the inaccuracies due to reduced formal accuracy of its linearized version. In other words, it should not be surprising that for simulation of nonlinear flows, criteria other than formal accuracy may be more important.

8. Summary of the main results and conclusions

In this work, three advanced fourth-order advection schemes have been developed. These schemes were derived from nonlinear energy cascade controlling second-order schemes that have been used for years for weather forecasting and climate studies on various scales. Switching from the second- to the fourth-order advection scheme is easy, and requires only a minor increase in computational cost.

The nonlinear fourth-order schemes are quadratic conservative and reduce to the Arakawa Jacobian form in case of nondivergent flow. In case of general flow, the conservation properties of the momentum advection schemes impose stricter constraint on the nonlinear energy cascade than the original energy conserving second-order schemes.

The fourth-order schemes were tested and compared to their second-order counterparts in models of increasing complexity. The linear solid-body rotation tests with advection of a cone confirm the advantage of the fourth-order version of the linear fourth-order scheme both visually and in terms of the rms error. However, already in case of linear advection with the cones in the initial conditions replaced by cones with superposed noisy patterns localized over the cone base area, the clear advantage of the fourth-order scheme in terms of the rms error disappears.

In nonlinear tests the fourth-order schemes ranked in the same order as the second-order schemes from which they were derived, and generally were competitive with the second-order schemes in preventing accumulation of energy at the small scales due to the nonlinear cascade throughout the experiments. Among the second-order schemes the best result was obtained with the scheme that conserved enstrophy as computed using the most accurate second-order approximation of Laplacian operating on the streamfunction points (Janjic 1984). This scheme was closely followed by the Arakawa (1972) scheme, while the scheme conserving enstrophy computed from the B-grid velocity components was distant third.

Furthermore, in order to examine the impact of increased formal accuracy a series of global medium-range forecasts were computed using the second- and fourth-order schemes that performed best in the nonlinear tests. The fourth-order advection scheme was not used for advection of temperature in the thermodynamic equation in order to preserve consistency between the pressure and temperature advections with the second-order scheme used for computation of the pressure gradient force in the momentum equation. Therefore, the fourth-order scheme was applied only to the horizontal momentum.

Contrary to what is generally expected, there was no significant difference between the scores for the forecasts using the second- and fourth-order momentum advection schemes. The average scores for the second- and fourth-order schemes were hardly distinguishable, and thus did not confirm the common expectation that higher formal order of accuracy would result in more accurate forecasts.

A reason for such result may be the presence of small-scale forcing and small-scale features produced by steep and high orography in real-data runs. As demonstrated in section 4, the rms errors of linear second- and fourth-order advection schemes can be similar if spectral components with significant amplitudes are present in spectral ranges with large computational dispersion errors. Also, one should have in mind that the second-order scheme used in the tests is very specific because of its conservation properties and the properties of differential operators it preserves. Note that in order to compute advection tendencies, the second-order scheme uses the same amount of information from neighboring grid points as do straightforward fourth-order schemes, although this information is obtained from points belonging to different gridpoint stencils. Considering that the real flows are nonlinear, favorable features of the scheme with respect to the nonlinear energy cascade may outweigh the inaccuracies due to reduced formal accuracy of its linearized version. In other words, it should not be surprising that for simulations of nonlinear flows, criteria other than formal accuracy may be equally or more important.

Finally, even though the specific example presented here demonstrated that second-order schemes can be as accurate as fourth-order schemes in real-data runs, this still does not mean that the same conclusion will be valid for any two second- and fourth-order schemes (i.e., that the formal order of accuracy is irrelevant). However, this example showed that the formal accuracy is not the only, and may not be the most important criterion in designing finite-difference schemes.

REFERENCES

REFERENCES
Abramopoulos
,
F.
,
1991
:
A new fourth-order enstrophy and energy conserving scheme
.
Mon. Wea. Rev.
,
119
,
128
133
.
Arakawa
,
A.
,
1966
:
Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part 1
.
J. Comput. Phys.
,
1
,
119
143
.
Arakawa
,
A.
,
1972
:
Design of the UCLA general circulation model
.
Tech. Rep. 7, Department of Meteorology, University of California, Los Angeles, Los Angeles, CA, 116 pp
.
Arakawa
,
A.
, and
V. R.
Lamb
,
1977
:
Computational design of the basic dynamical processes of the UCLA general circulation model
.
Methods in Computational Physics, J. Chang, Ed., Vol. 17, Academic Press, 173–265
.
Arakawa
,
A.
, and
V. R.
Lamb
,
1981
:
A potential enstrophy and energy conserving scheme for the shallow water equations
.
Mon. Wea. Rev.
,
109
,
18
36
.
Bochev
,
P. B.
, and
J. M.
Hyman
,
2006
:
Principles of mimetic discretizations of differential operators
.
Compatible Spatial Discretizations, D. N. Arnold et al., Eds., IMA Volumes in Mathematics and its Applications, Vol. 142, Springer, 89–119
.
Boris
,
J. P.
, and
D. L.
Book
,
1973
:
Flux corrected transport. I: SHASTA, a fluid transport algorithm that works
.
J. Comput. Phys.
,
11
,
38
39
.
Collela
,
P.
, and
P.
Woodward
,
1984
:
The piecewise parabolic method (PPM) for gas-dynamical simulations
.
J. Comput. Phys.
,
54
,
174
201
.
Crowley
,
W. P.
,
1968
:
Numerical advection experiments
.
Mon. Wea. Rev.
,
96
,
1
12
.
Eckermann
,
S.
,
2009
:
Hybrid σ–p coordinate choices for a global model
.
Mon. Wea. Rev.
,
137
,
224
245
.
Gavrilov
,
M. B.
,
1997
:
Integration of the shallow water equations in a plane geometry using semi-Lagrangian and Eulerian schemes
.
Meteor. Atmos. Phys.
,
62
,
141
160
.
Hyman
,
J. M.
, and
M.
Shashkov
,
1999
:
The orthogonal decomposition theorems for mimetic finite difference methods
.
SIAM J. Numer. Anal.
,
36
,
788
818
.
Hyman
,
J. M.
,
J. E.
Morel
,
M.
Shashkov
, and
S. L.
Steinberg
,
2002
:
Mimetic finite difference methods for diffusion equations
.
Comput. Geosci.
,
6
,
333
352
.
Jacobson
,
M. Z.
,
2001
:
GATOR-GCMM: A global- through urban-scale air pollution and weather forecast model 1. Model design and treatment of subgrid soil, vegetation, roads, rooftops, water, sea ice, and snow
.
J. Geophys. Res.
,
106
(
D6
),
5385
5401
.
Janjic
,
Z.
,
1977
:
Pressure gradient force and advection scheme used for forecasting with steep and small scale topography
.
Contrib. Atmos. Phys.
,
50
,
186
199
.
Janjic
,
Z.
,
1979
:
Forward-backward scheme modified to prevent two-grid-interval noise and its application in sigma coordinate models
.
Contrib. Atmos. Phys.
,
52
,
69
84
.
Janjic
,
Z.
,
1984
:
Nonlinear advection schemes and energy cascade on semi-staggered grids
.
Mon. Wea. Rev.
,
112
,
1234
1245
.
Janjic
,
Z.
,
2003
:
A nonhydrostatic model based on a new approach
.
Meteor. Atmos. Phys.
,
82
,
271
285
.
Janjic
,
Z.
,
2005
:
A unified model approach from meso to global scales
.
Geophysical Research Abstracts, Vol. 7, Abstract 05582. [Available online at http://meetings.copernicus.org/www.cosis.net/abstracts/EGU05/05582/EGU05-J-05582.pdf.]
Janjic
,
Z.
, and
T.
Black
,
2007
:
An ESMF unified model for a broad range of spatial and temporal scales
.
Geophysical Research Abstracts, Vol. 9, Abstract 05025. [Available online at http://meetings.copernicus.org/www.cosis.net/abstracts/EGU2007/05025/EGU2007-J-05025.pdf.]
Janjic
,
Z.
,
J. P.
Gerrity
Jr.
, and
S.
Nickovic
,
2001
:
An alternative approach to nonhydrostatic modeling
.
Mon. Wea. Rev.
,
129
,
1164
1178
.
Janjic
,
Z.
,
H.
Huang
, and
S.
Lu
,
2009
:
A unified atmospheric model suitable for studying transport of mineral aerosols from meso to global scales
.
IOP Conf. Ser.: Earth Environ. Sci.
,
7
(
1
),
012011
,
doi:10.1088/1755-1307/7/1/012011
.
Klemp
,
J. B.
,
W. C.
Skamarock
, and
O.
Fuhrer
,
2003
:
Numerical consistency of metric terms in terrain-following coordinates
.
Mon. Wea. Rev.
,
131
,
1229
1239
.
Kreiss
,
H.
, and
J.
Oliger
,
1973
:
Methods for the Approximate Solution of Time Dependent Problems
.
GARP Publications Series, Vol. 10, WMO and ICSU, 107 pp
.
Kurihara
,
Y.
,
1968
:
Note on finite difference expressions for the hydrostatic relation and pressure gradient force
.
Mon. Wea. Rev.
,
96
,
654
656
.
Rancic
,
M.
,
1988
:
Fourth-order horizontal advection schemes on the semi-staggered grid
.
Mon. Wea. Rev.
,
116
,
1274
1288
.
Rancic
,
M.
,
H.
Zhang
, and
V.
Savic-Jovcic
,
2008
:
Nonlinear advection schemes on the octagonal grid
.
Mon. Wea. Rev.
,
136
,
4668
4686
.
Sadourny
,
R.
,
1975a
:
The dynamics of finite-difference models of the shallow water equations
.
J. Atmos. Sci.
,
32
,
680
689
.
Sadourny
,
R.
,
1975b
:
Compressible model flows on the sphere
.
J. Atmos. Sci.
,
32
,
2103
2110
.
Shchepetkin
,
A. F.
, and
J. C.
McWilliams
,
1998
:
Quasi-monotone advection schemes based on explicit locally adaptive dissipation
.
Mon. Wea. Rev.
,
126
,
1541
1580
.
Simmons
,
A. J.
, and
D. M.
Burridge
,
1981
:
An energy and angular-momentum conserving vertical finite-difference scheme and hybrid vertical coordinates
.
Mon. Wea. Rev.
,
109
,
758
766
.
Thuburn
,
J.
,
1997
:
TVD schemes, positive schemes, and the universal limiter
.
Mon. Wea. Rev.
,
125
,
1990
1993
.
Tripoli
,
G. J.
,
1992
:
A nonhydrostatic mesoscale model designed to simulate scale interaction
.
Mon. Wea. Rev.
,
120
,
1342
1359
.
van Leer
,
B.
,
1985
:
Upwind-difference methods for aerodynamic problems governed by the Euler equations
.
Lectures in Applied Mathematics, Part 2, C. B. Vreugdenhil and B. Koren, Eds., Vieweg, 117–138
.
Zalesak
,
S. T.
,
1979
:
Fully multidimensional flux corrected transport for fluids
.
J. Comput. Phys.
,
31
,
335
362
.