Abstract

An unmeteorological oscillation in the forecast fields produced by a two-time-level semi-Lagrangian and semi-implicit model has been traced to the method used for finding the departure point position. Various alternatives to the traditional method for determining it are examined and tested in a forecast model in order to find the one that removes the unmeteorological oscillation most satisfactorily.

1. Background

A problem with noise adjacent to the jet stream arose during testing of version 2.7 of the HIRLAM forecast model using an approximately 30 km × 30 km mesh length in the horizontal and 24 layers in the vertical. The time step was 10 min. [The two-time-level semi-Lagrangian and semi-implicit (2TLSLSI) model described in McDonald and Haugen (1993) was being used.] This seemed a reasonably conservative choice of time step since the operational model used a time step of 20 min for an approximately 55 km × 55 km mesh length in the horizontal and 16 layers in the vertical. That it was not can be seen from Fig. 1, which shows the vertical velocity at 300 hPa after a 24-h integration. The date of the starting analysis was 1200 UTC 24 February 1997. As can be seen from Fig. 2, which shows the 24-h forecast of horizontal velocity at 300 hPa valid at the same time, the noise is associated with a very strong jet.

Fig. 1.

The 24-h forecast of 300-hPa vertical velocity valid on 25 Feb 1997 at 1200 UTC.

Fig. 1.

The 24-h forecast of 300-hPa vertical velocity valid on 25 Feb 1997 at 1200 UTC.

Fig. 2.

The 24-h forecast of 300-hPa horizontal velocity valid on 25 Feb 1997 at 1200 UTC.

Fig. 2.

The 24-h forecast of 300-hPa horizontal velocity valid on 25 Feb 1997 at 1200 UTC.

This disconcerting phenomenon, which was not seen in the coarser mesh operational forecast, was immune to the standard “semi-Lagrangian noise abatement measures.” (a) Increasing the decentering coefficient from 0.1 to 0.4 (Rivest et al. 1994) had no impact on the noise. (b) Increasing the constant temperature associated with the semi-implicit scheme from 300 to 350 K (Simmons and Temperton 1997) had no impact, and increasing the associated constant pressure from 800 to 1013 hPa had only a small positive impact. (c) Integrating the nonlinear terms implicitly by means of an iterative scheme did not solve the problem, either.

Three different methods were found for controlling the noise; all were unsatisfactory for different reasons:(a) reducing the time step to 5 min (this doubles the computer time), (b) increasing the horizontal diffusion coefficient (but to a level where it was damping four to six grid waves significantly), (c) switching off the condensation scheme (no precipitation).

A closer examination of the noise showed it had an approximately 8Δt structure in time and an approximately four to five grid structure in space. The latter makes the effectiveness of the large horizontal diffusion coefficient understandable. The initial impulse was to“blame the physics,” since the problem went away when the physics, and in particular the condensation scheme, was switched off. However, closer examination lead to the conclusion that the dynamics was actually generating a small oscillation that the condensation scheme was amplifying significantly. This is illustrated in Fig. 3, which shows the evolution of the temperature at grid point (146, 42, 4) (approximately 46°N, 18°W, 116 hPa) during a forecast with a 10-min time step. The dots result from a forecast with no physics and no horizontal diffusion. Notice the small oscillation setting in after about 6 h. The diamonds result from a repeat of this forecast, but with the condensation scheme (but no other physics) switched on. Notice the much larger oscillations. The periods of both oscillations are similar enough to allow one to argue that the condensation scheme is amplifying the more modest oscillations produced by the dynamics.

Fig. 3.

Evolution of the forecast temperature at the grid point situated at approximately 46°N, 18°W, 116 hPa over the 24 h starting from 1200 UTC on 24 Feb 1997.

Fig. 3.

Evolution of the forecast temperature at the grid point situated at approximately 46°N, 18°W, 116 hPa over the 24 h starting from 1200 UTC on 24 Feb 1997.

It was stated in the paragraph before last that the dynamics was initially deemed innocent of blame for the noise. This is made understandable by the solid line in Fig. 3, which results from a “no physics” forecast run with the default horizontal diffusion (fourth-order implicit with K = 3 × 1013). There is almost no hint of an oscillation.

2. Estimating the departure point position

Hortal (1998) has indicated the possible source of this noise. He has demonstrated that the computation of the departure point position in the 2TSLSI scheme is “noise sensitive.” Consider the following one-dimensional advection equation,

 
formula

The 2TLSLSI discretization of this equation is given by

 
ψ{IΔx, (n + 1)Δt} − ψ{x∗, nΔt} = 0,
(2.2)

where x = IΔx, t = nΔt, and x∗ is obtained by solving the equation dx/dt = u(x, t) as

 
formula

Since u(x, t) is unknown and the integral is along the trajectory of the parcel, an approximate method must be devised for finding the departure point, x(t) = x∗. [The arrival point, x(t + Δt) = IΔx, is known]. A number of Ot2) schemes have been devised to do this.

One approach can be thought of as approximating the integral in Eq. (2.3) using the “midpoint rule,”

 
x(t) = x(t + Δt) − Δtu{x(t + Δt) + x(t)]/2, t + Δt/2},
(2.4)

and to evaluate x(t) iteratively; see McDonald (1991) for details. In order to make comparison with other methods transparent, the first two iterations are given explicitly in Eqs. (2.6) and (2.7). First, u is extrapolated to time level n + ½: u{x, (n + ½)Δt} = 1.5u{x, nΔt} − 0.5u{x, (n − 1)Δt}. Then

 
formula

Another approach can be thought of as approximating the integral in Eq. (2.3) using the “trapezoidal rule,”

 
formula

and then making iterative estimates of u{x(t), t}. Bartello and Thomas (1996) suggests the following: first, u is extrapolated to time level n + 1: un+1I = 2unIun−1I. Then

 
formula

Two other alternative integrations of Eq. (2.3) have been suggested. Both can be interpreted as using a variant of the midpoint rule,

 
x(t) = x(t + Δt) − Δtu{x(t + Δt/2), t + Δt/2},
(2.12)

and both use the same first estimate,

 
x(1) = IΔxunIΔt,
(2.13)

but can be interpreted as differing in their next estimate of u{x(t + Δt/2), t + Δt/2}. Temperton and Staniforth [1987, see their Eq. (42)] extrapolate along the trajectory, yielding the following estimate:

 
formula

Hortal (1998), on the other hand, uses the following estimate:

 
formula

which has the virtue over Eq. (2.14) of needing one less interpolation to find the departure point value.

There are other combinations of arrival point and departure point quantities that combine to yield x(2) to Ot2) accuracy. In fact, any combination of the form

 
x(2) = IΔx − Δ(1),
(2.16)

with

 
formula

with a = 1 + d + e + 2f, b = −½ − df, and c = ½ − d − 2e − 2f, will do so. To show this, rewrite Eq. (2.2) as

 
ψ{IΔx, (n + 1)Δt} − ψ{IΔx − Δtû, nΔt} = 0,
(2.18)

and expand ψ about (IΔx, nΔt), giving, for a quadratic or any higher-order interpolation scheme,

 
formula

If

 
formula

then the scheme becomes O2) accurate, as can be seen by substituting Eq. (2.20) into Eq. (2.19), giving

 
formula

See McDonald (1987) for details. Expanding all the terms on the right-hand side of Eq. (2.17) about (IΔx, nΔt) and matching them to Eq. (2.20) yields the conditions on a, b, and c above.

Further iterations can be performed; however, the scheme remains O2) accurate:

 
formula

There are many alternative ways of extrapolating. An intruiging one that still only uses two time levels is

 
û(k+1) = 1.5u{IΔx(k)Δt, nΔt} − 0.5u{IΔx(k)Δt, (n − 1)Δt},
(2.23)

where g = 3h − 1 guarantees Ot2) accuracy: h = ½ corresponds to Eq. (2.5) if unI is used in Eq. (2.6); h = 1 corresponds to “extrapolation along the trajectory.” One can also consider combining the approach of Eq. (2.22) and Eq. (2.23), and also of introducing additional time levels.

A good way to see where this may lead is to look at the equation dx/dt = u(x, t) from a “Runge–Kutta point of view” by writing the solution as (e.g., see chapter 5 of Lambert 1991)

 
formula

The accuracy of the scheme is dictated by the choice of s, aij, bi, ci. For example, looked at from this point of view, Eq. (2.11) has s = 2, k1 = u(xn+1, tn+1), and k2 = u(xn+1 − Δtk1, tn+1 − Δt); b1 = b2 = ½. For higher-order schemes (s > 2), the number of options is enormous. Alsos (1998) has demonstrated improved accuracy by applying some of the “classic” higher-order Runge–Kutta-style schemes to trajectory computations in semi-Lagrangian integrations.

3. Experimental results

In this section, when the equations in section 2 are referred to, it is their obvious three-dimensional generalization that is actually being used in the integrations described.

In order to study the noise produced by the dynamics in more detail the forecast described in section 1 was repeated with all of the physical parameterization schemes switched off. Also, in order to amplify the noise, the time step was increased to 15 min and the horizontal diffusion was switched off. Then using Eq. (2.7) to compute the departure point positions leads to a forecast with serious noise problems in the 24-h forecast of the vertical velocity at 300 hPa, along the jet and adjacent to it; see Fig. 4.

Fig. 4.

The 24-h forecast of 300-hPa vertical velocity valid on 25 Feb 1997 at 1200 UTC using Δt = 15 min, no physics, and no horizontal diffusion.

Fig. 4.

The 24-h forecast of 300-hPa vertical velocity valid on 25 Feb 1997 at 1200 UTC using Δt = 15 min, no physics, and no horizontal diffusion.

A convenient crude measure of this noise is the square root of the vertical velocity squared at 300 hPa summed in the horizontal over all of the grid points and divided by N, the total number of grid points: 〈ω〉 = {Σi,jω2i,j}1/2/N (〈ω〉 = 0.364 Pa s−1 for Fig. 4, and 〈ω〉 = 0.147 Pa s−1 for Fig. 1).

Equation (2.7) does not exactly describe the scheme used for the forecasts discussed in section 1. For those, a different first guess was used, which is equivalent to replacing u{IΔxun+1/2IΔt/2, (n + ½)Δt} with u{IΔxûΔt/2, (n + ½)Δt}, where û was retained from the previous time step, a procedure first suggested by Robert (1981). This makes matters worse: 〈ω〉 = 0.508 Pa s−1.

Using Eq. (2.11) produces a 24-h forecast with even more noise; 〈ω〉 = 0.739 Pa s−1. The noise can be reduced in this case by not extrapolating in time for the first guess. That is, using u{IΔxunIΔt, nΔt} instead of u{IΔxun+1IΔt, nΔt} in Eq. (2.11) reduces 〈ω〉 to 0.548 Pa s−1.

On the other hand, Temperton and Staniforth’s (1987) scheme, Eq. (2.14), produces almost noise-free forecasts, 〈ω〉 = 0.138 Pa s−1. The three-dimensional generalization of Eq. (2.15), Hortal’s (1998) scheme, also reduces the noise, but not as impressively as Temperton and Staniforth’s: 〈ω〉 = 0.162 Pa s−1.

Can we do any better, while maintaining Ot2) accuracy? To answer this the coefficients d, e, and f in Eq. (2.17) were varied systematically in units of 0.25 and the forecasts were repeated, with 〈ω〉 being recorded for each 24-h forecast. Two clear messages emerged. (a) Extrapolations that use u{IΔx, (n − 1)Δt} with anything other than a small positive coefficient produce a noisy forecast. (b) Extrapolations that use u{IΔx, nΔt} with anything other than a small coefficient produce a noisy forecast. Using the somewhat arbitrary definition of 〈ω〉 ≥ 0.160 Pa s−1 as “too noisy,” then all the forecasts with 〈ω〉 less than 0.160 had 0.0 ⩽ b ⩽ 0.25, −0.5 ⩽ a ⩽ 0.5, and 0.0 ⩽ e ⩽ 1.0.

The forecasts with the 10 smallest values of 〈ω〉 are listed in Table 1 (see column 8). It must be emphasized that the value of 〈ω〉 is only a crude guide to the amount of noise on the chart. Nevertheless, examining the charts in detail confirms number 1 in Table 1 as the “winner,” from a subjective viewpoint also.

Table 1.

The value of 〈ω〉 for the 10 least noisy forecasts along with the corresponding values of a, b, c, d, e, f in Eq. (2.17); 〈ω15 refers to forecasts using a 15-min time step, 〈ω20 refers to forecasts using a 20-min time step.

The value of 〈ω〉 for the 10 least noisy forecasts along with the corresponding values of a, b, c, d, e, f in Eq. (2.17); 〈ω〉15 refers to forecasts using a 15-min time step, 〈ω〉20 refers to forecasts using a 20-min time step.
The value of 〈ω〉 for the 10 least noisy forecasts along with the corresponding values of a, b, c, d, e, f in Eq. (2.17); 〈ω〉15 refers to forecasts using a 15-min time step, 〈ω〉20 refers to forecasts using a 20-min time step.

To gain some insight into what makes a forecast less noisy, the different extrapolations were plotted on a“space–time diagram,” see Fig. 5, where time is plotted on the y axis and departure point distance on the x axis. In order to estimate û to Ot2) accuracy we must extrapolate to the point “O,” located at {IΔxûΔt/2, (n + ½)Δt}. To do this we are using combinations of the values of u located at the points A, B, C, D, E, F. Then, for the scheme of Eq. (2.14), which has û(1) = (3unun−1)/2, extrapolation along the trajectory is represented by the line FCO in Fig. 5. (In this section “*” is shorthand notation for {IΔxunIΔt} and “**” for {IΔx2unIΔt}.)

Fig. 5.

Space–time representation of the different extrapolations used to find the departure point position.

Fig. 5.

Space–time representation of the different extrapolations used to find the departure point position.

What is interesting is that all the combinations of coefficients that have smaller 〈ω〉 than this scheme have smaller slopes than it in the space–time diagram. Consider number 1 in Table 1 as an example:

 
formula

The three terms within the first square brackets describe an “interpolation” to the point T while the two in the second square brackets describe an extrapolation to Q in Fig. 5. Thus the extrapolation to the departure point is represented by QTO. Looked at from this point of view all the other eight extrapolations listed in Table 1 can be represented by PSO, QTO, or RUO; see column 10 of Table 1.

By the same token we can rewrite û in Eq. (2.11) as

 
formula

which can be viewed as the extrapolation BJO in Fig. 5. Recall that this extrapolation produced quite a noisy forecast (〈ω〉 = 0.548 Pa s−1). Also, the extrapolation described by Eq. (2.7), which was less noisy (〈ω〉 = 0.364 Pa s−1), can be viewed as LMO. The extrapolation described by Eq. (2.15), which is even less noisy (〈ω〉 = 0.162 Pa s−1), can be viewed as DKO.

It is of interest to ask whether these trends were also shown by the extrapolation of Eq. (2.23). Does h = 4/3 (QTO in Fig. 5) give the least noisy forecast? The answer is no (〈ω〉 = 0.223 Pa s−1). The extrapolation along the trajectory (h = 1) is the least noisy.

In order to see whether these results were dependent on the time step the forecasts were rerun for the 10 best cases with a smaller (Δt = 10 min) and larger time step (Δt = 20 min). For Δt = 10 min there was very little to choose from between the different schemes. All had approximately the same value of 〈ω〉. For Δt = 20 min, however, there was one striking result. The value of 〈ω〉 increased significantly for number 5 in Table 1 (see column 9). Looking at the 24-h forecast of 300-hPa vertical velocity confirms an outbreak of noise west of Spain. Aside from this exception the order of merit, in general, remains the same (no. 2 now the lowest value of 〈ω〉, and no. 1 the second lowest, but the difference is marginal).

Is the reduction of noise maintained when the physics is switched on? To answer this, the forecast described in section 1 was repeated with the departure point computation being replaced by Eq. (2.17) with the coefficients taken from row number 1 in Table 1. This resulted in the 24-h forecast of 300-hPa vertical velocity shown in Fig. 6, which has 〈ω〉 = 0.110 Pa s−1, contrasting with 〈ω〉 = 0.147 Pa s−1 for Fig. 1.

Fig. 6.

The 24-h forecast of 300-hPa vertical velocity valid on 25 Feb 1997 at 1200 UTC, using the new method for computing the departure point.

Fig. 6.

The 24-h forecast of 300-hPa vertical velocity valid on 25 Feb 1997 at 1200 UTC, using the new method for computing the departure point.

4. Stability analysis

The results described in sections 1 and 3 indicate the presence of an unstable, or at the very least an oscillating, mode generated by the departure point computation. In this section an attempt is made to gain some insight into its origin.

There is no obvious way to perform a stability analysis of the highly nonlinear equation dx/dt = u(x, t) combined with Eq. (2.1). Therefore, let us resort to the standard procedure of assuming that we can gain some useful insight by investigating the oscillation equation,

 
formula

where u0 is a constant advecting velocity, as a first-order linear representation of dx/dt = u(x, t).

Let us assume a simple discretization, which although only Ot) accurate, illustrates two interesting phenomena:

 
χ{IΔx, (n + 1)Δt} − χ{IΔxu0Δt, nΔt} = iνΔtχ{IΔxu0Δt, nΔt}.
(4.2)

For simplicity, to find the departure point values of χ, let us assume a linear interpolation using the surrounding grid points, (Ip − 1) and (Ip):

 
χn+1I = (1 + Δt)[(1 − α̂)χnIp + α̂χnIp−1]
(4.3)

(α̂ = u0Δtxp). Substituting the solution χnI = χ0λn exp(ikIΔx) gives

 
formula

The details are the same as in Bates and McDonald (1982).

The first factor in Eq. (4.6) is always greater than 1, even for zero advecting velocity (α̂ = p = 0), for which case the scheme reduces to an Euler forward scheme. Mesinger and Arakawa (1976) point out that, “although this satisfies the Von Neumann necessary condition for stability, experience shows that indiscriminate use of this scheme leads to amplification at an unacceptable rate.” This leads to conjecture number 1.

1) The underlying “forward in time” nature of the schemes used when estimating the departure point position, imposed by the 2TLSLSI approach, means that there will be a tendency toward instability in integrations that use it.

The second factor in Eq. (4.6) is always less than or equal to 1 and will thus tend to stabilize the scheme no matter how large the advecting velocity (|λ|2 is independent of p). As has been pointed out by D. Durran in a slightly different context (1998, personal communication, but see his forthcoming book) this is not a general result, but depends on the fact that the departure point in the term multiplied by ν in Eq. (4.2) is on the trajectory. To see this at its simplest, substitute iνχ{IΔx, nΔt} for iνχ{IΔxu0Δt, nΔt} in Eq. (4.2). Then the factoring of exp(−ikpΔx) that occurred in Eq. (4.5) is absent and |λ|2 depends on p. This leads to conjecture number 2.

2) In estimating the departure point position in the 2TLSLSI scheme only extrapolations along the trajectory will have stability criteria independent of the number of grid points between the arrival and departure points.

In section 2, therefore, only the scheme described by Eq. (2.14) fulfils this criterion. Looking at row 10 of Table 1 and comparing column 8 with column 9 in Table 1, we see the “measure of noise” is unchanged, possibly a confirmation of conjecture (2).

Because |λ| depends on p does not necessarily mean that the scheme becomes more amplifying as p increases, although we saw a possible example of this in section 3. See row 5 of Table 1. It could become more damping as p increases, and, in fact, with the exception of row 5, the measure of noise decreased with increasing time step for all the other cases in Table 1. This might be an advantage since increasing Δt will make the first factor in Eq. (4.6) more amplifying.

Assuming the most general form, Eq. (2.17), for the right-hand side of Eq. (4.1) leads to

 
λ2 − λ[ρ + iνΔt(a + cρ + eσ)] − iνΔt(b + dρ + fσ) = 0
(4.7)

with

 
σ = [(1 − â) + â exp(ikΔx)] exp(−iksΔx),
(4.8)

where (â = 2u0Δtxs). Solving for |λ| numerically it is possible to confirm conjecture 2 for a = b = d = e = 0.

5. Discussion

The most popular method used to find the departure point position has been shown to be a source of noise in a reasonably fine-mesh (30 km) integration of the HIRLAM forecast model using a dataset that contains a strong upper-air jet. An assortment of alternative methods for computing the departure point position have been tested in the forecast model and a number of candidates have been found that give noise-free (or at very least “very low noise”) forecasts. One of these, extrapolation along the trajectory, has the attractive feature that its stability is independent of the number of grid points between the departure and arrival points—albeit in an analysis that bears a somewhat tenuous connection to the highly nonlinear problem being examined. It is with dismay that one notes that Temperton and Staniforth (1987), in shallow water experiments, found that with this scheme, “results were consistently worse than those for the simpler method 2” [their method 2 corresponds to our Eq. (2.5)]. This is difficult to understand since both schemes are Ot2) accurate but it emphasizes the point that the experiments described in this paper concentrated on noise abatement, and not on forecast accuracy. (Comparing the full physics forecasts using Δt = 10 min with the equivalent Eulerian forecast with Δt = 2 min, the only accuracy test I have performed, showed no hint of the discrepency described by them.)

There is an additional expense involved in using a scheme that requires evaluation of the velocity fields at both {IΔxunIΔt} and {IΔx2unIΔt}. For the experiments described in section 3 one iteration was used to find the departure point; that is, the three-dimensional equivalent of Eq. (2.17) was used. This involved three additional trilinear interpolations. Thus, the increase in computational cost was minimal.

Extrapolating along the trajectory to find the departure point may well be dangerous when that trajectory passes over rapidly varying orography. Although no evidence of bad behavior over steep orography was seen in the experiments described in this paper, this might manifest itself as a problem for finer grids with many vertical layers.

The noise in the forecast is not a “classic instability,” which blows up after a few time steps. As seen in this dataset it might be more accurately described as an“unmeteorological oscillation.” There is a pressing need for a better attempt at analyzing this effect than the preliminary investigation presented in section 4.

Acknowledgments

Thanks to Ivar Lie and two anonymous reviewers whose constructive comments strengthened the original manuscript.

REFERENCES

REFERENCES
Alsos, T., 1998: Effective ODE-solvers for trajectory calculations in semi-Lagrangian methods used in weather forecasting. Norwegian University of Science and Technology, 53 pp. [Available from Department of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim, Norway.]
.
Bartello, P., and S. J. Thomas, 1996: The cost-effectiveness of semi-Lagrangian advection. Mon. Wea. Rev., 124, 2883–2897.
Bates, J. R., and A. McDonald, 1982: Multiply-upstream, semi-Lagrangian advective schemes: Analysis and application to a multi-level primitive equation model. Mon. Wea. Rev., 110, 1831–1842.
Durran, D. R., 1998: Numerical Methods for Wave Equations in Geophysical Fluid Dynamics. Springer-Verlag, 465 pp.
Hortal, M., 1998: Some recent advances at ECMWF. LAM News., 27, 32–36.
Lambert, J. D., 1991: Numerical Methods for Ordinary Differential Systems. John Wiley and Sons, 293 pp.
McDonald, A., 1987: Accuracy of multiply-upstream semi-Lagrangian advective schemes II. Mon. Wea. Rev., 115, 1446–1450.
——, 1991: Semi-Lagrangian methods. Proc. ECMWF Seminar on Numerical Methods in Atmospheric Models, Vol. I, Reading, United Kingdom, ECMWF, 257–271. [Available from A. McDonald, Irish Meteorological Service, Glasnevin Hill, Dublin 9, Ireland.]
.
——, and J. E. Haugen, 1993: A two-time level, three-dimensional semi-Lagrangian and semi-implicit grid point model of the primitive equations. Part II: Extension to hybrid coordinates. Mon. Wea. Rev., 121, 2077–2087.
Mesinger, F., and A. Arakawa, 1976: Numerical Methods Used in Atmospheric Models. GARP Publ. Series, Vol. 17, WMO, 64 pp.
Rivest, C., A. Staniforth, and A. Robert, 1994: Spurious resonant response of semi-Lagrangian discretizations to orographic forcing: Diagnosis and solution. Mon. Wea. Rev., 122, 366–376.
Robert, A., 1981: A stable numerical integration scheme for the primitive meteorological equations. Atmos.–Ocean, 19, 35–46.
Simmons, A. J., and C. Temperton, 1997: Stability of a two-time level semi-implicit integration scheme for gravity wave motions. Mon. Wea. Rev., 125, 600–615.
Temperton, C., and A. Staniforth, 1987: An efficient two-time-level semi-Lagrangian semi-implicit integration scheme. Quart. J. Roy. Meteor. Soc., 113, 1025–1039.

Footnotes

Corresponding author address: Dr. Aidan McDonald, Irish Meteorological Service, Glasnevin Hill, Dublin 9, Ireland.