The one-dimensional advection equation and the one-dimensional advection adjustment equation with rotation are used as test beds for the design of well-posed boundary conditions for the initial-boundary-value problem using semi-Lagrangian discretization. Three options are found to be stable in experimental tests: trajectory truncation, time interpolation, and a well-posed buffer zone. Stability is proved for all three for the one-dimensional advection equation when linear interpolation is used for the interpolation associated with the semi-Lagrangian discretization.
In a limited area model the lateral borders are not physical boundaries to the flow but, instead, artificial constructs imposed by our desire to model what is happening in a subvolume of the global atmosphere. In order to solve the equations describing the evolution of the atmosphere in that subvolume we must supply values for the fields at these artificial boundaries. Doing this correctly is a rather subtle matter. As Oliger and Sundström (1978) have shown, for each inward pointing characteristic a field must be supplied at the lateral boundary. This implies that only a certain subset of all the variables should be imposed. If this subset has been correctly chosen the initial-boundary-value problem is “well posed.”
Constructing a well-posed system for the primitive meteorological equations is a daunting task. For that reason, and because an easy-to-implement alternative solution to the boundary problem exists, numerical weather prediction models with well-posed lateral boundaries have never gotten farther than the “drawing board” stage of existence; see, for example, Elvius and Sundström (1973). The above-mentioned alternative is to overspecify the boundaries and damp the resulting noise with a relaxation scheme; see Davies (1976). This gives stable forecasts and is almost trivial to implement. Evidence of its popularity is its almost universal use in operational limited area numerical weather prediction models. It is also used in many research models.
Despite this there are a number of reasons for taking a fresh look at well-posed boundaries. (i) There is evidence that the flow relaxation scheme can cause mass loss or gain; see McDonald (1998). (ii) At points at which the characteristic velocity is pointing out of the domain of integration we are imposing the external model solution when we do not need to if we use the flow relaxation scheme. Thus we may introduce unnecessary errors when these fields are incorrect (as they will be in an operational environment). (iii) Consider the situation where both the “host model” (the model supplying the boundary fields for the limited area “guest model”) and guest model fields are each in geostrophic balance. Then, as Cats and Åkesson (1983) demonstrate in their appendix the relaxation scheme can destroy this balance throughout the boundary relaxation zone. (iv) In the future 4DVAR (Le Dimet and Talagrand 1986) will be used in the forecast–analysis cycle. This has implications for our boundary treatment. Because of the iterative nature of 4DVAR the influence of the boundary conditions spreads farther into the limited area domain with each iteration. When the adjoint model is integrated backward in time, what were formerly boundary points with incoming characteristic velocities may become points describing error gradient fields with outgoing characteristic velocities. These must now pass out of the limited area without reflection. If our boundary treatment is only partially successful in doing this then one may speculate that this error could build on itself as we continue the forecast–analysis cycle. See Gustafsson et al. (1998). The Davies relaxation scheme may not be sufficiently subtle to accommodate these requirements. (v) In the past the operational limited area philosophy could reasonably be described as follows. The boundaries were placed far enough away from the region of interest that the boundary errors would not have time to corrupt the forecast. The source of these boundary errors was twofold: (a) the boundary formulation, and (b) the inaccuracy of the host model fields (the discretization mesh was coarse, the fields were “old,” they were available only every 6 h, possibly). For the future, with increasing computing power, a new scenario can be envisaged. The host model fields will be sufficiently accurate to describe the broad-scale features of the atmosphere accurately (highs, lows, fronts). The purpose of the operational limited area model will be to forecast the fine details of the meteorological phenomena associated with these broad-scale features (rainfall, surface temperatures and winds, fog, cloud information, pollution dispersion, etc.). In a nesting situation, for example, the finest mesh area could be so small that a front can cross the area during the forecast period. Thus in the later stages the forecast will be dominated by the boundary conditions. Now, if there are errors in our boundary formulation they will corrupt the forecast.
For these reasons the High Resolution Limited Area Model Group has set itself the goal of building a numerical weather prediction model that has well-posed boundaries. This paper takes a first step in that direction by looking at the problems associated with the “multiply upstream well-posedness” problem. That is, how do we invent well-posed boundary conditions when using semi-Lagrangian semi-implicit discretizations? This question is addressed for two one-dimensional systems:the advection equation, and the advection–adjustment equation. Three possible answers are examined: trajectory truncation, time interpolation, and a well-posed buffer zone. Examining the one-dimensional advection equation (see section 2) allows us to ask the following question: how do we maintain stability when the departure point is outside the integration area, in isolation from the complications introduced by the C grid and the semi-implicit scheme? These latter complications are considered in section 3. In sections 2 and 3 we demonstrate stability experimentally. To complement this, in section 4 we make a first attempt at establishing stability theoretically. We do so by applying the energy method to the one-dimensional advection equation for our three postulated solutions for the multiply upstream well-posedness problem. Last, section 5 contains a discussion of the strengths and weaknesses of these different approaches to the boundary value problem.
Naïve extrapolations from the interior to the boundary are risky; see section 6 of Gustafsson et al. (1972). Therefore, as a design strategy, they are introduced only when unavoidable in what follows.
2. The one-dimensional advection equation
The one-dimensional advection equation is ideal for looking at the semi-Lagrangian discretization in isolation. Well-posedness is transparent for the analytical system, and complications introduced by the C grid, the semi-implicit scheme, and higher spatial dimensions are absent. Therefore let us start with it:
where 0 ⩽ x ⩽ L, t ≥ 0, and u is a constant velocity. An initial condition is needed: ϕ(x, 0) = f(x). Also, a boundary condition is needed: ϕ(0, t) = g(t), if u > 0. If u < 0, the boundary condition must be supplied at L.
a. Discretization of the equations
If the solution enters through the boundary at x = 0 and exits through the boundary at x = L without distortion or reflection, the discretization is well posed. For that discretization let us use the semi-Lagrangian scheme. The solution to Eq. (2.1) is
where . The arrival point is at iΔx and the departure point, x∗, is at iΔx − uΔt. With a Lagrange cubic interpolation,
where α = uΔt/Δx; p = [α] if α ≥ 0 and p = [α − 1] if α < 0; and [α] is the integer part of α; p + α̂ = α. To maintain stability p must be chosen such that (i − p − 1)Δx ⩽ x∗ ⩽ (i − p)Δx or equivalently 0 ⩽ α̂ ⩽ 1. See Bates and McDonald (1982) for more details.
If the departure point is outside the area we cannot use Eq. (2.3). We examine various methods for addressing this problem below. Nor can we use Eq. (2.3) when the departure point lies in the grid box immediately next to the boundary. In that case let us use a quadratic Lagrange interpolation:
when 0 ⩽ x∗ < Δx and
when L − Δx ⩽ x∗ < L.
Let us consider the outflow boundary first. Recall that for well posedness we must not impose a value of ϕ at this boundary; we must extrapolate from the interior. The semi-Lagrangian scheme provides us with an ideal way of doing this; see Robert and Yakimiw (1986). When |α| < 1, then we use Eq. (2.4). It yields an O(Δx3/Δt) accurate solution. If |α| ≥ 1, then we use Eq. (2.3), which is O(Δx4/Δt) accurate.
It is interesting to note that if, for example, we were using the leapfrog Eulerian scheme to integrate Eq. (2.1) we would have to switch to either Eq. (2.4) or some other upwind type of discretization in order to compute ϕ at the outflow boundary point. See McDonald (1999a), for example.
Next consider the inflow boundary and points in its immediate vicinity. On the boundary line itself we must impose a value on the field. If |α| > 1, then updating ϕ at the grid points adjacent to the inflow boundary causes difficulties, which we examine in sections 2d–2f. If |α| ⩽ 1, then p = 0, and we proceed as in section 2c.
b. Setup for the numerical tests
To test this let us model a bell curve of “width” L/10 being advected in the positive direction along the x axis with a velocity u, starting with the maximum of the bell curve at x = xs. At any subsequent time t the position of that maximum will be at x(t) = xs + ut. The initial condition is
Consider u > 0 without loss of generality. Then we must supply a boundary condition at x = 0:
The following parameters will be used in the tests of our discretization: Δx = 10 km, u = 20 m s−1, L = 1000 km, and ϕ̂ = 10. We will perform the following two experiments.
Test 1: start with the bell curve centered at x = L/2 and choose the time such that at the end of the integration it will be centered at 3L/2 in order to see whether it passes through the boundary at x = L without false reflections.
Test 2: start with the bell curve centered at x = −L/2 and choose the time such that at the end of the integration it will be centered at L/2 in order to see whether it passes through the boundary at x = 0 without distortion.
c. Testing with |α| < 1
With α = 0.5 it passes both tests. Figure 2 shows the results of test 1. The bell curve has exited through the boundary accurately. There are no false reflections (the dots coincide exactly with the diamonds). Figure 3 shows the results of test 2. The bell curve has entered the area accurately (the dots coincide almost exactly with the diamonds).
Let us return to the problem near the inflow boundary when |α| > 1. We now have |p| > 0, and when, for example, u > 0, the departure point associated with ϕ1 is outside the area. What are the options?
d. Option 1, when |α| > 1: Trajectory truncation
If the departure point is beyond the boundary, then truncate the trajectory at the boundary. Thus, for example, if the departure point associated with ϕ2 were at x = −Δx/2 we would use for Eq. (2.2)
This method is stable but slightly inaccurate, as can be seen from Fig. 4, which shows the result of repeating test 2 with α = 2.5
e. Option 2, when |α| > 1: Time interpolation
Let us try to improve the accuracy of the trajectory truncation method. Recall that at inflow we know the field value at every time step on the boundary line, in principle, anyway. Let us use this information by invoking an idea of Skålin and Lie (1995). See also Falcone and Ferretti (1998). Consider the following argument. The solution to Eq. (2.1) can be written as
Choose uδt = μΔx, where μ is a nonzero integer. Also choose t + δt = (n + 1)Δt. Substituting in Eq. (2.8) we get
For the following argument consider u > 0. There is no loss of generality. Let i = μ; then this equation expresses the value of the field at time (n + 1)Δt at the interior grid point, μ, in terms of the field on the boundary line, i = 0, at time (n + 1 − μ/α)Δt. The field is not known at this time. However, we can “interpolate in time” because, in principle, on the boundary line we know ϕ(x0, mΔt) for all integer m > 0. For example, using a linear Lagrange interpolation for the right-hand side of Eq. (2.9) yields
This is stable from a von Neumann point of view, anyway, because 0 < μ/α < 1. When |α| < μ, we do not use this approach to finding the departure point value of the field; we use Eq. (2.3) or (2.4). Equation (2.9) is only used when i < [α], where [α] is the integer part of α. Thus, in Eq. (2.10), where we have chosen i = μ, μ/α ⩽ 1.
If, instead, we have inflow at the x = L boundary we choose i = μ + iL. Both μ and α are now negative but the arguments are the same and instead of Eq. (2.10) we use
Obviously, O(Δt2) accuracy can be achieved by using a Lagrange quadratic interpolation,
Repeating test 2 with α = 2.5, using quadratic interpolation in time at the inflow boundary, we see from Fig. 5 that full accuracy is restored.
f. Option 3, when |α| > 1: Well-posed buffer zone
From a semi-Lagrangian point of view it would be ideal if we had a buffer zone external to the boundary. Then if the departure point were located there, we could interpolate in the same way as if it were in the interior. In this section we generate such a buffer zone, which gives stable integrations by using only information from the inflow boundary line. Again, for the purpose of demonstration consider the situation when u > 0. Then we wish to compute ϕ−1, ϕ−2, etc., in Fig. 1.
Naïve extrapolation from the interior to the inflow buffer zone will cause instability. Thus we cannot, for example, use Eq. (2.3) to extrapolate, abandoning the restriction 0 < α̂ ⩽ 1. [A von Neumann stability analysis shows that for α̂ = −2 the amplification factor is 15 for the two-grid wave; see Eq. (23) of McDonald (1984).]
Since for well-posedness we must supply ϕ(0, t) for all t we know ∂ϕ(0, t)/∂t, ∂2ϕ(0, t)/∂t2, and so on. Thus, from Eq. (2.1) we know
and so on, enabling us to estimate ϕ−1, ϕ−2, etc., using a Taylor expansion and substituting from Eq. (2.14):
When u < 0, the same procedure can be used to compute ϕ(L + iΔx).
Repeating test 2 with α = 2.5, using this procedure to compute the necessary fields in the boundary zone restores full accuracy. The figure (not shown) looks identical to Fig. 5.
3. The one-dimensional advection-adjustment equation
We saw in section 2 that we could integrate the one-dimensional advection equation stably for large |α| by using time interpolation, a well-posed buffer zone, or, less accurately, trajectory truncation. In the real atmosphere adjustment and rotation also play significant roles. Thus in this section we add adjustment and rotation terms to the advection equation in order to see whether these schemes can still be implemented in this more realistic environment.
a. Discretization of the equations
The equations we wish to discretize are
where the advection velocity u, the Coriolis parameter f, and the mean geopotential height ϕ are constants.
The characteristic variables of these equations (and their associated characteristics) are
For each inward pointing characteristic a field must be imposed on the boundary; see, for example, Sundström and Elvius (1979).
Substituting ϕ = ϕ̂ exp[i(kx − ωt)] into Eqs. (3.1)–(3.3) yields three solutions; a “slow” solution with ω = ku and two “fast” solutions with ω = k(u ± ck), where
where Xni,* means interpolate X to the departure point associated with the grid point i. Having computed Rϕ, Ru, and Rυ, one approach to solving this system of equations is as follows. Equations (3.8)–(3.10) can be rearranged to yield a tridiagonal system for un+1. Having solved this, ϕn+1 and υn+1 emerge from Eqs. (3.8) and (3.10), respectively.
If we were to proceed in this way for the two-dimensional shallow water equations we would generate the doubly averaged terms and , each of which has a nine-point template. Near the boundary these terms would further complicate the already difficult problem of inventing a stable discrete system of equations. With this in mind let us use the following iterative approach to solving these equations; its two-dimensional generalization will cause minimal complications on the boundary:
Having solved this equation for ϕ with appropriate boundary conditions, u(k+1)i+1/2 can be computed from Eq. (3.15), and υ(k+1)i from Eq. (3.16). The procedure can be iterated until sufficient accuracy is attained.
The indices for which Eqs. (3.14)–(3.17) are valid depend on the grid points at which the boundaries are imposed. Examples are given in Table 1 for four different boundary choices: u1/2, uiL+1/2 imposed; ϕ0, ϕiL imposed; p0, qiL imposed; and p1/2, qiL+1/2 imposed. As we shall see in section 3b, for example, it will be necessary to approximate Eq. (3.3) on the boundary line. Because Eq. (3.10) cannot be used we must replace it with an equation that uses a stable and accurate extrapolation from the interior of the required unknown fields.
b. Trajectory truncation
To clarify the discussion let us impose our boundary values on ϕ0 and ϕiL and υ at inflow. They fit naturally into Eq. (3.17). There is a complication: we need to know υ on the outflow boundary line in order to solve Eq. (3.17). Also, as we discuss below, υ at inflow is problematical. Thus, we are forced into making extrapolations from the interior.
Consider the problem of computing υ at outflow. We can guess that “semi-Lagrangian extrapolation” will be stable. That is, we assume an approximation to Eq. (3.3) is obeyed at the outflow boundary point. Notice we must extrapolate unknown fields from the interior; let us use the simplest possible approximation:
When the departure point is within two grid points of the boundary, we will need (Xυ)n0 or (Xυ)niL for the computation of Rυ. For consistency we define them as
In fact, as is shown in McDonald (2000) by experiment, the integration is inaccurate unless we make sure to use Xυ0 and XυiL at inflow as well. Thus trajectory truncation for these boundary conditions will mean using Ru = (Xu)n1/2, when x∗ < Δx/2; Ru = (Xu)niL−1/2, when x∗ > L − Δx/2; Rϕ = (Xϕ)n1, when x∗ < Δx; Rϕ = (Xϕ)niL−1, when x∗ > L − Δx; Rυ = (Xυ)n0, when x∗ < 0; and Rυ = (Xυ)niL, when x∗ > L.
Extrapolations from the interior are also necessary if we impose boundary conditions on u, p, or q. See McDonald (2000) for details.
c. Time interpolation
Here b represents a boundary line, and μ is the integer number of points the ϕ point of interest is from the boundary. Also, ν is a half odd-integer, which measures the distance of the u point of interest from the boundary.
Assume c > u. Then well-posedness implies that we impose two fields on the inflow boundary and one on the outflow boundary. To keep the arguments clear let us assume we are imposing ϕ and υ at inflow and ϕ at outflow. Then as can be seen, most of the fields required to compute X̂ϕ, X̂u, and X̂υ are not known on the boundary, even at time level n. How can we proceed?
The terms u, ∂u/∂x, and ∂ϕ/∂x are not known on the boundary and, thus, must be extrapolated. The simplest possible extrapolations give, at x = 0,
and at x = L,
How do we solve Eqs. (3.23)–(3.25)? We have a number of choices. (i) Use X̂n and X̂n−1 to extrapolate to time X̂n+1−μ/α (X̂n+1−ν/α). I have not been able to invent a stable method for doing this. (ii) Incorporate the time interpolation into the semi-implicit scheme. That is, use Eq. (2.12) to rewrite Eqs. (3.23)–(3.25), resulting in the following:
Even using the simplest of extrapolations given in Eqs. (3.26) and (3.27) the discretization of Eqs. (3.28)–(3.30) leads to a set of equations that are unattractive from the point of view of programming efficiency. The complications are as follows. (a) Both the left- and right-hand sides of the discretized versions Eqs. (3.28)–(3.30) can formally differ for each half-integer increment of α. As is shown in the appendix of McDonald (1999b) there are three alternative equations for u1/2, depending on the size of α. If we use the next highest order of approximation for u0, there will be five. (The programming complication being that we need “if” statements or their equivalents.) (b) The number of floating point operations increases significantly. (c) Because the semi-Lagrangian and semi-implicit systems have become intertwined we no longer have a tridiagonal system of equations to solve.
From a computing point of view these drawbacks may not be too serious in a single-processor scalar environment. We would simply pay a possibly moderate computational price for getting the boundaries right. However, with vector processors in a multiprocessor environment these complications seem very serious. For these reasons it is interesting to consider an alternative that avoids these pitfalls, that of iterating.
(iii) Move the troublesome terms in Eqs. (3.28)–(3.30) to the right-hand side and compute them iteratively:
When discretized, these equations give rise to a set that is formally the same as Eqs. (3.14)–(3.16) and is thus very easy to program. (The tridiagonal system now has variable coefficients and there are a few extra terms on the right side to be evaluated during each iteration.)
d. Well-posed buffer zone
How can we fill a well-posed buffer zone when some of the fields needed for the extrapolation cannot be imposed externally? For discussion purposes assume c > u. Then only two fields at inflow will be imposed externally. In order to compute fields for a buffer zone at inflow we need to know the third field, which to maintain well-posedness we must not impose externally. How can we estimate it stably using only fields from the interior?
Gustafsson et al. (1972) showed that, for the Lax–Wendroff scheme, many extrapolations cause instability on 0 ⩽ x ⩽ L; t > 0. We can infer that such extrapolations will also be unstable for the semi-Lagrangian schemes since the Lax–Wendroff scheme corresponds to the latter with quadratic interpolation when |α| < 1. A scheme that they proved to be stable uses Eqs. (3.3)–(3.5) to extrapolate the field corresponding to the outgoing characteristic to the boundary. Let us put a semi-Lagrangian twist to this idea.
For u > 0 we need a buffer zone for the region x < 0. Thus we extrapolate q0 from the interior using Eq. (3.5). For u < 0 we need a buffer zone for the region x > L. Thus we extrapolate piL from the interior using Eq. (3.4). To maintain stability we solve Eqs. (3.5) and (3.4) by finding the value of the field at the departure points x∗q = x0 − (u − c)Δt and x∗p = xiL − (u + c)Δt, respectively:
Say we are imposing ϕ at both boundaries and υ at inflow. For u > 0, υn+10 is imposed. Thus we know qn+1 at x0 from Eq. (3.34), meaning we can approximate (∂q/∂t)n and (∂2q/∂t2)n there. We can also estimate (∂p/∂t)n and (∂2p/∂t2)n there, using pn+1 = 2ϕn+1 − qn+1. In the same way, for u < 0, υn+1iL is imposed, and we can approximate (∂p/∂t)n and (∂2p/∂t2)n there, using Eq. (3.35). We can also estimate (∂q/∂t)n and (∂2q/∂t2)n there, using qn+1 = 2ϕn+1 − pn+1. Thus we can use Eqs. (3.4) and (3.5) to compute ∂p/∂x and ∂q/∂x:
e. Numerical testing (1): Advection
Let us perform numerical tests using a slow solution for which u = 0, and fυ = ∂ϕ/∂x, and the analytical solution transports these fields with a velocity u without changing their shape.
With this solution we mirror what we would be attempting with the full three-dimensional system of equations in an operational context. There, the analysis and initialization produces a well-balanced initial state; the slow initial state mirrors this. On the boundary we want the meteorological fields to enter the domain without corruption from the gravity waves; here we model this by imposing the slow solution on the boundary.
For ϕ we use the same bell curve as in section 2; thus
To crudely mirror the meteorological situation of a “front” entering the area through the western boundary and exiting through the eastern boundary we will begin the integration by placing the bell curve ϕ solution with its maximum on the western boundary. We will then impose boundary conditions consistent with the bell curve moving in the positive x direction with a velocity u.
We use L = 1000 km and Δx = 10 km. In addition, f = 10−4. In the tests that follow the length of the integration (T) will always be chosen such that in the analytical solution the bell curve moves exactly a distance L.
First, the forecast was run using trajectory truncation. That is, Eqs. (3.14)–(3.17) were solved with Ru, Rυ, and Rϕ computed as described in section 3b when the departure point was adjacent to the boundary. The result is shown in Fig. 7. The positive aspects of this forecast are (a) that the front has entered and exited the area reasonably accurately, (b) that there are no false reflections at the boundary, and (c) that there is no short-wave noise generated. The negative aspects are, of course, (a) the “tail,” which the ϕ solution has attained, and (b) the slowing down of the ϕ and υ solutions, just as in section 2.
Compare this with the forecast obtained using time interpolation on the boundary. That is, Eqs. (3.14)–(3.17) were solved when the departure point was within the boundary, and Eqs. (3.31)–(3.33) when the departure point was outside the boundary. The result is shown in Fig. 8. The positive aspects of the solution are retained and the negative ones eliminated. The solution is almost as good as the analytical one.
The well-posed buffer zone also gives stable forecasts, better than trajectory truncation. In the region x < 800 km the forecast is not quite as good as the time interpolation forecast; for x > 800 km it is marginally better; see Fig. 9.
f. Numerical testing (2): Adjustment
The analysis of Oliger and Sundström (1978) does not derive a unique set of fields required for well posedness. There are many choices. Obviously, “transparent” boundaries, that is, those that allow waves to pass through them without reflection, are the most desirable. In this section we demonstrate that imposing p/q on the boundary allows the gravity waves to pass through with little or no reflection, making it a more attractive option than imposing ϕ, which causes the gravity waves to be reflected.
Let us model the steady-state solution to the Rossby adjustment problem, described, for example, in Gill (1982, section 7.2.2). The initial state is given by
If we choose f = 10−4 and = 300 m s−1, then L = 30 000 km is 10 times the Rossby deformation radius . Thus, the geostrophically balanced large-scale flow at the boundaries will remain steady with time. We set u = 1 m s−1 so that we know we must impose υ at x = 0. Last, Δx = 100 km, Δt = 10 min, and ϕ̂ = 10. Because the boundary fields are constant and α is minuscule we will use trajectory truncation for both tests.
First let us examine what happens if we impose ϕ on both boundaries and υ at inflow:
The 10-day forecast is shown in Fig. 10. It is badly corrupted by the gravity waves, which have been reflected from the boundaries.
Next, let us impose our boundary values on the fields corresponding to the ingoing characteristics valid at the ϕ points, that is, p0 and qiL. How can we proceed? Let us assume that the following equations are valid at i = 0 and i = iL: (3.11), (3.13), and (3.14); also Eq. (3.16) at outflow, (i = iL). Now use the imposed boundary conditions to determine u−1/2 and uiL+1/2:
The trajectories are truncated as follows: Ru = (Xu)n1/2, when , when x∗ > L − Δx/2; Rϕ = (Xϕ)n0, when , when x∗ > L; Rυ = (Xυ)n0, when x∗ < 0; and , when x∗ > L.
We are now ready to do our second test. We impose υ and p at x = 0, and q at x = L,
When the departure point is outside the integration area an added complication is introduced into the discretized initial-boundary-value problem. We have seen that our three solutions to this complication, trajectory truncation, a well-posed buffer zone, and time interpolation, yielded stable integrations in practical tests. In this section we make a first attempt at confirming these results analytically by using the energy method to demonstrate the stability of the one-dimensional advection equation with constant velocity u > 0:
where 0 ⩽ x ⩽ L and t ≥ 0.
Using a linear interpolation to find the departure point value of the field the semi-Lagrangian discretization of Eq. (4.1) gives
Squaring both sides and rearranging the terms on the right-hand side results in
a. Trajectory truncation
This assumption implies that for the first p points Eq. (4.2) is replaced with
Therefore we get the following when we sum over the grid points:
Recall that the stability of the semi-Lagrangian scheme is dictated by the condition 0 ⩽ α̂ ⩽ 1. Therefore 0 ⩽ (α̂ − α̂2) ⩽ 1, enabling us to write
Using this equation again to substitute for ϕni on the right-hand side we get
Repeating this process we finally see that the solution is bounded by the initial condition and the boundary condition, thus establishing stability:
b. Time interpolation
For the first p points,
Squaring this, summing from 1 to p, and rearranging the terms we arrive at
Since we only “time interpolate” when α > p ≥ 1, we can write
Once the departure point is inside the integration area we switch to the standard space interpolation scheme, which means the remainder of the derivation is the same as above, resulting in
Thus the solution is bounded by the initial condition and the boundary condition, establishing stability.
c. Well-posed buffer zone
For the points outside the area we define the fields in terms of the imposed boundary condition:
Thus, in Eq. (4.3) the fields outside the area are defined in terms of the boundary condition and when we sum over all the grid points we get
Since (α̂ − α̂2) ⩽ 1 we can write
As before we can continue this procedure until we arrive at
Using Eq. (4.13), we finally see that the solution is bounded by the initial condition and the boundary condition, thus establishing stability, provided g is sufficiently smooth.
As we have seen, inventing stable boundaries even for simple one-dimensional systems is a nontrivial task. When we move to two (or three) dimensions the complications multiply. As is discussed in Durran (1999, section 8.2), there are fundamental difficulties encountered when we try to extend to two dimensions the idea of controlling the gravity waves as we did in section 3f. Nevertheless, Engquist and Majda (1977) have established a method constructing boundary conditions that are well posed and theoretically perfectly absorbing at least for waves striking at angles sufficiently close to normal incidence. For waves impinging at other angles of incidence (except tangential) reflections at the boundary can be reduced to a minimum by their method. The plan is to use their ideas and those discussed in this paper to construct a two-dimensional semi-implicit semi-Lagrangian shallow water model with boundaries that are “as transparent as possible.”
Extrapolations from the interior were found to be unavoidable. However, they did not lead to instabilities in our tests. This is consistent with the findings of Gustafsson et al. (1972), who proved that first-order extrapolations are stable for the Lax–Wendroff scheme (which is the same as a semi-Lagrangian scheme with quadratic interpolation and |α| < 1).
We have seen that both the time interpolation and well-posed buffer zone schemes overcome the accuracy limitation of trajectory truncation for large |α|. Looking ahead, it is interesting to ask what are the limitations of the former schemes in the context of an operational numerical weather prediction model? There is a limitation that applies to both; that is, their accuracy depends on how often the boundary is updated. For instance, if this is done every 3 h, neither scheme will be much more accurate than trajectory truncation. Each scheme has a particular limitation, at least as we have implemented them here. Convergence may be slow for an iterative implementation of the time interpolation scheme. Also, for our particular implementation of the well-posed buffer zone, looking at Eqs. (3.36) and (3.37) we see that it cannot be used when |u| = |c|, which can happen for internal gravity waves. Last, solving the two-dimensional equivalent of Eq. (3.17) will require the inversion of a sparse matrix. Thus any iterative procedures will be potentially expensive. This does not condemn either method but warns us that these preliminary implementations are not without possible weaknesses.
Although these are early days it is still nice to have a “fallback” position. One is to use trajectory truncation and accept the associated inaccuracy when |α| is large on the boundary. The second is based on an idea of J. E. Haugen (1999, personal communication). If uΔt/Δx < 0.5 on the western boundary, say, we need never invoke the time interpolation nor well-posed buffer zone schemes. To accomplish this we surround the uniformly discretized area with a buffer zone in which we use stretched coordinates, each Δx being, say, 15% bigger than its neighbor until this criterion is fulfilled for the maximum advecting velocity we expect perpendicular to the boundary. We choose this as our boundary and apply trajectory truncation there as an insurance against instability for unanticipated large inflow velocities.
Finally, although we have concentrated on looking at well-posed boundaries for semi-Lagrangian schemes in this paper, Eulerian systems are equally problematic. Consider an Eulerian approach using fourth-order spatial differencing. This will, in principle, require knowledge of fields outside the boundaries. Thus, stable extrapolations will be required to update the fields in the vicinity of the boundary. A possible solution is to adopt the upstream ideas inherent in the semi-Lagrangian approach. See also, however, Elvius and Sundström (1973) for an alternative approach.
Suggestions for improving the original manuscript by an anonymous reviewer and Dale Durran are gratefully acknowledged. In particular, the latter’s comment that the Rossby adjustment problem would provide an interesting challenge inspired section 3f.
Corresponding author address: Aidan McDonald, Met Éireann, Glasnevin Hill, Dublin 9, Ireland.