Abstract

Transparent boundary conditions for the linearized shallow water equations are constructed by incorporating the boundary conditions into equations that describe unidirectional waves. The shallow water equations are then discretized using a semi-Lagrangian approach and the transparency of the boundaries is demonstrated for three scenarios: adjustment waves radiating out of the area, a geostrophically balanced disturbance being advected in through the boundaries, and a geostrophically balanced disturbance being advected out of the area.

1. Introduction

There is a unique question associated with limited area modeling: what is the best way to treat the lateral boundaries? The answer can easily be stated: these boundaries should be “transparent.” That is, first, all waves approaching them from the interior of the limited area should exit without reflection; and second, in a nested environment, where the boundary data are being supplied by a “host” model, all meteorologically important waves impinging on the boundaries from the exterior of the limited area should enter without their amplitude or phase being changed and without exciting spurious high-frequency waves or noise.

Full transparency is probably unattainable in most practical integrations, but there exist a variety of methods, see Tsynkov (1998), for making the boundaries as transparent as possible. Of these, the method of Engquist and Majda (1977) is particularly appropriate for meteorological modeling; see Durran et al. (1993). The basic idea is that we incorporate the boundary conditions into equations that describe unidirectional waves traveling out of the area (and into the area when nested). In section 2 we will use this idea to arrive at one-way wave equations for the boundaries that are correct to a well-defined level of approximation.

Robert and Yakimiw (1986) point out that, since most boundary strategies contain difficulties that cannot be easily identified when they are considered in the framework of realistic models, it is a good idea to first test them on simple problems with known solutions. (If our boundary treatment will not work for a simplified subsystem it certainly will not work for the full nonlinear system.) That is the philosophy we adopt. Therefore, in section 4 we will take the linearized shallow water equations and compare the forecasts using our discretization and boundary treatment with the known solution for two cases. The first is an adjustment case for which an asymptotic solution is known. The second is an advection case whose analytical solution is known at all times.

2. The shallow water model

a. Solution to the linearized shallow water equations on an f-plane

The linearized shallow water equations are as follows:

 
formula

The x and y components of the winds are u and υ and Φ = gz is the geopotential height; g = 9.81 m s−2. The Coriolis parameter f, the advecting velocities u0 and υ0, and Φ0, the mean geopotential height, are all constants.

Substituting (u, υ, Φ) = (û, υ̂, Φ̂) exp[i(kx + lyωt)], gives rise to the equations

 
formula

which have solutions provided the following dispersion relation is valid:

 
formula

Thus there are three types of waves; advection waves with frequency

 
ωa = ku0 + lυ0,
(2.6)

and two adjustment waves with frequencies

 
formula

and the solutions to Eq. (2.4) are

 
formula

Here

 
formula

The following combination of fields is a pure advection wave:

 
formula

The following combinations of fields project out the two adjustment waves:

 
formula

The group velocity of the advective waves is

 
(Cga)x = u0;  (Cga)y = υ0,
(2.15)

and group velocities of the adjustment waves can be written as follows:

 
formula

b. Transparent boundary conditions

In this section we discuss the problem of designing transparent boundary conditions for the problem of modeling Eqs. (2.1)–(2.3) over the limited area 0 ≤ xLx; 0 ≤ yLy. We will assume Φ0 > u20 + υ20, which means we must impose two fields at inflow and one field at outflow; see Oliger and Sundström (1978). The arguments, however, can also be applied when this condition does not hold.

Consider the coefficient f/(KΦ0). We take f = 10−4 s−1. What would be a typical range of its values in the context of limited area meteorological modeling? For the external mode (Φ0 ∼ 300 m s−1) this coefficient is ≤0.05, assuming the longest wavelength to be ∼1000 km. The internal modes can have much smaller values of Φ0, of course. When the wavelength is 1000 km and Φ0 = 16 m s−1 then f/(KΦ0) ∼ 1. Therefore, if we are modelling over a 1000 km × 1000 km area f/(KΦ0) ≤ 1 for the majority of modes, and most importantly, f/(KΦ0) ≪ 1 for the fast-moving short wavelength modes that we particularly do not want to be reflected from the boundaries.

We examine first of all the boundaries at x = Lx, and x = 0. Looking at the group velocity in Eq. (2.15) we see that the energy of the advection waves propagates in the same direction as the advection velocity. From Eq. (2.16) we see that the energy of the ω± wave always propagates in the (±x) direction as long as

 
formula

From our arguments in the previous paragraph this condition will hold for a large fraction of the waves in a meteorological environment, provided l2/k2 does not become large. The ratio l2/k2 is small for waves almost perpendicular to these boundaries.

We will assume u0 > 0 in what follows. The arguments are quite general, however. We will subsequently state the results for u0 < 0 without repeating all the details.

The x = Lx boundary is an outflow boundary when u0 > 0. If we choose to extrapolate u and υ, how should we impose a third field? At this boundary there will be ωa wave packets and ω+ wave packets impinging from the interior. We would like our boundary treatment to facilitate the transmission of these waves without reflection and, of course, without exciting any ω waves. If we drop all the terms of order l2/k2, f2/(k2Φ0), and fl/(k2Φ0), Eq. (2.14) tells us that

 
formula

Thus, if we impose ikΦ̂Φ0(ikû + ilυ̂) − fυ̂ = 0 on this boundary we will be describing the outgoing ωa and ω+ waves correctly to this level of approximation. What about the incoming ω waves? From Eq. (2.14) we see that this condition imposes Φ̂ = 0, to this level of approximation, which is exactly what we want in a meteorological context where such a wave is unwanted “noise.”

Replacing ik and il with ∂/∂x and ∂/∂y respectively, this condition becomes ∂Φ/∂xfυΦ0(∂u/∂x + ∂υ/∂y) = 0. However, from Eqs. (2.1) and (2.3) we can rewrite it as

 
formula

which is probably easier to implement; certainly it is from a semi-Lagrangian point of view.

Next consider the boundary at x = 0. Here we must impose two fields and extrapolate the third. We would like to do this in such a way that the host model advective waves enter the area accurately without stimulating any ω+ waves, while simultaneously allowing the ω waves impinging from the interior to exit without reflection. Assume the host model fields are in geostrophic balance, that is, ikΦ̂hfυ̂h = 0, ilΦ̂h + fûh = 0, and ikûh + ilυ̂h = 0. (We are using the superscript h to designate the host model fields.) In these circumstances Eq. (2.12) tells us that the advective wave is dominated by υ̂h when l2/k2 ≪ 1 and f2/(k2Φ0) ≪ 1. Therefore we take υh as one of the imposed fields. Consider the equation

 
formula

From Eq. (2.13) we see that the ωa and ω waves obey this equation to the order of accuracy we are considering. Simultaneously it forces Φ̂+ to be zero; exactly what we want. Thus we must impose a second field on the boundary such that this equation holds. The simplest way to accomplish this is to impose (Φh + Φ0uh) such that Φh and uh obey ∂Φh/∂xfυh = 0 and ∂uh/∂x + ∂υh/∂y = 0. The third field, u, is extrapolated from the interior.

By the same argument we arrive at the conditions when u0 < 0. At x = Lx extrapolate u from the interior and impose υh and (ΦhΦ0uh), making sure these host model fields are geostrophically balanced. At x = 0 extrapolate u and υ from the interior and impose ∂Φ/∂xfυ + Φ0(∂u/∂x + ∂υ/∂y) = 0. Again, the most attractive way to impose this condition is probably to solve the following:

 
formula

The same logic leads to the following conditions at the other boundaries. When υ0 > 0, then at y = Ly extrapolate u and υ from the interior and use

 
formula

as the third condition. At y = 0 extrapolate υ from the interior and impose uh and (Φh + Φ0υh), making sure these host model fields are geostrophically balanced.

When υ0 < 0, then at y = Ly extrapolate υ from the interior and impose uh and (ΦhΦ0υh), making sure these host model fields are geostrophically balanced. At y = 0 extrapolate u and υ from the interior and impose

 
formula

In what follows we will write ∂/∂t + u0∂/∂x + υ0∂/∂y as d/dt.

c. Semi-Lagrangian discretization of the equations

We discretize such that xi = iΔx, yj = jΔy, and tn = nΔt; xI = Lx yJ = Ly. A semi-Lagrangian and semi-implicit discretization of Eqs. (2.1)–(2.3) on an Arakawa C-grid that is Ot2) accurate is as follows:

 
formula

where

 
formula

where

 
formula

where

 
formula

In these equations, the subscript i, j, * means evaluate this field by interpolating to the departure point associated with the grid point i, j. Defining α = u0Δtx and β = υ0Δty then Yni,j,∗ = Y(iα, jβ, nΔt). We use a Lagrangian bicubic interpolation to estimate Y at this point.

Equation (2.25) is valid for the indices i = 0, I − 1; j = 1, J − 1. Equation (2.28) is valid for the indices i = 1, I − 1; j = 0, J − 1. Equation (2.31) is valid for the indices i = 1, I − 1; j = 1, J − 1. We may need to estimate Yu, Yυ, or YΦ outside these ranges of indices. In that case, we will use carefully chosen approximations, as will be discussed in section 3.

In order to make our discussion of the boundary conditions in section 3 transparent let us solve Eqs. (2.25), (2.28), and (2.31) iteratively:

 
formula

d. Solving the equations

When we are imposing fields corresponding to the ingoing characteristics, then Eqs. (A.3), (A.6), (A.10), and (A.11) must be used in the vicinity of the boundary (see the appendix). This introduces a positional dependency to the coefficients multiplying Φ(k+1), u(k+1), and υ(k+1) in Eqs. (2.34)–(2.36). Taking this into account the set of equations we need to solve can formally be written as

 
formula

Substituting u from Eq. (2.37) and υ from Eq. (2.38) in Eq. (2.39) results in

 
formula

For the tests described in section 4 we use a successive overrelaxation (SOR) solver to solve Eq. (2.40) for Φ, and recover u and υ from Eqs. (2.37) and (2.38), respectively.

3. Discretizations at the boundaries

We impose the tangential velocity (υT) at inflow. We assume Φ0 > |v0|. Thus we must impose another field at inflow and one at outflow. We will consider three options below: 1) imposing Φ at all boundary points; 2) imposing the field corresponding to the ingoing characteristic, Φ − Φ0υN, as recommended by Elvius and Sundström (1973) (υN is the normal velocity); 3) imposing Φ − Φ0υN at inflow, and d(Φ − Φ0υN)/dt = 0 at outflow, as we discussed in section 2b. The normal velocities at both inflow and outflow emerge naturally when we solve the equations because of the C-grid staggering. [They are defined half a grid length in from the boundary and can be deduced from Eqs. (2.34) and (2.35).] There remain the tangential velocities at outflow. These we must extrapolate from the interior.

We update the tangential velocity at outflow as follows. At x = 0 and x = Lx we assume Eq. (2.2) is valid and use Eq. (2.35) to compute υ there. At y = 0 and y = Ly we assume Eq. (2.1) is valid and use Eq. (2.34) to compute u there.

As we discuss below, Yu and Yυ may be required on the boundary lines. We use different strategies depending on whether we are imposing Φ or Φ − Φ0υN on the boundary.

If the Coriolis terms are required on the boundary we use the following extrapolations. At y = 0 and y = Ly

 
formula

and at x = 0 and x = Lx

 
formula

The departure point can be beyond the boundary. We must devise a method for taking this into account. In this paper we will restrict ourselves to modest time steps. In that case trajectory truncation suffices. For large time steps it is necessary to use time interpolation; see McDonald (2001).

The trajectories are truncated as follows. At x = 0 we truncate by using (Yu)ni+(1/2),j,∗ = Yu(1/2, jβi+(1/2),j, nΔt) when x∗ < Δx/2; (Yυ)ni,j+(1/2),∗ = Yυ(0, j + 1/2 − βi,j+(1/2), nΔt) when x∗ < 0; and (YΦ)ni,j,∗ = YΦ(1, jβi,j, nΔt) when x∗ < Δx.

At x = Lx we truncate by using (Yu)ni+(1/2),j,∗ = Yu(I − 1/2, jβi+(1/2),j, nΔt) when x∗ > (I − 1/2)Δx; (Yυ)ni,j+(1/2),∗ = Yυ(I, j + 1/2 − βi,j+(1/2), nΔt) when x∗ > IΔx; and (YΦ)ni,j,∗ = YΦ(I − 1, jβi,j, nΔt) when x∗ > (I − 1)Δx.

At y = 0 we truncate by using (Yu)ni+(1/2),j,∗ = Yu(i + 1/2 − αi+(1/2),j, 0, nΔt) when y∗ < 0; (Yυ)ni,j+(1/2),∗ = Yυ(iαi,j+(1/2), 1/2, nΔt) when y∗ < Δy/2; and (YΦ)ni,j,∗ = YΦ(iαi,j, 1, nΔt) when y∗ < Δy.

At y = Ly we truncate by using (Yu)ni+(1/2),j,∗ = Yu(i + 1/2 − αi+(1/2),j, J, nΔt) when y∗ > JΔy; (Yυ)ni,j+(1/2),∗ = Yυ(iαi,j+(1/2), J − 1/2, nΔt) when y∗ > (J − 1/2)Δy; and (YΦ)ni,j,∗ = YΦ(iαi,j, J − 1, nΔt) when y∗ > (J − 1)Δy.

  1. Imposing Φ on the boundary fits quite naturally into the system of Eqs. (2.34)–(2.36). The only extrapolations required are those of Eqs. (3.1)–(3.4).

  2. To impose fields corresponding to the ingoing characteristics we proceed as above, except at the lines half a grid in from the boundary. There we proceed as described in the appendix with aui+(1/2),j = Δt/(2Δx) and aυi,j+(1/2) = Δt/(2Δy).

  3. Using d(Φ − Φ0υN)/dt = 0 at outflow is a minor variation on (2) in a semi-Lagrangian discretization; see the appendix.

Consider the corner defined by x = 0, y = 0. To compute (Yυ)0,1/2 and (Yu)1/2,0, we need Φ(0, 0), an unknown when we are using options 2 or 3. In the tests described in section 4 we use (Yυ)0,1/2 = υ0,1/2 and (Yu)1/2,0 = u1/2,0, and of course their equivalents in the other three corners. Similarly, for the balanced quantities in Yu and Yυ, whenever Φ is required on the boundary we use the following for options 2 or 3:

 
formula

where b + 1 means one grid point removed perpendicular to the boundary toward the interior. This is consistent with our arguments in section 2b.

Any other extrapolations from the interior to compute Φ on the boundary gave less accurate forecasts when the departure point was outside the area. When the departure point was inside the area, using Eqs. (A.2), (A.5), (A.8), and (A.9) to generate Φb was just as accurate in the above-mentioned tests as using Eq. (3.5).

4. Numerical testing

First, we examine an adjustment case whose asymptotic solution is known, and compare our longer forecasts with this asymptotic solution. Second, we examine an advective case that has an analytical solution with which we can compare our forecasts at all times.

For the following demonstrations Δx = Δy = 100 km, and Lx = Ly = 10 000 km (there are 101 grid points in each direction). Below we model features of size Lx/10. Thus we will generate waves with a typical wavelength of 1000 km or less. In all the following tests we use Φ0 = (5000 m)g, Φ = (500 m)g and f = 0.729 × 10−4 s−1. We use Δt = 15 min.

a. Testing adjustment

We would like to investigate the permeability of the boundaries to adjustment waves. Consider an initial state with ∇Φ(x, y, 0) ≠ 0 but with u(x, y, 0) = 0 and υ(x, y, 0) = 0. It will not be in geostrophic balance. Because of this, the system will radiate adjustment waves as it adjusts to a balanced state, which asymptotically is given by the field Ψ(x, y) that satisfies the equation

 
formula

See section 7.2.2 of Gill (1982) for details.

The rms differences between our extended forecasts and Ψ(x, y) will be an accurate measure of the effectiveness of our boundary treatment.

We will start from a bell shape at the center of the area:

 
formula

We also set the advecting velocities to zero. (In fact, both are assigned minuscule positive values in order to define inflow and outflow boundaries unambiguously.) The adjustment process consists of adjustment waves radiating away radially from the center of the area. We would like them to pass through the boundary without reflection.

The asymptotic balanced state arrived at by solving Eq. (4.1) with Φ(x, y, 0) given by Eq. (4.2) and Ψ = (5000 m)g on the boundary is not displayed but it is identical to the eye to the 48-h forecast shown in Fig. 2 below. In our tests below we will compare the 48-h forecasts with this balanced state to measure how transparent are the boundaries.

Fig. 2.

The adjustment of a bell shape using option 2: imposing Φ − Φ0υN on the boundaries

Fig. 2.

The adjustment of a bell shape using option 2: imposing Φ − Φ0υN on the boundaries

1) Φ imposed

We impose Φ(0, y, t) = Φ(Lx, y, t) = Φ(x, 0, t) = Φ(x, Ly, t) = (5000 m)g and υT = 0 on the inflow boundaries at all times. The forecasts for the first 12 h are displayed every 3 h in Fig. 1, as is the 48-h forecast. On the positive side, the forecast is stable (we have run it out to 10 days). So our boundary discretization is “experimentally well posed.” On the negative side, the adjustment waves are being strongly reflected from the boundaries. Measured against the asymptotic balanced state the errors for the 48-h forecast are huge. See column 2 of Table 1.

Fig. 1.

The adjustment of a bell shape using option 1: imposing Φ on the boundaries

Fig. 1.

The adjustment of a bell shape using option 1: imposing Φ on the boundaries

Table 1.

Rms difference between the 48-h forecast and the asymptotic solution given by Eq. (4.1). Column headings correspond to imposed conditions (see text)

Rms difference between the 48-h forecast and the asymptotic solution given by Eq. (4.1). Column headings correspond to imposed conditions (see text)
Rms difference between the 48-h forecast and the asymptotic solution given by Eq. (4.1). Column headings correspond to imposed conditions (see text)

2) Φ − Φ0υN imposed

In section 4a(1) we saw that imposing Φ produced stable forecasts but caused the adjustment waves to be reflected at the boundary. In this section we investigate to what extent imposing the fields associated with ingoing characteristics can reduce those reflections. Thus we are implementing the scheme described in the appendix and, as previously, imposing υT = 0 at the inflow boundaries. Now, however, instead of imposing Φ = (5000 m)g we impose px/2, y, t) = q(Lx − Δx/2, y, t) = s(x, Δy/2, t) = r(x, Ly − Δy/2, t) = (5000 m)g; (p, q, r, and s are defined in the appendix). The forecasts for the first 12 h are displayed every 3 h in Fig. 2, as is the 48-h forecast. The latter is identical to the eye to the asymptotic balanced state described by Eq. (4.1). Now, the boundaries are almost transparent to adjustment waves. The system is also experimentally well posed: the 10-day forecast shows no sign of instability. Measured against the asymptotic balanced state the errors are now small for the 48-h forecast. See column 3 of Table 1.

3) Φ − Φ0υN imposed at inflow and d(Φ − Φ0υN)/dt = 0 imposed at outflow

From Eqs. (2.8)–(2.9) we see that for adjustment waves this scheme is not more accurate than option 2). So this is just a test of stability. The forecasts are almost identical to those displayed in Fig. 2, which is confirmed by the rmse's, displayed in column 4 of Table 1.

b. Testing advection; bell exiting the area

In this section and the next we investigate the permeability of the boundaries to advective solutions. Do they enter and exit the area without generating any instabilities or unwanted adjustment waves? Do they enter and exit without refraction or reflection? Does option 3 reduce reflection at the outflow boundary as our analysis suggests?

The following is an analytical solution of Eqs. (2.1)–(2.3) that describes the advection of a bell shape with a constant velocity (u0, υ0) starting from a position (xs, ys). Thus we have an exact answer with which we can compare our integrations. Also, the system is in geostrophic balance and the analytical divergence is always zero, providing us with an additional useful test of the efficacy of our discretization:

 
formula

With our choices of the various parameters the bell can be thought of as a geostrophically balanced sharp pseudometeorological feature. (The maximum geostrophic wind is 57.5 m s−1.)

In this section we address the question, when the bell exits the area, how big are the reflections from the boundary? In particular, if the imposed boundary conditions are wrong, how badly behaved is the system? [This models an unavoidable (hopefully only occasional) operational situation in which the host model fields are inaccurate at outflow.] In order to make this a truly two-dimensional test let us start with the bell at the center of the area and advect it so that it exits through a corner; thus (xs, ys) = (Lx/2, Ly/2) and (u0, υ0) = (50, 50), so that it leaves the area through the corner defined by x = Lx, y = Ly.

In section 4c, we address the question, how accurately does the bell enter the area when we impose accurate boundary conditions?

1) Φ imposed

We impose Φ(0, y, t) = Φ(Lx, y, t) = Φ(x, 0, t) = Φ(x, Ly, t) = (5000 m)g. We impose υT = 0 on the inflow boundaries at all times. We extrapolate υT at outflow boundaries as explained in section 3. The initial state and the 18-, 24-, 30-, 36-, and 48-h forecasts are displayed in Fig. 3. As can be seen, the bell disappears almost without a trace. There is no sign of instability nor of two-grid noise. Looking carefully at the 36- and 48-h forecasts we can just see a feature with a dominant wavelength of about Lx/2. This feature is not in geostrophic balance (as can be seen from “x”s in the graph of mean absolute divergence shown in Fig. 5 below). The rms difference between the 48-h forecast and the analytical solution is shown in column 2 of Table 2.

Fig. 3.

The advection of a bell shape out of the area using option 1: imposing Φ on the boundaries

Fig. 3.

The advection of a bell shape out of the area using option 1: imposing Φ on the boundaries

Fig. 5.

Graph of the mean absolute divergence multiplied by 108. “x”s: Φ imposed; circles: Φ − Φ0υN imposed; and dots: Φ − Φ0υN imposed at inflow and d(Φ − Φ0υN)/dt = 0 imposed at outflow

Fig. 5.

Graph of the mean absolute divergence multiplied by 108. “x”s: Φ imposed; circles: Φ − Φ0υN imposed; and dots: Φ − Φ0υN imposed at inflow and d(Φ − Φ0υN)/dt = 0 imposed at outflow

Table 2.

Rms difference between the 48-h forecast and the analytical solution; bell exiting. Column headings correspond to imposed conditions (see text)

Rms difference between the 48-h forecast and the analytical solution; bell exiting. Column headings correspond to imposed conditions (see text)
Rms difference between the 48-h forecast and the analytical solution; bell exiting. Column headings correspond to imposed conditions (see text)

2) Φ − Φ0υN imposed

We have seen in section 4a that adjustment wave reflections can be reduced by imposing fields associated with ingoing characteristics. Will they have a similar beneficial effect in this advective scenario? In particular, will the geostrophic balance be maintained as the bell passes through the boundary? As before, we impose px/2, y, t) = q(Lx − Δx/2, y, t) = s(x, Δy/2, t) = r(x, Ly − Δy/2, t) = (5000 m)g. We impose υT = 0 at the inflow boundaries. When there is outflow, υT is computed by extrapolation as described in section 3.

We repeat experiment of section 4b(1) with these new boundary conditions. The forecasts are shown in Fig. 4. Looking at these charts we conclude that there is no instability and almost no reflection at the boundary. In fact, there is a small amount of error. The rms difference between the 48-h forecast and the analytical solution is shown in column 3 of Table 2. The errors have been reduced significantly, but there remains some reflection. The mean absolute divergence, shown as circles in Fig. 5, still shows a big increase as the bell reaches the boundary, but subsequent to that the geostrophic balance is much better than for section 4b(1).

Fig. 4.

The advection of a bell shape out of the area using option 2: imposing Φ − Φ0υN on the boundaries

Fig. 4.

The advection of a bell shape out of the area using option 2: imposing Φ − Φ0υN on the boundaries

3) Φ − Φ0υN imposed at inflow and d(Φ − Φ0υN)/dt = 0 imposed at outflow

We now impose Φ − Φ0υN = (5000 m)g at inflow and d(Φ − Φ0υN)/dt = 0 at outflow, treating υT in the same manner as always. Is the boundary is now transparent to the advective solution as we argued in section 2? We repeat the experiment of section 4b(2) with this new condition at outflow. The forecasts (not displayed) show greater transparency. This is demonstrated by the reduction in rmse (see column 4 of Table 2) and good geostrophic balance throughout the forecast; see the dots in the graph of mean absolute divergence, Fig. 5.

c. Testing advection; bell entering the area

In order to make this a truly two-dimensional test let us start with the bell at (xs, ys) = (5Lx/4, 5Ly/4) and choose (u0, υ0) = (−50, −50) so that it enters through the corner defined by x = Lx, y = Ly.

1) Φ imposed

We impose Φ on all boundaries and υT on the inflow boundaries at every time step so that these fields agree with the analytical description of the bell entering the area as we described above. The initial state and the 18-, 24-, 30-, 36-, and 48-h forecasts are displayed in Fig. 6. As can be seen, the bell enters the area stably but is accompanied by “trailing waves.” The stability of the 10-day forecast is confirmed by the plot of mean absolute divergence represented by the “x”s in Fig. 7, which also shows that geostrophic balance has still not been established after 10 days, long after the bell has exited the area. The rms difference between the 48-h forecast and the analytical solution is quite large; see column 2 of Table 3.

Fig. 6.

The advection of a bell shape into the area using option 1: imposing Φ on the boundaries

Fig. 6.

The advection of a bell shape into the area using option 1: imposing Φ on the boundaries

Fig. 7.

Graph of the mean absolute divergence multiplied by 108. “x”s: Φ imposed; dots: Φ − Φ0υN imposed at inflow; and d(Φ − Φ0υN)/dt = 0 imposed at outflow

Fig. 7.

Graph of the mean absolute divergence multiplied by 108. “x”s: Φ imposed; dots: Φ − Φ0υN imposed at inflow; and d(Φ − Φ0υN)/dt = 0 imposed at outflow

Table 3.

Rms difference between the 48-h forecast and the analytical solution; bell entering. Column headings correspond to imposed conditions (see text)

Rms difference between the 48-h forecast and the analytical solution; bell entering. Column headings correspond to imposed conditions (see text)
Rms difference between the 48-h forecast and the analytical solution; bell entering. Column headings correspond to imposed conditions (see text)

2) Φ − Φ0υN imposed

We know from section 4b(2) that imposing Φ − Φ0υN can improve the forecasts even for balanced solutions exiting the area. Will it have a similar beneficial effect with a balanced solution entering the area? In particular, will the geostrophic balance be maintained as the bell enters the area? Now we impose px/2, y, t), q(Lx − Δx/2, y, t), s(x, Δy/2, t), and r(x, Ly − Δy/2, t) and υT at the inflow boundaries at every time step such that these fields agree with the analytical fields for the bell entering as described above. We repeat the experiment of section 4c(1) with these new boundary conditions. The forecast is shown in Fig. 8. There are no trailing waves visible. The rmse for the 48-h forecast is significantly reduced in comparison with that of section 4c(1), see column 3 of Table 3. Also, the geostrophic balance has improved dramatically; the graph of mean absolute divergence (not shown to avoid clutter) coincides exactly with the dots in Fig. 7. The divergence is small throughout the 10 days. There is still a small increase in mean absolute divergence as the bell enters the area. The integration is still stable after 10 days.

Fig. 8.

The advection of a bell shape into the area using option 2: imposing Φ − Φ0υN on the boundaries

Fig. 8.

The advection of a bell shape into the area using option 2: imposing Φ − Φ0υN on the boundaries

3) Φ − Φ0υN imposed at inflow and d(Φ − Φ0υN)/dt = 0 imposed at outflow

We impose Φ − Φ0υN and υT at inflow as in section 4c(2) but now impose d(Φ − Φ0υN)/dt = 0 at outflow. Again this is mainly a test of stability, since all the activity is at inflow. The forecasts are almost identical to those of the experiment of section 4c(2), as can be seen from column 4 of Table 3 and Fig. 7, where it is inseparable from section 4c(2).

d. Time step dependence

For the tests described in section 4b the Courant number was less than 1 (α = β = 0.45). All these tests except that of section 4b(3) have been repeated in McDonald (2001) with α = β = 1.8 with satisfactory results. (Time interpolation is needed at inflow to maintain accuracy.)

It is important to emphasize that the test results described in this section do not depend on the semi-Lagrangian implementation. See McDonald (2001), where an Eulerian implementation is described and tested for every boundary condition except option 3.

To illustrate these two points let us repeat the “bell exiting” experiments described in section 4b. For the semi-Lagrangian discretization we use Δt = 60 min (α = β = 1.8) and for the three time level Eulerian discretization we use Δt = 10 min. For the latter the outflow boundary condition d(Φ − Φ0υN)/dt = 0 is implemented using the simplest possible upstream scheme; see, for example, Eq. (3.35) of Durran (1999). The results are displayed in Table 4.

Table 4.

Rms difference between the 48-h forecast and the analytical solution; bell exiting. SL refers to the semi-Lagrangian discretization; EU refers to the Eulerian discretization

Rms difference between the 48-h forecast and the analytical solution; bell exiting. SL refers to the semi-Lagrangian discretization; EU refers to the Eulerian discretization
Rms difference between the 48-h forecast and the analytical solution; bell exiting. SL refers to the semi-Lagrangian discretization; EU refers to the Eulerian discretization

Comparing Table 4 with Table 3 we see that for the semi-Lagrangian discretization increasing the time step is causing a small deterioration in accuracy for boundary condition 2 only. We also see that the Eulerian and semi-Lagrangian implementations are in reasonable agreement.

e. Bell exiting the area through the center of a boundary

One of the reviewers made the point that a common type of glancing condition can be inhibited in the corners. To address this the bell exiting test described in section 4b was repeated but now starting from a new position: (xs, ys) = (3Lx/4, Ly/4). Since (u0, υ0) = (50, 50), the apex of the bell passes through the boundary at (Lx, Ly/2) making a glancing angle of 45°. All other parameters are as described in section 4b. The results are displayed in Table 5, where we can now see that the relative improvement due to the option 3 boundary condition is more modest than previously.

Table 5.

Rms difference between the 48-h forecast and the analytical solution; bell exiting at (Lx, Ly/2). Column headings correspond to imposed conditions (see text)

Rms difference between the 48-h forecast and the analytical solution; bell exiting at (Lx, Ly/2). Column headings correspond to imposed conditions (see text)
Rms difference between the 48-h forecast and the analytical solution; bell exiting at (Lx, Ly/2). Column headings correspond to imposed conditions (see text)

5. Discussion

In this section the phrase “experimentally well posed” will be used to describe a discretization of the linearized shallow water equations that produced forecasts that had no visible two-grid noise on the 10-day chart and for which the mean absolute divergence was not increasing significantly with time at day 10.

The three sets of boundary conditions that we tested were experimentally well posed. They showed varying degrees of transparency in agreement with our analysis of section 2. It is perhaps worth emphasizing that well-posed boundaries of themselves, although attractive from a theoretical point of view, may not be particularly useful. They guarantee stable solutions but, as we have seen, they do not guarantee accurate solutions. When we imposed Φ on all boundaries we saw that adjustment waves were reflected at the boundary in our adjustment experiments, and that trailing waves developed in our advection experiments. When, instead, we imposed Φ − Φ0υN on all boundaries we obtained not only stable forecasts, but their accuracy was also improved dramatically. The adjustment wave reflection was reduced significantly, and the trailing waves phenomenon was eliminated. When we imposed Φ − Φ0υN at inflow boundaries and d/(Φ − Φ0υN)/dt = 0 at outflow we saw that the transmission of the advection solution improved significantly.

Although the rmse's are very small for the “bell entering” experiment there is a modest “ageostrophic shock” seen in Fig. 7. Perhaps the treatment of the corners is lacking in sophistication, or possibly we need to include higher-order terms than O(l/k) and O[f/(kΦ0)]. Looking at Eq. (2.12) we see that to all higher orders we would argue that the potential vorticity should be imposed at inflow. It is interesting that Charney (1960) proposed such a boundary condition for the barotropic primitive equations. However, Oliger and Sundström (1978, p. 441), have made a statement that gives us reason to pause: “[Charney's] approach has the liability that a small error committed initially or at any later time t0 will influence the boundary values at all later times. These errors will spread into the region of integration, contaminating the solution.”

We have demonstrated the technical feasibility of transparent boundary conditions for the linearized shallow water equations. An obvious next question to be addressed is, will these results hold up when the nonlinear terms, with their potential for explosive growth, are included? Assuming they do, can we apply these ideas directly to a multilevel model? Oliger and Sundström (1978) established necessary conditions for the well-posedness of the linearized hydrostatic equations by doing a normal mode decomposition. For each vertical mode this projects out a shallow water equation for which the well-posed boundaries are well known. (There is no guarantee that these will be sufficient for well-posedness.) This gives us an idea of how to proceed for hydrostatic models. Each mode has an associated Φ0, which enables us to apply the transparent boundary conditions presented in this paper to each of the projected shallow water equations. Once the boundaries have been updated we transform back to physical space. (For many semi-implicit models this concept is familiar. We solve the Helmholtz equation for the transformed fields.)

Acknowledgments

Thanks to two anonymous reviewers for their helpful criticism of the original manuscript, and to Jim Hamilton for assistance with the graphics.

REFERENCES

REFERENCES
Charney
,
J. G.
,
1960
:
Integration of the primitive and balance equations.
Proc. Int. Symp. on Numerical Weather Prediction, Tokyo, Japan, Meteorological Society of Japan, 131–152
.
Durran
,
D. R.
,
1999
:
Numerical Methods for Wave Equations in Geophysical Fluid Dynamics.
Springer-Verlag, 465 pp
.
Durran
,
D. R.
,
M-J.
Yang
,
D. N.
Slinn
, and
R. G.
Brown
,
1993
:
Toward more accurate wave-permeable boundary conditions.
Mon. Wea. Rev
,
121
,
604
620
.
Elvius
,
T.
, and
A.
Sundström
,
1973
:
Computationally efficient schemes and boundary conditions for a fine mesh barotropic model based on the shallow water equations.
Tellus
,
25
,
132
156
.
Engquist
,
B.
, and
A.
Majda
,
1977
:
Absorbing boundary conditions for the numerical simulation of waves.
Math. Comput
,
31
,
629
651
.
Gill
,
A. E.
,
1982
:
Atmosphere–Ocean Dynamics.
Academic Press, 662 pp
.
McDonald
,
A.
,
2001
:
Well-posed boundary conditions for semi-Lagrangian schemes: The two-dimensional case.
HIRLAM Tech. Rep. 47, 38 pp. [Available from A. McDonald, Met Éireann, Glasnevin Hill, Dublin 9, Ireland.]
.
Oliger
,
J.
, and
A.
Sundström
,
1978
:
Theoretical and practical aspects of some initial boundary value problems in fluid dynamics.
SIAM J. Appl. Math
,
35
,
419
446
.
Robert
,
A.
, and
E.
Yakimiw
,
1986
:
Identification and elimination of an inflow boundary computational solution in limited area model integrations.
Atmos.–Ocean
,
24
,
369
385
.
Tsynkov
,
S. V.
,
1998
:
Numerical solutions of problems on unbounded domains. A review.
Appl. Numer. Math
,
27
,
465
532
.

APPENDIX

Characteristic Boundaries

In this section we describe the procedure for solving the equations when the fields corresponding to the ingoing characteristics are imposed. To make the derivation general enough to fit all our various options consider the following generic u equation:

 
formula

We impose p at x = Δx/2:

 
formula

which we use to eliminate Φ(k+1)0,j from Eq. (A.1):

 
formula

where

 
formula

Similarly, we impose q at x = Lx − Δx/2:

 
formula

which when substituted in Eq. (A.1) results in

 
formula

We employ the same procedure to impose fields corresponding to the ingoing characteristics at the other two boundaries. Consider the generic υ equation:

 
formula

At y = Δy/2 and y = Ly − Δy/2,

 
formula

respectively. When substituted into Eq. (A.7) these give

 
formula

At all other grid points the u and υ equations are derived as in the body of the text and we proceed as described in section 2d.

For example, when u0 > 0 we proceed as follows. If we are using option 2, pn+1(1/2),j = (ph)n+1(1/2),j, and qn+1I−1/2,j = (qh)n+1I−(1/2),j. If we are using option 3, pn+11/2,j = (ph)n+11/2,j, and qn+1I−(1/2),j = qnI−(1/2),j,∗.

For an alternative approach see section 5.3.2 of Elvius and Sundström (1973), who impose Φ − Φ0υN on the boundary line itself instead of half a grid in.

Footnotes

Corresponding author address: Aidan McDonald, Met Éireann, Glasnevin Hill, Dublin 9, Ireland. Email: aidan.mcdonald@met.ie