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

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 *u*_{0} and *υ*_{0}, and Φ_{0}, the mean geopotential height, are all constants.

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

which have solutions provided the following dispersion relation is valid:

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

and two adjustment waves with frequencies

and the solutions to Eq. (2.4) are

Here

The following combination of fields is a pure advection wave:

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

The group velocity of the advective waves is

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

### 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 ≤ *x* ≤ *L*_{x}; 0 ≤ *y* ≤ *L*_{y}. We will assume Φ_{0} > *u*^{2}_{0} + *υ*^{2}_{0}, 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* = *L*_{x}, 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

From our arguments in the previous paragraph this condition will hold for a large fraction of the waves in a meteorological environment, provided *l*^{2}/*k*^{2} does not become large. The ratio *l*^{2}/*k*^{2} is small for waves almost perpendicular to these boundaries.

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

The *x* = *L*_{x} boundary is an outflow boundary when *u*_{0} > 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 *l*^{2}/*k*^{2}, *f*^{2}/(*k*^{2}Φ_{0}), and *fl*/(*k*^{2}Φ_{0}), Eq. (2.14) tells us that

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 ∂Φ/∂*x* − *f**υ* − Φ_{0}(∂*u*/∂*x* + ∂*υ*/∂*y*) = 0. However, from Eqs. (2.1) and (2.3) we can rewrite it as

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*Φ̂^{h} − *f**υ̂*^{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 *l*^{2}/*k*^{2} ≪ 1 and *f*^{2}/(*k*^{2}Φ_{0}) ≪ 1. Therefore we take *υ*^{h} as one of the imposed fields. Consider the equation

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} + Φ_{0}*u*^{h}) such that Φ^{h} and *u*^{h} obey ∂Φ^{h}/∂*x* − *f**υ*^{h} = 0 and ∂*u*^{h}/∂*x* + ∂*υ*^{h}/∂*y* = 0. The third field, *u*, is extrapolated from the interior.

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

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

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

When *υ*_{0} < 0, then at *y* = *L*_{y} extrapolate *υ* from the interior and impose *u*^{h} and (Φ^{h} − Φ_{0}*υ*^{h}), making sure these host model fields are geostrophically balanced. At *y* = 0 extrapolate *u* and *υ* from the interior and impose

In what follows we will write ∂/∂*t* + *u*_{0}∂/∂*x* + *υ*_{0}∂/∂*y* as *d*/*dt*.

### c. Semi-Lagrangian discretization of the equations

We discretize such that *x*_{i} = *i*Δ*x*, *y*_{j} = *j*Δ*y*, and *t*_{n} = *n*Δ*t*; *x*_{I} = *L*_{x }*y*_{J} = *L*_{y}. A semi-Lagrangian and semi-implicit discretization of Eqs. (2.1)–(2.3) on an Arakawa C-grid that is *O*(Δ*t*^{2}) accurate is as follows:

where

where

where

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 *α* = *u*_{0}Δ*t*/Δ*x* and *β* = *υ*_{0}Δ*t*/Δ*y* then *Y*^{n}_{i,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 *Y*_{u}, *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:

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

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

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} > |**v**_{0}|. 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* = *L*_{x} we assume Eq. (2.2) is valid and use Eq. (2.35) to compute *υ* there. At *y* = 0 and *y* = *L*_{y} we assume Eq. (2.1) is valid and use Eq. (2.34) to compute *u* there.

As we discuss below, *Y*_{u} 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* = *L*_{y}

and at *x* = 0 and *x* = *L*_{x}

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 (*Y*_{u})^{n}_{i+(1/2),j,∗} = *Y*_{u}(1/2, *j* − *β*_{i+(1/2),j}, *n*Δ*t*) when *x*∗ < Δ*x*/2; (*Y*_{υ})^{n}_{i,j+(1/2),∗} = *Y*_{υ}(0, *j* + 1/2 − *β*_{i,j+(1/2)}, *n*Δ*t*) when *x*∗ < 0; and (*Y*_{Φ})^{n}_{i,j,∗} = *Y*_{Φ}(1, *j* − *β*_{i,j}, *n*Δ*t*) when *x*∗ < Δ*x*.

At *x* = *L*_{x} we truncate by using (*Y*_{u})^{n}_{i+(1/2),j,∗} = *Y*_{u}(*I* − 1/2, *j* − *β*_{i+(1/2),j}, *n*Δ*t*) when *x*∗ > (*I* − 1/2)Δ*x*; (*Y*_{υ})^{n}_{i,j+(1/2),∗} = *Y*_{υ}(*I*, *j* + 1/2 − *β*_{i,j+(1/2)}, *n*Δ*t*) when *x*∗ > *I*Δ*x*; and (*Y*_{Φ})^{n}_{i,j,∗} = *Y*_{Φ}(*I* − 1, *j* − *β*_{i,j}, *n*Δ*t*) when *x*∗ > (*I* − 1)Δ*x*.

At *y* = 0 we truncate by using (*Y*_{u})^{n}_{i+(1/2),j,∗} = *Y*_{u}(*i* + 1/2 − *α*_{i+(1/2),j}, 0, *n*Δ*t*) when *y*∗ < 0; (*Y*_{υ})^{n}_{i,j+(1/2),∗} = *Y*_{υ}(*i* − *α*_{i,j+(1/2)}, 1/2, *n*Δ*t*) when *y*∗ < Δ*y*/2; and (*Y*_{Φ})^{n}_{i,j,∗} = *Y*_{Φ}(*i* − *α*_{i,j}, 1, *n*Δ*t*) when *y*∗ < Δ*y*.

At *y* = *L*_{y} we truncate by using (*Y*_{u})^{n}_{i+(1/2),j,∗} = *Y*_{u}(*i* + 1/2 − *α*_{i+(1/2),j}, *J*, *n*Δ*t*) when *y*∗ > *J*Δ*y*; (*Y*_{υ})^{n}_{i,j+(1/2),∗} = *Y*_{υ}(*i* − *α*_{i,j+(1/2)}, *J* − 1/2, *n*Δ*t*) when *y*∗ > (*J* − 1/2)Δ*y*; and (*Y*_{Φ})^{n}_{i,j,∗} = *Y*_{Φ}(*i* − *α*_{i,j}, *J* − 1, *n*Δ*t*) when *y*∗ > (*J* − 1)Δ*y*.

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

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

*a*^{u}_{i+(1/2),j}= Δ*t*/(2Δ*x*) and*a*^{υ}_{i,j+(1/2)}= Δ*t*/(2Δ*y*).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 (*Y*_{u})_{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 (*Y*_{u})_{1/2,0} = *u*_{1/2,0}, and of course their equivalents in the other three corners. Similarly, for the balanced quantities in *Y*_{u} and *Y*_{υ}, whenever Φ is required on the boundary we use the following for options 2 or 3:

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 *L*_{x} = *L*_{y} = 10 000 km (there are 101 grid points in each direction). Below we model features of size *L*_{x}/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

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:

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.

#### 1) Φ imposed

We impose Φ(0, *y*, *t*) = Φ(*L*_{x}, *y*, *t*) = Φ(*x*, 0, *t*) = Φ(*x*, *L*_{y}, *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.

#### 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 *p*(Δ*x*/2, *y*, *t*) = *q*(*L*_{x} − Δ*x*/2, *y*, *t*) = *s*(*x*, Δ*y*/2, *t*) = *r*(*x*, *L*_{y} − Δ*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 (*u*_{0}, *υ*_{0}) starting from a position (*x*_{s}, *y*_{s}). 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:

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 (*x*_{s}, *y*_{s}) = (*L*_{x}/2, *L*_{y}/2) and (*u*_{0}, *υ*_{0}) = (50, 50), so that it leaves the area through the corner defined by *x* = *L*_{x}, *y* = *L*_{y}.

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*) = Φ(*L*_{x}, *y*, *t*) = Φ(*x*, 0, *t*) = Φ(*x*, *L*_{y}, *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 *L*_{x}/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.

#### 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 *p*(Δ*x*/2, *y*, *t*) = *q*(*L*_{x} − Δ*x*/2, *y*, *t*) = *s*(*x*, Δ*y*/2, *t*) = *r*(*x*, *L*_{y} − Δ*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).

#### 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 (*x*_{s}, *y*_{s}) = (5*L*_{x}/4, 5*L*_{y}/4) and choose (*u*_{0}, *υ*_{0}) = (−50, −50) so that it enters through the corner defined by *x* = *L*_{x}, *y* = *L*_{y}.

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

#### 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 *p*(Δ*x*/2, *y*, *t*), *q*(*L*_{x} − Δ*x*/2, *y*, *t*), *s*(*x*, Δ*y*/2, *t*), and *r*(*x*, *L*_{y} − Δ*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.

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

### 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: (*x*_{s}, *y*_{s}) = (3*L*_{x}/4, *L*_{y}/4). Since (*u*_{0}, *υ*_{0}) = (50, 50), the apex of the bell passes through the boundary at (*L*_{x}, *L*_{y}/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.

## 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 *t*_{0} 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

**,**

**,**

**,**

**,**

**,**

**,**

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

We impose *p* at *x* = Δ*x*/2:

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

where

Similarly, we impose *q* at *x* = *L*_{x} − Δ*x*/2:

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

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

At *y* = Δ*y*/2 and *y* = *L*_{y} − Δ*y*/2,

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

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 *u*_{0} > 0 we proceed as follows. If we are using option 2, *p*^{n+1}_{(1/2),j} = (*p*^{h})^{n+1}_{(1/2),j}, and *q*^{n+1}_{I−1/2,j} = (*q*^{h})^{n+1}_{I−(1/2),j}. If we are using option 3, *p*^{n+1}_{1/2,j} = (*p*^{h})^{n+1}_{1/2,j}, and *q*^{n+1}_{I−(1/2),j} = *q*^{n}_{I−(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