## 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.

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.

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 × 10^{13}). 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,

The 2TLSLSI discretization of this equation is given by

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

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 *O*(Δ*t*^{2}) 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,”

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.5*u*{*x, n*Δ*t*} − 0.5*u*{*x,* (*n* − 1)Δ*t*}. Then

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

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: *u*^{n+1}_{I} = 2*u*^{n}_{I} − *u*^{n−1}_{I}. Then

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

and both use the same first estimate,

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:

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

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 *O*(Δ*t*^{2}) accuracy. In fact, any combination of the form

with

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

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

If

then the scheme becomes *O*(Δ^{2}) accurate, as can be seen by substituting Eq. (2.20) into Eq. (2.19), giving

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 *O*(Δ^{2}) accurate:

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

where *g* = 3*h* − 1 guarantees *O*(Δ*t*^{2}) accuracy: *h* = ½ corresponds to Eq. (2.5) if *u*^{n}_{I} 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)

The accuracy of the scheme is dictated by the choice of *s, a*_{ij}, *b*_{i}, *c*_{i}. For example, looked at from this point of view, Eq. (2.11) has *s* = 2, *k*_{1} = *u*(*x*_{n+1}, *t*_{n+1}), and *k*_{2} = *u*(*x*_{n+1} − Δ*tk*_{1}, *t*_{n+1} − Δ*t*); *b*_{1} = *b*_{2} = ½. 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.

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} *ω*^{2}_{i,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*Δ*x* − *u*^{n+1/2}_{I}Δ*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*Δ*x* − *u*^{n}_{I}Δ*t, n*Δ*t*} instead of *u*{*I*Δ*x* − *u*^{n+1}_{I}Δ*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 *O*(Δ*t*^{2}) 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.

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 *O*(Δ*t*^{2}) 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)} = (3*u*^{n}_{∗} − *u*^{n−1}_{∗})/2, extrapolation along the trajectory is represented by the line FCO in Fig. 5. (In this section “*” is shorthand notation for {*I*Δ*x* − *u*^{n}_{I}Δ*t*} and “**” for {*I*Δ*x* − 2*u*^{n}_{I}Δ*t*}.)

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:

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

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.

## 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,

where *u*_{0} 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 *O*(Δ*t*) accurate, illustrates two interesting phenomena:

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

(*α̂* = *u*_{0}Δ*t*/Δ*x* − *p*). Substituting the solution *χ*^{n}_{I} = *χ*_{0}*λ*^{n} exp(*ikI*Δ*x*) gives

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*Δ*x* − *u*_{0}Δ*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

with

where (*â* = 2*u*_{0}Δ*t*/Δ*x* − *s*). 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 *O*(Δ*t*^{2}) 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*Δ*x* − *u*^{n}_{I}Δ*t*} and {*I*Δ*x* − 2*u*^{n}_{I}Δ*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

## Footnotes

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

Email: aidan.mcdonald@met.ie