The Robert–Asselin time filter is widely used in numerical models of weather and climate. It successfully suppresses the spurious computational mode associated with the leapfrog time-stepping scheme. Unfortunately, it also weakly suppresses the physical mode and severely degrades the numerical accuracy. These two concomitant problems are shown to occur because the filter does not conserve the mean state, averaged over the three time slices on which it operates. The author proposes a simple modification to the Robert–Asselin filter, which does conserve the three-time-level mean state. When used in conjunction with the leapfrog scheme, the modification vastly reduces the impacts on the physical mode and increases the numerical accuracy for amplitude errors by two orders, yielding third-order accuracy. The modified filter could easily be incorporated into existing general circulation models of the atmosphere and ocean. In principle, it should deliver more faithful simulations at almost no additional computational expense. Alternatively, it may permit the use of longer time steps with no loss of accuracy, reducing the computational expense of a given simulation.
From a functional perspective, the task of predicting future weather and climate may be reduced to the following iterative procedure. First, given the state of the atmosphere, ocean, and other Earth-system components at any time (the input), use the governing equations to compute the state at a slightly later time (the output). Then, repeat the loop as many times as required, always using the previous output as the next input.
The above prediction framework presents three main challenges, each of which potentially degrades the reliability of the forecast. First, Earth observations, which always contain measurement errors, are required to serve as the initial state. Second, the vast array of active physical processes and interactions is incompletely known and imperfectly represented in the spatially truncated governing equations. Third, the discrete stepping from one time level to the next is merely an approximation to the exact time-continuous evolution. This paper presents a possible avenue for progress with the third of these three challenges, which has received scant attention compared to the extensive research efforts devoted to the first two.
Pfeffer et al. (1992) have assessed the sensitivity of an atmospheric general circulation model to the time stepping, for fixed spatial discretization and physical parameterizations. They find that two different methods of time discretization—the leapfrog and Matsuno schemes—result in significant quantitative differences in the simulated climate. For example, the leapfrog scheme gives much more precipitation over the western tropical Pacific Ocean and less precipitation over the western North Atlantic Ocean. Therefore, time stepping appears to be an important contributor to model error.
Many different time-stepping methods have been proposed, including the leapfrog scheme (or centered difference scheme), the Matsuno scheme (e.g., Pfeffer et al. 1992), the Adams–Bashforth family of schemes (e.g., Durran 1991) and the Runge–Kutta family of schemes (e.g., Kar 2006). The leapfrog scheme has emerged as the method of choice in weather and climate models, despite related disciplines choosing differently (e.g., Runge–Kutta schemes are widely used in computational fluid dynamics but hardly ever used in numerical weather prediction) and despite evidence that other methods may be superior (e.g., the third-order Adams–Bashforth scheme is more accurate; Durran 1991). The leapfrog scheme is used so widely in weather and climate models probably because it is easy to implement, computationally inexpensive, and has low run-time storage requirements. Indeed, today’s widespread use of the leapfrog scheme in general circulation models is perhaps merely a legacy of computer memory having been such a severe constraint when the models were first developed.
A major problem with the leapfrog scheme is that it admits spurious computational modes (e.g., Mesinger and Arakawa 1976; Haltiner and Williams 1980; Durran 1999). In general, a differential equation that is first order in time has one degree of freedom, but an n-time-level numerical approximation to it constitutes an (n − 1)th-order difference equation with n − 1 degrees of freedom. Of these n − 1 modes, one is the physical mode and the remaining n − 2 are computational modes. The leapfrog is a three-time-level scheme, so one computational mode arises in it, in addition to the physical mode, because a second-order difference equation is used to approximate a first-order differential equation. The computational mode (or parasitic mode) is manifest as a spurious oscillation between even and odd time steps, which is referred to as time splitting.
One possible solution to time splitting is to periodically reinitialize the leapfrog scheme by applying a single step of a two-time-level scheme, which does not admit any computational modes. For example, Pfeffer et al. (1992) apply a single Matsuno step after every 11 leapfrog steps. This approach does not remove the computational mode, but merely resets its amplitude to zero periodically so that it never becomes large enough to be problematic.
The far more widely used solution to time splitting is to apply a time filter during the time-stepping procedure. Robert (1966) designed such a filter for the leapfrog scheme and Asselin (1972) showed that it selectively suppresses the computational mode but leaves the physical mode relatively undamped at low frequencies. The filter is now referred to as the Robert filter, the Asselin filter, or the Robert–Asselin filter. The behavior of the filter has been investigated not only for simple equation sets, with no space dependence, but also for the shallow-water equations (Schlesinger et al. 1983) and the hydrostatic primitive equations (Cordero and Staniforth 2004).
As testament to the filter’s success, Asselin (1972),1 has been cited over 450 times according to one citation database, mostly in journals of meteorology and the atmospheric sciences (around 300 citations) but also in journals of oceanography (around 100 citations) and fluid mechanics (around 50 citations). Examples include the use of the filter in models of regional climate (e.g., Caya and Laprise 1999), palaeoclimate (e.g., Fraedrich et al. 2005), ocean circulation (e.g., Griffies et al. 2001), geophysical fluid dynamics (e.g., Ford 1994; Bartello 2002), rotating laboratory fluids (e.g., Williams et al. 2009), and the atmosphere of Mars (e.g., Hartogh et al. 2005). André Robert’s contributions to numerical modeling, including the time filter, have been reviewed by Staniforth (1997) following a memorial symposium held at the University of Québec in 1994.
Currently, the Robert–Asselin filter is used in
operational numerical weather prediction models, including the Mesoscale Model (MSM) of the Japan Meteorological Agency (JMA), the global model of the Australian Bureau of Meteorology (BOM) and Research Centre (BMRC), the global model (GME) and regional model (COSMO-EU) of Deutscher Wetterdienst (DWD), and the Royal Netherlands Meteorological Institute (KNMI) model;
atmospheric general circulation models for climate simulation, including the ECHAM5 model of the Max-Planck-Institut für Meteorologie (MPI-M) and the Community Atmosphere Model (CAM) of the National Center for Atmospheric Research (NCAR);
ocean general circulation models, including Océan Parallélisé (OPA), the Nucleus for European Modeling of the Ocean (NEMO), the oceans of Met Office Hadley Centre climate models [i.e., the Hadley Centre Coupled Climate Model version 3 (HadCM3), the Hadley Centre Global Environmental Model (HadGEM), the High Resolution Global Environmental Model (HiGEM), and the Fast Met Office U.K. Universities Simulator (FAMOUS)], the Hybrid Coordinate Ocean Model (HYCOM), and (as an option) the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model (MOM); and
models of the fluids in rotating annulus laboratory experiments, including the Quasi-Geostrophic Model for Investigating Rotating fluids Experiments (QUAGMIRE) and the Met Office/Oxford Rotating Annulus Laboratory Simulation (MORALS).
Despite its unquestioned success, the Robert–Asselin-filtered leapfrog scheme suffers from two related problems. First, in addition to suppressing the computational mode, the scheme also weakly suppresses the physical mode. Therefore, physical quantities (e.g., energy) that are conserved by the time-continuous equations are not necessarily conserved by the time-discretized equations when the filter is activated. The damping and nonconservation may be benign for sufficiently short integrations, but possibly not for longer ones.
Second, the Robert–Asselin filter severely degrades the leapfrog scheme’s numerical accuracy, measured as the rate at which the error tends to zero as the time step is progressively refined. Specifically, a numerical scheme is defined to be nth-order accurate if, after a given time interval, the difference between the numerical solution of the time-discretized equations with time step Δt, and the exact solution of the time-continuous equations, scales as (Δt)n as Δt → 0. Higher-order schemes are generally preferred to lower-order schemes, because they may permit the use of longer time steps with no loss of accuracy, reducing the computational expense of a given simulation. The Robert–Asselin-filtered leapfrog scheme is only first-order accurate for amplitude errors (although higher-order contributions may dominate for very small values of the filter parameter).
Because the Robert–Asselin filter is used so widely, simple-to-implement modifications that deliver more faithful simulations are very attractive. The author proposes such a modification in this paper. When used in conjunction with the leapfrog scheme, the modification vastly reduces the impacts on the physical mode and increases the numerical accuracy for amplitude errors by two orders, yielding third-order accuracy. Section 2 motivates the modified filter from a geometrical perspective. Section 3 derives analytically the amplification factor and numerical accuracy for the modified filtered leapfrog scheme, and compares them with the corresponding results for the standard Robert–Asselin-filtered leapfrog scheme. Section 4 concludes the paper with a summary and discussion.
2. The Robert–Asselin filter and proposed modification
The standard Robert–Asselin filter, and the modified filter proposed in this paper, are illustrated graphically in Fig. 1. Suppose that the values of a dependent variable, x, are given at three successive and equally spaced times, tn−1, tn, and tn+1. Then, from a geometrical perspective, the standard filter (Fig. 1a) operates by moving the inner point, with coordinates [tn, xn], a fraction ν toward the midpoint, [tn, (xn−1 + xn+1)/2], of the two outer points. Therefore, the displacement of the inner point under the influence of the filter is
The filter parameter, ν, is usually chosen to be O(0.01–0.2).
Two relevant properties of the three points are their mean,
By displacing xn through the amount d, the standard Robert–Asselin filter reduces the magnitude of the curvature of the three points, |Cn|. When used in conjunction with the leapfrog scheme, this feature of the filter strongly suppresses the computational mode, as desired. But, crucially, the application of the filter does not conserve the three-time-level mean, Mn. The theoretical analysis of section 3 will show that, when used in conjunction with the leapfrog scheme, it is this feature of the filter that severely degrades the numerical accuracy.
In an attempt to include the possibility of conserving the three-time-level mean, the modified filter proposed in this paper (Fig. 1b) acts on the right outer point as well as the inner point. Specifically, for any α satisfying 0 ≤ α ≤ 1, the modified filter displaces xn through the amount αd and xn+1 through the amount (α − 1)d, where d is given by (1). All members of this family of modified filters reduce the magnitude of the curvature of the three points, |Cn|, with α controlling the relative contributions to the reduction from the displacements of the inner and right outer points. The special case α = 1 yields the standard Robert–Asselin filter discussed above, which displaces the inner point only. The special case α = 0 displaces the right outer point only. The special case α = 1/2 will be of particular interest in this paper, because it displaces the inner and right outer points equally and oppositely, conserving the three-time-level mean, Mn.
Before embarking upon the theoretical analysis, we briefly demonstrate the improvement that may be achieved by the proposed modification, when used in conjunction with the leapfrog scheme. We numerically integrate the equations of simple harmonic motion,
by alternately applying a leapfrog step and the modified filter with either (in two separate integrations) α = 1 (i.e., the standard Robert–Asselin filter) or α = 1/2. The numerical solutions so obtained are compared with each other, and with the exact solution, in Fig. 2, at the parameter values given in the caption. Amplitude errors are clearly much smaller with the modified filter than with the standard filter. One consequence is that X2 + Y2, which is conserved by the continuous equations and corresponds to the energy of the oscillation, decreases by 89% using the standard filter, but is approximately conserved using the modified filter, between the beginning and end of the integration shown in the figure.
3. Theoretical analysis
Amplitude and phase errors of time-stepping schemes are traditionally examined by analyzing solutions to the oscillation equation (e.g., Durran 1999), which, for the complex variable F(t), is
where i = −1 and ω is a given (real) angular frequency. Equation (6) is related to (4) and (5) by F = X + iY. Using the modified filter proposed in section 2 to control the computational mode, the leapfrog scheme for (6), with time step Δt, is
In this three-stage method, (7) implements the basic leapfrog scheme and (8) and (9) implement the modified filter, with 0 ≤ α ≤ 1. Here F denotes a provisional value, obtained by applying (7) during the current time step; F denotes another (singly filtered) provisional value, obtained by applying (9) during the current time step; and denotes the definitive (doubly filtered) value, obtained by applying (8) during the next time step. The occurrence of filtered values on the right sides of (7)–(9) makes the scheme recursive: F is overwritten with F as soon as it is calculated, and so is F with .
from which the numerical amplification factor is found to be
Assuming 1 − ν/2 > 0, and taking the output of the square root operator to be the branch with nonnegative real part, then the positive sign (A+) corresponds to the physical mode and the negative sign (A−) to the computational mode. For the special case α = 1, (12) reduces to the amplification factor derived by Asselin (1972) for the standard Robert–Asselin-filtered leapfrog scheme, as expected. The exact solution to (6) is F(t) = F(0) exp(iωt), from which the exact amplification factor is found to be
for comparison with (12).
Figure 3 compares the amplification factor for the numerical solution, in the three cases α = 0, α = ½, and α = 1, with the amplification factor for the exact solution. The exact amplification factor (Fig. 3a) lies on the unit circle in the first quadrant and rotates anticlockwise as ωΔt increases from 0 to 1. The numerical amplification factors for the physical mode (Figs. 3b–d) also rotate anticlockwise in the first quadrant, but depart slightly from the unit circle as ωΔt increases from 0. The growing radii for α = 0 (Fig. 3b) and α = ½ (Fig. 3c), and the shrinking radius for α = 1 (Fig. 3d), correspond respectively to an artificial amplification and suppression of the physical mode. The numerical amplification factors for the computational mode (Figs. 3b–d) rotate clockwise in the second quadrant as ωΔt increases from 0 to 1. They each remain inside the unit circle, corresponding to a suppression of the computational mode, as desired.
As suggested by Fig. 3, the standard Robert–Asselin filter (α = 1) behaves qualitatively differently from all other filters in the modified family (α ≠ 1). Only the standard filter (Fig. 3d), exhibits a point in the complex plane at which the amplification factors for the physical and computational modes meet. The singularity occurs because α = 1 is the only case for which the imaginary term within the square root of (12) vanishes, allowing A+ = A− at ωΔt = 1 − ν/2. For all other values of α, the presence of the imaginary term ensures that there is no value of ωΔt for which A+ = A−, and the singular behavior of the standard Robert–Asselin filter is avoided.
Figure 4 shows in more detail how the magnitudes of the amplification factors depend upon ωΔt. The qualitatively different behavior between the cases α ≠ 1 and α = 1 is clearly visible. The singularity for the case α = 1 (Fig. 4c), which renders the curves nondifferentiable at ωΔt = 1 − ν/2, is replaced for the cases α = 0 (Fig. 4a) and α = ½ (Fig. 4b) with a smooth transition from the small ωΔt regime to the large ωΔt regime. The consequence for the physical mode (A+) is that, for the case α = 1 only, |A+| − 1 changes sign as ωΔt increases from 0 to 1, corresponding to a transition from artificial suppression to artificial amplification. The consequence for the computational mode (A−) is that, for the case α = 1 compared to the cases α ≠ 1, the suppression is much less uniform as ωΔt increases from 0 to 1.
Figure 5 shows an enlarged view of how, as ωΔt → 0, the magnitudes of the numerical amplification factors for the physical mode approach the magnitude of the exact amplification factor. The limiting value, unity, appears to be approached much more rapidly for the case α = ½ than for the cases α = 0 and α = 1, suggesting a higher numerical accuracy. To confirm this suggestion, we Taylor expand the square root in (12) to obtain, for the physical mode,
The coefficients, cn(α, ν), of the terms in the series expansion are real. The first few are given by
For the special case α = 1, which corresponds to the standard Robert–Asselin filter, the coefficient of the quadratic term in (20) is negative, yielding |A+| − 1 ∼ −(ωΔt)2 as ωΔt → 0, in agreement with the leading-order behavior shown in Fig. 5. Hence, the amplitude error per time step varies as (Δt)2 and the amplitude error per unit time varies as Δt. Therefore, the numerical scheme with α = 1 is first-order accurate for amplitude errors and unconditionally stable. Alternatively, for the special case α = 0, the coefficient of the quadratic term in (20) is positive, yielding |A+| − 1 ∼ +(ωΔt)2 as ωΔt → 0, also in agreement with the leading-order behavior shown in Fig. 5. Therefore, the numerical scheme with α = 0 is first-order accurate for amplitude errors and unconditionally unstable.
In contrast, for the special case α = ½, the coefficient of the quadratic term in (20) vanishes and the coefficient of the quartic term is positive, yielding |A+| − 1 ∼ +(ωΔt)4 as ωΔt → 0, in agreement with the leading-order behavior shown in Fig. 5. Hence, the amplitude error per time step varies as (Δt)4 and the amplitude error per unit time varies as (Δt)3. Therefore, the numerical scheme with α = ½ is third-order accurate for amplitude errors and unconditionally unstable. In summary, the filter that conserves the three-time-level mean gives a numerical scheme that is two orders more accurate for amplitude errors than the standard Robert–Asselin filter. The increased accuracy may be exploited, by using it either to decrease the error for a given time step, or to increase the time step without increasing the error.
The unconditional instability of the case α = ½ may be avoided by instead choosing α ≳ ½, which almost conserves the three-time-level mean. The resulting filter is effectively a weighted blend of the third-order filter with α = ½ and the first-order filter with α = 1, the weighting of the latter being comparatively tiny. For this case, the coefficient of the quadratic term in (20) is negative, but very small, and the coefficient of the quartic term is positive. The negative quadratic term dominates for small ωΔt and the positive quartic term dominates for larger ωΔt, in agreement with the behavior shown in Fig. 5 for the case α = 0.53. Therefore, the numerical scheme with α ≳ ½ is conditionally stable. The finite stable range, for which |A+| ≤ 1, may be estimated from (20) by approximating the quartic term (but not the quadratic term) by its value when α = ½, to give
or, for ν ≪ 1,
The numerical scheme with α ≳ ½ is strictly only first-order accurate for amplitude errors. However, it is clear from Fig. 5 that the error for the case α = 0.53 is much smaller than the error of the first-order scheme with α = 1. Indeed, the error for the case α − ½ ≪ 1 is comparable in magnitude to the error of the third-order scheme with α = ½ across much of the finite stable range. Therefore, for practical purposes, the numerical scheme with α ≳ ½ is as good as third-order accurate for amplitude errors.
For all cases of α and ν in the practical range, the coefficient of the cubic term in (23) is positive, yielding arg(A+) − ωΔt ∼ +(ωΔt)3 as ωΔt → 0. Hence, the phase error per time step varies as (Δt)3 and the phase error per unit time varies as (Δt)2. Therefore, all numerical schemes in the modified family are second-order accurate for phase errors.
Table 1 summarizes the conservation, stability, and accuracy properties of the modified filter, when used in conjunction with the leapfrog scheme, for various values of α.
4. Summary and discussion
In the decades that have elapsed since the first general circulation models were developed, there have been major advances in the Earth observation systems from which initial conditions are derived, and in techniques for parameterizing unresolved physical processes. These advances have helped to improve the fidelity of weather and climate simulations. However, many general circulation models still use the same Robert–Asselin-filtered leapfrog time-stepping scheme as when they were first developed, despite evidence that time stepping appears to be an important contributor to model error.
This paper proposes a simple modification to the Robert–Asselin filter. The modified filter displaces the state at the future time slice as well as the current time slice. The modification yields a generalized family of filters, with a parameter controlling the relative sizes of the two displacements. The standard Robert–Asselin filter is a special case.
The behavior of the family of modified filters, when used in conjunction with the leapfrog scheme, is analyzed. The standard Robert–Asselin filter is shown to behave qualitatively differently from all other filters in the family. Each filter reduces the magnitude of the curvature at the three time slices operated upon, suppressing the computational mode. For the physical mode, each filter yields second-order accuracy for phase errors. But only the filter that conserves the three-time-level mean yields third-order accuracy for amplitude errors, whereas all other filters in the family (including the standard Robert–Asselin filter) yield only first-order accuracy. The filter that conserves the three-time-level mean yields an unconditionally unstable scheme, but conditional stability is recovered by adding a tiny amount of the standard Robert–Asselin filter, to yield a scheme that is as good as third-order accurate.
The modified filter proposed in this paper could improve weather and climate models. For example, it may permit the use of longer time steps with no loss of accuracy, reducing the computational expense of a given simulation. Alternatively, if the time step cannot be lengthened because it is constrained more strongly by other conditions (e.g., the CFL criterion) than by accuracy requirements, then the modified filter may permit an increase in accuracy at almost no additional computational expense. The modified filter would be extremely easy to implement in an existing computer model: the Robert–Asselin-filtered leapfrog routine could be upgraded by changing only a few lines of code. There may be a slight increase in the computational expense—the standard scheme is a two-stage method and the modified scheme is a three-stage method—but no extra function evaluations are required.
There are alternative methods for controlling the computational mode of the leapfrog scheme, which do not involve the application of a time filter. Kurihara (1965) proposed the leapfrog–trapezoidal method, which consists of obtaining a provisional value by applying a leapfrog predictor and then improving it by recursively applying a corrector. Being a predictor–corrector method, however, this scheme is iterative and potentially computationally demanding. Magazenkov (1980) proposed the alternate application, from one time step to the next, of a leapfrog step and a second-order Adams–Bashforth step. The need to execute a different algorithm at even- and odd-numbered time steps is cumbersome, however. In contrast, the modified filtered leapfrog scheme proposed in this paper is noniterative and nonalternating, yet still suppresses the computational mode and achieves third-order numerical accuracy for amplitude errors.
Déqué and Cariolle (1986) have shown that, despite the demonstrated ability of the standard Robert–Asselin filter to stabilize numerical solutions to atmospheric motion equations for certain combinations of temporal differencing and physical processes, in some other cases even a very weak filter may lead to an instability that can only be suppressed by a severe reduction of the time step. It remains to be seen whether the modified filter proposed in this paper also exhibits this unexpected behavior. Finally, it is possible that the modified filter could also improve the Robert–Asselin-filtered Adams–Bashforth schemes (e.g., Tandon 1987), but the exploration of this possibility is left for future work.
The author is funded through a Fellowship from the U.K. Natural Environment Research Council (NE/D009138/1). This study was supported in part by the National Science Foundation (PHY05-51164) and by a travel grant from the Royal Astronomical Society. The suggestions of two anonymous reviewers are gratefully acknowledged. The carbon footprint of this study—due to international air travel and office power consumption—is estimated by the author to be 3100 kg of CO2.
Corresponding author address: Paul D. Williams, Department of Meteorology, University of Reading, P.O. Box 243, Earley Gate, Reading RG6 6BB, United Kingdom. Email: email@example.com
Robert (1966) predates the standard citation databases.
* Kavli Institute for Theoretical Physics Report 08-95.