A nonlinear time-dependent radiation condition for simulations of internal gravity waves in geophysical fluid flows

Applied Numerical Mathematics 110 (2016) 75–92

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A nonlinear time-dependent radiation condition for simulations of internal gravity waves in geophysical fluid flows V. Nijimbere, L.J. Campbell ∗ a r t i c l e

i n f o

Article history: Received 20 March 2016 Received in revised form 21 July 2016 Accepted 9 August 2016 Available online 12 August 2016 Keywords: Nonreflecting boundary condition Radiation condition Internal gravity waves

a b s t r a c t This paper examines the development of a time-dependent nonreflecting boundary condition (or radiation condition) for use in simulations of the propagation of internal gravity waves in a two-dimensional geophysical fluid flow configuration. First, a linear radiation condition, originally derived by Campbell and Maslowe, is implemented in some linear test cases. It involves the computation of a Laplace convolution integral which is nonlocal in time and thus requires values of the dependent variable at all previous time levels. An approximation for the integral is implemented here to reduce the expense of the computation and the results obtained are shown to be more accurate than those obtained using steady boundary conditions. For larger amplitude waves, nonlinear equations are required and the application of the linear radiation condition gives rise to instabilities. A new nonlinear time-dependent nonreflecting boundary condition is introduced which takes into account wave mean flow interactions in the vicinity of the outflow boundary by including a component corresponding to the vertical divergence of the horizontal momentum flux. This prevents the development of numerical instabilities and gives more accurate results in a nonlinear test problem than the results obtained using the linear radiation condition. © 2016 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Internal gravity waves in the atmosphere and ocean arise from the competing effects of the buoyancy and gravitation forces and play an important role in the transport of momentum and energy within the general circulation [17,22]. Gravity waves are relatively small-scale waves with wavelengths of up to 100 km, and are not generally resolved by large-scale general circulation models. However, in spite of their relatively small size, they contribute to large-scale phenomena that profoundly affect weather and climate, such as the quasi-biennial oscillation in the equatorial stratosphere, and it is thus important to represent their effects accurately in general circulation models. The drag forces that result from gravity wave interactions are included in these models by means of parameterization schemes which are based on information obtained from the mathematical theory of gravity waves and from numerical solutions of simplified models. Numerical and analytical studies of internal gravity waves are often based on the simplified two-dimensional configuration comprising a rectangular region defined in terms of horizontal and vertical coordinates. It is assumed that the waves are generated by a localized source representing a mechanism such as convection or topography [2,22]. Gravity waves gen-

*

Corresponding author. E-mail address: [email protected] (L.J. Campbell).

http://dx.doi.org/10.1016/j.apnum.2016.08.001 0168-9274/© 2016 IMACS. Published by Elsevier B.V. All rights reserved.

76

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

erated in the lower atmosphere generally propagate upwards into the stratosphere; thus, the source level of the waves can be considered to be the lower boundary of the rectangular domain and in theory the upper boundary should be at infinity. In practice, in numerical simulations, the waves are limited to a finite domain and a nonreflecting boundary condition may be needed at the upper boundary to prevent wave reflection and instabilities. Ideally, nonreflecting boundary conditions should take into account nonlinearity and time-dependence, but they are generally based on linearized steady equations. The issue of implementing an appropriate accurate nonreflecting boundary condition is not unique to internal gravity wave simulations but arises, in general, in numerical simulations of different types of waves in fluids, for example, planetary-scale Rossby waves which are generated in geophysical fluid flows by the effects of the Coriolis force [3]. This issue arises in simulations of wave propagation in other contexts besides fluids, e.g. for the scalar wave equation [11], for Maxwell’s equations [12] and for elastic wave equations [9]. A general technique called the artificial boundary method [13] has been applied to a variety of problems involving partial differential equations, some of which are described in the survey article of Han and Wu [14] and numerous references therein. In the geophysical fluid flow context, a number of different methodologies have been employed to simulate the propagation of waves in an infinite domain using a finite computational domain. One method, frequently used in atmospheric modeling, consists of surrounding the outflow boundary of the computational domain with a so-called sponge layer, a layer of high viscosity generally in the form of an additional term added to the model equations, with the goal of dissipating the energy of incident waves near the boundary. The viscosity in the sponge layer is generally represented using either a Rayleigh friction term of the form −γ ( z)ψ or a second-order viscous term of the form γ ( z)∇ 2 ψ , where in each case ψ is the unknown function representing the wave and γ is some specified function of the spatial variable z. The function γ ( z) is typically chosen to be either zero or very small in the interior of the domain and relatively-large near the outflow boundary. However, the addition of the artificial viscosity tends to introduce spurious nonphysical effects that affect the evolution of the waves in the rest of the computational domain. Another approach is to modify the flow characteristics in such a way as to create an evanescent region near the outflow boundary, thereby removing the waves from the computational domain near the boundary and preventing reflection. This idea was used by Geisler and Dickinson [8] for a simulation of Rossby waves and is an example of the buffer zone technique, variants of which have been used in a variety of different contexts in fluid dynamics over the past 50 years. The idea of using a radiation condition was described by Pearson [20]. The amplitude, wavelength and direction of propagation of the waves are specified at the outflow boundary by making use of wavenumber and group velocity information obtained from linear theory. For example, if the wave is defined by a function ψ which near the outflow boundary takes the ∂ψ form e ±imz , where m is the wavenumber in the z-direction, then a boundary condition of the form ∂ z = ±imψ is imposed. The sign of the wavenumber is chosen to give a wave with outward energy propagation. This type of implementation has the advantage of correctly representing the wave properties at the outflow boundary without damping or other nonphysical effects. However, it does not take into account the temporal evolution of waves. In the geophysical fluid context, a time-dependent radiation condition was first developed by Béland and Warn [3] for a configuration involving the propagation of Rossby waves on a two-dimensional horizontal beta-plane. The method involves linearizing the governing equations and making use of a Laplace transform to obtain a boundary condition in terms of a Laplace convolution integral of the time variable. Campbell [5] used this linear radiation condition for simulations of Rossby wave packets. A time-dependent radiation condition for gravity wave propagation also based on a Laplace convolution integral obtained from linearized equations was derived by Campbell and Maslowe [6]. One drawback of using this type of linear radiation condition, however, is that the Laplace convolution integral, in each case, is nonlocal in time and expensive to compute, requiring values of the dependent variable at all previous time levels. Also, when these linear time-dependent radiation conditions are applied to the original nonlinear problems, instabilities sometimes develop at late time as a result of having neglected the nonlinear effects near the outflow boundary. The numerical studies of Béland and Warn [3], Campbell and Maslowe [6] and Campbell [5], all dealt with configurations where there is a critical line in the flow. A critical line or a critical level is a position in space where the velocity of the wave is equal to the speed of the background shear flow and a critical layer is the region adjacent to the critical line. In both linear and nonlinear simulations, the waves are absorbed or reflected at the critical line and this reduces the extent of wave activity reaching the outflow boundary. In spite of this, Béland and Warn [3] noted some instabilities and observed that the rate of degradation of their simulations depended on both the intensity of the nonlinearities and the scale of motion near the computational boundary. Campbell and Maslowe [6] and Campbell [5] concluded that, with the level of spatial and temporal resolution used in their computations, it was enough simply to apply a zero boundary condition at the outflow boundary even in the nonlinear case, at least within the time frame required to examine the nonlinear critical layer phenomena that they were investigating. However, in a configuration where there is no critical line, the waves do reach the outflow boundary and the implementation of a radiation condition becomes essential. In particular, in a nonlinear simulation, nonlinear effects would need to be included near the outflow boundary. This paper re-visits the development of a time-dependent radiation condition for gravity wave simulations that was introduced by Campbell and Maslowe [6] and extend it to include nonlinearity in order to allow simulations involving larger amplitude waves. An approximation is proposed for the convolution integral so as to make the computations local in time and hence less expensive. Numerical simulations are then carried out using the time-dependent linear boundary condition in a configuration without a critical level where the waves are forced at the lower boundary of the rectangular domain and propagate without attenuation to the upper boundary of the domain. The results are compared with the results

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

77

obtained by using a steady radiation condition or by simply setting the fluid variables to zero at the outflow boundary. The special case of the long-wave approximation is examined, where the vertical-to-horizontal aspect ratio is taken to be zero in the governing equations. An exact analytical solution was derived by Nadon and Campbell [16] for this long-wave configuration. The accuracy of the time-dependent radiation condition is assessed by comparing the numerical results for this configuration with the analytical solution. It is seen that, with a steady radiation condition or a zero boundary condition, inaccurate results are obtained and instabilities develop, while the simulations with the time-dependent radiation condition more accurately represent the exact solution and are able to maintain their stability for a prolonged period of time. Following these preliminary linear experiments, we proceed with the primary goal of our investigation which is the development and implementation of a radiation condition for the nonlinear equations that can be used for simulations of larger amplitude waves. To this end, a new correction term based on weakly-nonlinear perturbation theory is added to the time-dependent radiation condition to include the effects of nonlinearity. This term represents the wave–mean-flow interactions near the outflow boundary and leads to a more accurate representation of the nonlinear evolution of the waves. 2. Nonlinear time-dependent radiation condition We develop and implement a nonlinear time-dependent boundary condition for internal gravity waves propagating in a density-stratified horizontal flow on a rectangular domain in a vertical plane defined by horizontal coordinate x and vertical coordinate z. The linear radiation condition of Campbell and Maslowe [6] is defined for a configuration based on the Boussinesq approximation [21]. The fluid density (x, z, t ) is approximated by the background density ρ¯ ( z), a function of z only, everywhere in the governing equations except in the term arising from the gravitational force, and the variation of dρ¯ the background density with altitude is considered to be small enough that we can neglect terms proportional to dz relative to the terms proportional to ρ¯ ( z). However, the nonlinear correction made here is not dependent on this approximation and can be applied to any configuration where weakly-nonlinear theory can be justified. Under the Boussinesq approximation, the energy conservation equation can be written in terms of the density and the mass continuity equation can be written in its incompressible form. The latter approximation allows us to define a streamfunction (x, z, t ) for the flow, in terms of the horizontal and vertical components of the fluid velocity u (x, z, t ) and w (x, z, t ), by

∂ = −u , ∂z

∂ = w. ∂x

(1)

In nondimensional form, the total streamfunction and the total density are both written as the sum of a background flow component and a relatively small-amplitude perturbation as

¯ z) + ε ψ(x, z, t ) (x, z, t ) = ψ(

(2)

(x, z, t ) = ρ¯ (z) + ερ (x, z, t ),

(3)

and

where ε is a small parameter. The background flow quantities ψ¯ ( z) and ρ¯ ( z) are, respectively, the horizontal mean values of ¯ dz (x, z, t ) and (x, z, t ) at the initial time t = 0. The horizontal component of the background flow velocity is u¯ = −dψ/ and the vertical component is zero. In this atmospheric context, the effect of viscosity and thermal conduction can be considered small and are neglected in the derivation of the radiation condition. However, these effects could potentially be included in the numerical simulations in the interior of the domain as done, for example, by Campbell and Maslowe [6]. Without viscosity and heat conduction, the nonlinear equations for the streamfunction and density perturbations are [6]

∇ 2 ψt + u¯ ∇ 2 ψx − u¯ zz ψx + g (ρ¯ )−1 ρx + ε (ψx ∇ 2 ψz − ψz ∇ 2 ψx ) = 0

(4)

ρt + u¯ ρx + ρ¯ z ψx + ε(ψx ρz − ψz ρx ) = 0,

(5)

and

where g is the acceleration due to gravity and the subscripts denote differentiation with respect to x, z and time t. In 2

2

these nondimensional equations, the Laplacian operator takes the form ∇ 2 = α 2 ∂∂x2 + ∂∂z2 , where α = L z / L x is the verticalto-horizontal aspect ratio. The vertical scale L z is of the order of the dimensional scale height of the atmosphere while the horizontal scale L x is of the order of the typical horizontal extent of the physical phenomena, such as topography and convection, that generate atmospheric gravity waves. It is reasonable therefore to consider L z to be much smaller than L x and thus consider α to be a small parameter in the nondimensional equations. Setting α to zero gives the long-wave limit and corresponds to the case of hydrostatic flow. In this special case, an analytical solution has been derived for the linearized time-dependent equations [16], which is used here validating the results of the numerical computations. Equations (4) and (5) are examined in a semi-infinite rectangular domain given by −∞ < x < ∞, z1 < z < ∞. At the lower boundary z = z1 , a boundary condition of the form

ψ(x, z1 , t ) = e ik(x−ct ) + c.c.

(with k and c being constants)

(6)

78

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

Fig. 1. A schematic diagram of the configuration described in Section 2.

is applied for t > 0. This is a linearized representation of a small-amplitude perturbation to the background flow streamfunction such as could arise from a forcing mechanism such as topography or convection in the lower atmosphere. The perturbation density and the perturbation vorticity and its z derivative are zero at this boundary z = z1 . The forcing (6) generates a gravity wave with horizontal wavenumber k and horizontal wavelength 2π /k, which propagates upwards into the domain. It is sufficient to consider a finite horizontal interval defined by x1 < x < x2 , with length |x2 − x1 | equal to a multiple of 2π /k, and apply periodic boundary conditions in the x direction. We require that there be no incoming energy from above, so that, in the absence of wave reflections, only waves with positive group velocity are present. This requirement serves as an upper boundary condition for the problem. We also impose the initial conditions that the perturbation density and the perturbation vorticity and its time derivative are zero at time t = 0. A schematic diagram of this initial–boundary-value problem is given in Fig. 1. g Applying the differential operator ∂∂t + u¯ ∂∂x to equation (4) and the differential operator ρ¯ ∂∂x to equation (5) and combining the resulting equations gives



∂ ∂ + u¯ ∂t ∂x

2

 ∇ 2 ψ − u¯ zz

   ∂ ∂ ∂ ∂ ψx + N 2 ψxx + ε + u¯ + u¯ J (ψ, ∇ 2 ψ) = 0, ∂t ∂x ∂t ∂x

(7)

where

J (ψ, ∇ 2 ψ) = ψx ∇ 2 ψz − ψz ∇ 2 ψx and N is the Brunt–Väisälä or buoyancy frequency, given by

N2 = −

g ∂ ρ¯ . ρ¯ ∂ z

(8)

In general, the background flow quantities u¯ and ρ¯ are continuous functions of z and, hence, N is as well. In order to solve this problem numerically, the semi-infinite interval z1 < z < ∞ must be replaced by a finite interval z1 < z < z2 with an appropriate nonreflecting boundary condition or radiation condition applied at the upper boundary z = z2 . The linear radiation condition of Campbell and Maslowe [6] requires that u¯ and N both be constant, but it can be applied in an approximate sense to situations where either of these quantities is a slowly-varying function of z near the upper boundary of the computational domain. For example, Campbell and Maslowe [6] considered the case where u¯ ( z) = tanh( z − 5), 0 < z < 10; this is a hyperbolic tangent profile centered at z = 5 and with u¯ ≈ 1 near the upper boundary z = 10. In the current investigations, we set both u¯ and N to constants. This means that the background density, according to (8), takes the form ρ¯ ( z) = e −z/ H 0 , where H 0 is the scale height of the atmosphere, and N 2 = g / H 0 . With u¯ constant, u¯ zz is zero and so the second term in (7) is omitted. For small-amplitude perturbations, a linear approximation can be obtained by setting ε to zero in (7). The radiation condition of Campbell and Maslowe [6] is based on this linear model. Here ε is allowed to be nonzero but small, and ψ is written as a perturbation series,

ψ ∼ ψ (0) + ε ψ (1) + O (ε 2 ),

(9)

where ψ (0) is the solution of the linearized equation. Substituting this into (7) and grouping terms according to the order of

ε gives the linear equation

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92



∂ ∂ + u¯ ∂t ∂x

2

( 0)

∇ 2 ψ (0) + N 2 ψxx = 0

79

(10)

at leading order. At this order, the solution satisfies the specified boundary condition (6) and thus takes the form

ˆ z, t )e ikx + c.c. ψ (0) (x, z, t ) = ψ(

(11)

ˆ z, t ) complex. For a more general configuration where the function specified at the boundary defines a wave packet with ψ( comprising a discrete spectrum of horizontal wavenumbers, the solution can be written as a complex Fourier series ∞ 

ψ (0) (x, z, t ) =

ψˆ (0) (κ , z, t )e i κ x .

(12)

κ =−∞

The linear time-dependent radiation condition of Campbell and Maslowe [6] is obtained by substituting (11) or (12) into (10), taking the Laplace transform of each term and applying the zero initial conditions. This gives ( 0) − ψ˜ zz



N 2κ 2

(s + i κ u¯ )2

 + α 2 κ 2 ψ˜ (0) = 0,

(13)

˜ κ , z, s) is the Laplace transform of ψ( ˆ κ , z, t ). This is the Taylor–Goldstein equation [10,23] for a wave with complex where ψ( phase speed c = is/κ and can be written as



  ∂ ∂ − H (κ , s) + H (κ , s) ψ˜ (0) (κ , z, s) = 0, ∂z ∂z

(14)

where

 H (κ , s) = ακ

(s + i κ u¯ )2 + N 2 /α 2 (s + i κ u¯ )2

1/2 .

(15)

If H (κ , s) is defined by choosing the branch of the square root in (15) so that Re[ H (κ , s)] > 0, then the solution of (14) which is bounded as z → ∞ is the one that satisfies ( 0) ψ˜ z + H (κ , s)ψ˜ (0) = 0.

(16)

The inverse Laplace transform of H (κ , s) is [1]

h(κ , t ) =

N 2κ

t

α

a

τ

¯ ¯ J1 (aτ )dτ e −i κ ut + ακ e −i κ ut δ(t − 0),

(17)

0

where J1 is the Bessel function of the first kind of order 1, δ(t − 0) is the Dirac delta function and a(κ ) = N 2 κ /α 2 . Thus, ψˆ (0) (κ , z, t ) satisfies the radiation condition

∂ ( 0) ψˆ (κ , z, t ) + ακ ψˆ (0) (κ , z, t ) = − ∂z

t

ψˆ (0) (κ , z, τ ) g (κ , t − τ )dτ ,

(18)

0

with

g (κ , t ) = ακ e

¯ −i κ ut

t

a

η

J1 (aη)dη.

(19)

0

This is the linear time-dependent radiation condition of [6]. In numerical computations, it is applied at the upper boundary of the computational domain (z = z2 ). For the special case of long waves where the aspect ratio α  1, the approximations

H (κ , s) ≈

Nκ s + i κ u¯

and ¯ h(κ , t ) ≈ N κ e −i κ ut

lead to an approximate linear radiation condition

80

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

∂ ( 0) ¯ ψˆ (κ , z, t ) = − N κ e −i κ ut ∂z

t

ψˆ (0) (κ , z, τ )e i κ u¯ τ dτ ,

(20)

0

which is applied at z = z2 . The O (ε ) terms in the perturbation equation give a nonhomogeneous equation



∂ ∂ + u¯ ∂t ∂x

2

 (1 ) ∇ 2 ψ (1) + N 2 ψxx =−

 ∂ ∂ + u¯ r (x, z, t ), ∂t ∂x

(21)

where

r (x, z, t ) = J (ψ (0) , ∇ 2 ψ (0) ), and the solution ψ (1) can be written as a Fourier series

ψ (1) (x, z, t ) =

∞ 

ψˆ (1) (κ , z, t )e i κ x .

(22)

κ =−∞

Taking the Fourier–Laplace transform of (21) gives



  r˜ (κ , z, s) ∂ ∂ − H (κ , s) + H (κ , s) ψ˜ (1) (κ , z, s) = − , ∂z ∂z s + i κ u¯

(23)

where ψ˜ (1) (κ , z, s) and r˜ (κ , z, s) are, respectively, the Fourier–Laplace transforms of ψ (1) (x, z, t ) and r (x, z, t ). With Re[ H (κ , s)] > 0, the solution of (23) which is bounded as z → ∞ is the one that satisfies



 z ∂ r˜ (κ , ξ, s) (1 ) H (κ , s ) z ˜ + H (κ , s) ψ (κ , z, s) = −e e − H (κ ,s)ξ dξ. ∂z s + i κ u¯

(24)

z1

With the boundary condition (6), the only nonzero components in the Fourier spectrum of ψ (0) correspond to the wavenumber ±k and they give rise to O (ε ) components corresponding to the wavenumbers zero and ±2k. For κ = ±2k, there is no exact formulation for the inverse Laplace transform of the expression on the right-hand side of (24), although a numerical approximation is possible. For the zero wavenumber, H (0, s) = 0 and so equation (24) is simply

∂ (1 ) ψ˜ (0, z, s) = − ∂z

z

r˜ (0, ξ, s) s

(25)

dξ.

z1

Thus, the zero wavenumber component of the solution satisfies the radiation condition

∂ (1 ) ψˆ (0, z, t ) = − ∂z

t z rˆ (0, ξ, τ )dξ dτ .

(26)

0 z1

This result could of course have been obtained directly without the use of a Laplace transform, by setting integrating (21) with respect to z and t. Noting that (Appendix B)

rˆ (0, z, t ) =

k

2π /k

κ r (x, z, τ )dx = 2π

2π 0

κ = 0 in (22) and

2π /κ

J (ψ (0) , ∇ 2 ψ (0) )dx 0

=

k

2π /k

2π

∂2 ( 0) ( 0) (ψx ψz )dx, ∂ z2

(27)

0

we observe that the expression on the right-hand side of (26) corresponds to the leading-order component of the horizontalmean momentum flux divergence at the outflow boundary. The total change in the horizontal-mean horizontal velocity component as a result of the momentum flux divergence is u¯ 0 ( z, t ) = −ε ∂∂z ψˆ (1) (0, z, t ) and equation (26) can be written as

t u¯ 0 ( z, t ) = −ε

2 0

∂F ( z, τ )dτ ∂z

(28)

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

81

or

∂ u¯ 0 ∂F ( z, t ) = −ε 2 ( z, t ), ∂t ∂z

(29)

where

F ( z, t ) = −

k

2π /k

( 0)

( 0)

ψx ψz dx

2π

(30)

0

is the leading-order horizontal-mean momentum flux or Reynolds stress. This is the usual equation for the acceleration of the mean flow that results from nonlinear gravity wave interactions. Since N and u¯ are constant, the leading-order linear equation (10) has a solution with steady (constant) amplitude of the form

ψ(x, z, t ) = Ae ikx+imz−i ωt + c.c.,

(31)

where k, m and ω satisfy the gravity wave dispersion relation (A.6). Substituting the expression (31) into the integral in (30) gives F equal to a constant. This result is the Eliassen–Palm theorem of geophysical fluid dynamics [7]; it states that for steady-amplitude gravity waves satisfying the linear equation (10), the horizontal-mean momentum flux or Reynolds stress is independent of altitude. Thus, if the leading-order solution takes the form (31), then ∂ F /∂ z is zero and so, according to (29), there is no mean flow acceleration at O (ε ). With appropriate choices of the constants A, m and ω , the steady-amplitude solution (31) can be made to satisfy the boundary condition (6) at z = z1 , as well as the condition of positive group velocity. Applying these conditions to (31) gives

ψ(x, z, t ) = e ik(x−ct )+im(z−z1 ) + c.c.,

(32)

where m is the vertical wavenumber and satisfies

 m=

N2

(u¯ − c )2

1/2 − α 2 k2

.

(33)

Following the discussion in Appendix A, m is defined to be positive (negative) in order to have a solution with upward group velocity for k > 0 and (u¯ − c ) > 0 ((u¯ − c ) < 0). However, the steady-amplitude solution (32) does not satisfy the zero initial conditions specified in the time-dependent initial–boundary-value problem. The solution of the linear time-dependent problem takes the form (11) with an ampliˆ z, t )| that depends on both time and altitude. An asymptotic approximation for this solution, valid for large t, tude |ψ( can be obtained by taking a Laplace transform of the linear time-dependent equation (10), taking into account the zero initial conditions, and then applying the boundary condition (6) at z = z1 , as well as the condition of positive group velocity. For the special case of zero aspect ratio (the long-wave limit), an exact analytical solution was obtained by Nadon and Campbell [16]. The solution can be expressed in two different but equivalent forms as infinite series involving Bessel functions [19]. One form is

ψ(x, z, t ) = e ik(x−ct ) e −e

iN ( z− z1 ) (u¯ −c )

−ik(x−(u¯ −c )t )

√ n  ∞  i N ( z − z1 ) J n 2 Nk( z − z1 )t , √ (u¯ − c ) kt n =1

(34)

where the functions J n are Bessel functions of order n. As t → ∞, all the terms in the infinite series approach zero and so this time-dependent solution approaches the steady-amplitude solution (32) (with α = 0). In the general case of nonzero aspect ratio as well, the solution of the linear time-dependent initial–boundary value problem approaches the steadyamplitude solution (32) as t → ∞. For finite time, the leading-order momentum flux (30) also depends on time and altitude and thus, in the nonlinear problem, there is an O (ε ) acceleration of the mean flow velocity, according to (29). The radiation condition (26) for the zero-wavenumber component of the time-dependent solution can be written as

∂ (1 ) ψˆ (0, z, t ) = ∂z

t

∂F ( z, τ )dτ . ∂z

(35)

0

Equations (18) and (35) can be paired for use as a nonlinear time-dependent radiation condition in numerical simulations involving the nonlinear time-dependent problem given by equations (4) and (5) with the specified initial and boundary conditions. The implementation of the radiation condition and some numerical experiments are described in sections 3 and 4.

82

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

3. Numerical implementation of the radiation condition For the numerical simulations, equation (4) is written as

ζt + u¯ ζx − u¯ zz ψx + g (ρ¯ )−1 ρx + ε (ψx ζz − ψz ζx ) = 0,

(36)

ζ = ∇2ψ

(37)

where

is the vorticity perturbation. Equations (36) and (5) are solved in a rectangular domain defined by x1 < x < x2 , z1 < z < z2 (Fig. 1) subject to the boundary and initial conditions given in section 2. The gravity wave perturbation is generated at z = z1 and a radiation condition, either linear or nonlinear, is implemented at z = z2 . A Fourier spectral representation is employed in the x-direction to take advantage of the horizontally-periodic structure of the problem [6]. For the discretization in time and in the vertical direction, previous studies by Campbell and Maslowe [6] and references therein, and McLaren, Campbell and Vaillancourt [15] describe various methods tailored to the solution of related problems involving configurations with shear, i.e., with u¯ being a function z. In the present problem where there is no shear, the discretization is more straightforward and any of a number of standard methods could be used. For the results shown here, a predictor–corrector method is used with the explicit second-order Adams–Bashforth scheme for the predictor step and the implicit second-order Adams–Moulton scheme for the corrector step. This method is stable and allows the use of relatively large time increments with small errors (section 4). In the vertical direction, a second-order finite difference approximation is used. The interval z1 < z < z2 is divided into subintervals of uniform size  z. The nonlinear terms are evaluated using a pseudo-spectral method; this means that the nonlinear terms are computed by inverting the Fourier transform numerically at each time level. Once ζ is obtained, equation (37) is solved for ψ using an elliptic solver and then (36) and (5) are used to compute ζ and ρ at the next time level. This numerical implementation procedure is described in detail in [18]. In implementing the linear time-dependent radiation condition, we note that the convolution integral in (18) is nonlocal in time; computing this integral at a given time level t = tn requires storing and summing up information from all previous levels. A less expensive method of computation is to make the discretization local in time by expressing the integral from t = 0 to t = tn as a sum of two integrals t n −1  tn ˆ ˆ κ , z, τ ) g (κ , tn − τ )dτ . Bn = − ψ(κ , z, τ ) g (κ , tn − τ )dτ − ψ(

(38)

t n −1

0

By approximating the first integral using the value of g known from the time level tn−1 , we obtain the approximation t n −1  tn ˆ ˆ κ , z, τ ) g (κ , tn − τ )dτ = B n−1 +  B n . Bn ≈ − ψ(κ , z, τ ) g (κ , tn−1 − τ )dτ − ψ(

(39)

t n −1

0

This allows us to simply compute an increment  B n and add it to B n−1 at each time level tn−1 to obtain an approximation for B n at the next time level tn . Starting with B 0 = 0 at t = 0, we approximate  B n , for n = 1, 2, . . . using the trapezoidal rule by

Bn ≈ −

t ˆ ˆ κ , z, tn ) g (κ , tn − tn )). (ψ(κ , z, tn−1 ) g (κ , tn − tn−1 ) + ψ( 2

Since g (κ , 0) = 0 from (19),

Bn ≈ −

t 2

ˆ κ , z, tn−1 ) g (κ , t ) = −bακ ψ(

t 2

ˆ κ , z, tn−1 )e −i κ u¯ t . ψ(

(40)

The real constant b is evaluated using the first few terms of the ascending series for the Bessel function J1 [1] as

t b=

a

σ 0

t J1 (aσ )dσ ∼

a2 2

 1−

a2 σ 2 8

+

a4 σ 4 192

   a2 a2 a4 + O (a6 σ 6 ) dσ ∼ 1− t 3 + t 5 + O (t 7 ) . 2

24

960

(41)

0

Equation (18) is discretized using the one-sided finite-difference approximation

(1 + ακ  z)ψˆ κn, jmax − ψˆ κn, jmax −1 =  zB n ,

(42)

ˆ κ , j  z, tn ) for j = jmax − 1, jmax (where jmax is the index of the grid point at the upper boundary). where ψˆ κn, j ≈ ψ( Since the computation of the integral  B n is only required once at each time level, the algorithm is fast and easy to implement. A small error is introduced at each time step due to the approximation of g (κ , tn − τ ) by g (κ , tn−1 − τ ) in (39).

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

83

Table 1 A list of the different configurations shown in the graphs in Section 4. Configuration

Upper boundary condition

Analytical solution (34): Linear with α = 0 (long-wave limit)

Positive group velocity for all z: z1 < z < ∞

Linear with α = 0 (long-wave limit)

Time-dependent radiation condition (20) at z = z2 Steady radiation condition (A.8) with (A.9) at z = z2 Zero boundary condition ψ = 0 at z = z2

Linear with α = 0 (nonzero aspect ratio)

Time-dependent radiation condition (18) at z = z2 Steady radiation condition (A.8) with (A.4) at z = z2 Zero boundary condition ψ = 0 at z = z2

Nonlinear with α = 0 (nonzero aspect ratio)

Time-dependent radiation condition (18) at z = z2 Time-dependent radiation condition (18) with nonlinear correction (35) at z = z2

The error that accumulates from t = 0 to t = tn is O (t  z tn−1 ). With the choice of t = 10−2 and  z = 10−2 used to obtain the results shown in section 4, the error is O (10−2 ) by t = 100. For the case where α = 0, both integrals in (38) can be evaluated directly, so the relation in (39) is exact and thus there is no error. For nonlinear simulations, the nonlinear correction (35) derived in section 2 is implemented. It involves evaluating the integral with respect to time of the vertical gradient of the horizontal-mean momentum flux at the outflow boundary and using this to evaluate the zero wavenumber component of the solution at the boundary. Since the integral in (35) is nonlocal in time, it is approximated as

tn Cn =

∂F ( z, τ )dτ = ∂z

0

t n −1 

∂F ( z, τ )dτ + ∂z

tn

∂F ( z, τ )dτ = C n−1 + C n ∂z

(43)

t n −1

0

and the increment C n is evaluated at each time level tn , starting with C 0 = 0 at t = 0. 4. Results of the numerical simulations The results of some representative linear and nonlinear simulations with the linear radiation condition (18) and the nonlinear correction (35) are presented here and compared with an analytical solution and with numerical results obtained using a zero boundary condition and a steady radiation condition (Appendix A). The simulations are carried out on a rectangular domain given by 0 ≤ x ≤ 2π and 0 ≤ z ≤ 15. The wave source is at z1 = 0 and is given by the boundary condition

ψ(x, z1 , t ) = cos kx =

e ikx 2

+ c.c.

(44)

This corresponds to the case where the phase speed of the forced wave c has been set to zero in (6). We note that choosing a nonzero phase speed c would give solutions of the same form as those presented here, but with the dependence on x replaced by dependence on (x − ct ). The horizontal wavenumber at the source is set to k = 2 and the √ mean flow speed is kept constant, u¯ = 0.5. The scale height is set to H 0 = 4.9 so that the Brunt–Väisälä frequency N = 2. In all the results shown here, the viscosity and the thermal diffusivity are set to zero. The various configurations presented here are summarized in Table 1. In the special case where the long-wave limit is taken in the linear problem (α = 0, ε = 0), the analytical solution (34) was obtained by Nadon and Campbell [16]. Figs. 2–5 show numerical results obtained for this case with different upper boundary conditions and compared with the analytical solution (34). Contour plots of the analytical solution and of the numerical solution obtained with the linear time-dependent radiation condition (20) are presented in Fig. 2. The corresponding results obtained with a zero upper boundary condition and with a steady radiation condition (Appendix A) ˆ z, t )] = ψ(0, z, t ) is shown in Fig. 4 are shown in Fig. 3. The real part of the complex amplitude of the streamfunction Re[ψ( as a function of altitude z at time t = 20 for the analytical solution and for each of the three different boundary conditions. In each of case, we also calculate the relative error of the computation, i.e. the absolute value of the difference between the numerical and analytical solution at each value of t and z divided by the maximum absolute value of the analytical solution over the vertical extent of the computational domain, 0 ≤ z ≤ 15. The absolute value of the analytical solution attains a maximum value of 1 at the lower boundary z = z1 = 0; thus, the relative error of each computation is equal to the absolute error. For each boundary condition, the relative error E ( z, t ) is shown in Fig. 5 as a function of time t at the upper boundary z = z2 = 15. Figs. 2–5 show that there is close agreement between the analytical solution and the numerical solution obtained with the time-dependent radiation condition (20). The relative error is shown by the solid curve in Fig. 5; it is barely visible

84

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

Fig. 2. Linear simulation of gravity wave propagation in the long-wave limit (ε = 0, α = 0): contour plots of the wave streamfunction ψ(x, z, t ). Figures (a), (b) and (c) show the analytical solution [16] at times t = 5, 10, 20, respectively. Figures (d), (e) and (f) show the numerical solution at times t = 5, 10, 20, respectively, obtained by applying the time-dependent radiation condition (20) at the upper boundary.

because it is close to the t axis and is much smaller than the relative errors for the other two boundary conditions. At t = 10, it is of the order of 10−3 . However, there are significant qualitative and quantitative differences between the analytical solution and the results obtained with the zero boundary condition and with the steady radiation condition. With the zero boundary condition, the waves are prevented from propagating beyond the upper boundary and this results in an incorrect wave structure, particularly in the upper two-thirds of the computational domain. With the steady radiation condition (Appendix A), the direction of propagation of the waves is artificially fixed at the upper boundary and the wave amplitude grows with time in the interior of the computational domain. The time-dependent radiation condition, on the

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

85

Fig. 3. Linear simulation of gravity wave propagation in the long-wave limit (ε = 0, α = 0): contour plots of the wave streamfunction ψ(x, z, t ). Figures (a), (b) and (c) show the numerical solution at times t = 5, 10, 20, respectively, obtained with a zero upper boundary condition. Figures (d), (e) and (f) show the numerical solution at times t = 5, 10, 20, respectively, obtained by applying the steady radiation condition (A.8) with (A.9) at the upper boundary.

other hand, allows the direction of propagation of waves to change correctly with time and gives an accurate representation of the wave amplitude at all altitudes. √ Fig. 6 shows results obtained for the general linear configuration with a nonzero aspect ratio, set to α = 0.4, with the other input parameters being the same as in Figs. 2–5. Contour plots of the streamfunction ψ(x, z, t ) at times t = 10, 20, 45 are shown for a simulation with the linear time-dependent radiation condition (18). In this case, we do not have an exact time-dependent solution with which to compare the numerical solution. However, from the discussion in section 2, we expect the numerical solution to approach the steady-amplitude solution (32) as t → ∞. This is seen in Fig. 7 which shows the numerical solution at time t = 45, along with the exact steady-amplitude solution (32).

86

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

Fig. 4. Linear simulation of gravity wave propagation in the long-wave limit (ε = 0, α = 0) with different boundary conditions at the outflow (upper) ˆ z, t )] = ψ(0, z, t ) as a function of altitude z at time t = 15. The solid curve shows the analytical boundary: the real part of the complex wave amplitude Re[ψ( solution [16], and the dashed, dashed–dotted and dotted curves, respectively, represent numerical solutions obtained by applying the time-dependent boundary condition (20), the steady boundary condition (A.8) with (A.9), and the zero boundary condition at the upper boundary.

Fig. 5. Linear simulation of gravity wave propagation in the long-wave limit (ε = 0, α = 0) with different boundary conditions at the outflow (upper) boundary: the error E ( z, t ) in the numerical approximation for the wave streamfunction ψ(0, z, t ), relative to the analytical solution [16], as a function of time t at altitude z = 15. The solid curve shows the error for the results obtained with the time-dependent radiation condition (20), the dashed curve shows the error for the results obtained with the steady radiation condition (A.8) with (A.9), and the dotted curve shows the error for the results obtained with a zero upper boundary condition.

Figs. 8–11 show the results of nonlinear simulations with the amplitude parameter set to ε = 0.05 and with the other input parameters being the same as in Figs. 6 and 7. In the case shown in Figs. 8(a)–(c), the linear time-dependent radiation condition (18) is applied at the upper boundary and instabilities develop at large time; this is also seen from the dashed curve in Fig. 9. These instabilities result from a spurious growth of the zero wavenumber component of the solution (the wave-induced mean flow). Fig. 10 shows the Fourier spectrum of the streamfunction amplitude |ψˆ (κ , z, t )| as a function of wavenumber κ at fixed time and at three different levels within computational domain. The results obtained with the linear time-dependent radiation condition (18) are shown in the left column, Figs. 10(a)–(c). By this time, the zero wavenumber component of the streamfunction perturbation at z = 15 has grown to an order of magnitude larger than the component corresponding to the forced wavenumber κ = 2. This is also seen in Fig. 8(c), which shows only positive (solid) contours and no negative (dotted) contours in the upper part of the computational domain, as well as from the dashed curve in Fig. 11, which shows ˆ 0, z, t ) as a function of altitude z at this time. the wave-induced mean streamfunction ψ( The corresponding change in the mean flow velocity is u¯ 0 ( z, t ) = −ε ψˆ z (0, z, t ). It is close to zero in the upper part of the ˆ 0, z, t ) is large but approximately constant, and it is large and negative in the lower part of the domain. domain where ψ(

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

87

√

Fig. 6. Linear simulation of gravity wave propagation (ε = 0, α = 0.4): contour plots of the wave streamfunction ψ(x, z, t ) at times (a) t = 10, (b) t = 20, (c) t = 45 obtained by applying the linear time-dependent radiation condition (18) at the upper boundary.

This means that there is a spurious wave-induced mean flow vorticity ζ¯0 ( z, t ) = −u¯ 0z in the lower part of the domain. This appears within the nonlinear Jacobian term in equation (36) and acts to increase the magnitude of the higher-wavenumber components of the solution, as seen in Fig. 10(c). If the computation shown in Figs. 8(a)–(c) and 10(a)–(c) is continued further, these small-scale instabilities spread though the domain. In contrast, the addition of the nonlinear correction (35) to the radiation condition in the nonlinear simulations prevents the spurious growth of the zero wavenumber component, as well as the high-wavenumber components and thus allows the nonlinear simulations to be continued to late time without the development of significant small-scale instabilities. This is seen in Figs. 8(d)–(f) and the solid curve in Fig. 9. In Fig. 10(f) we observe that, even at late time, the high-wavenumber

88

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

√

ˆ z, t )] = ψ(0, z, t ) as a function Fig. 7. Linear simulation of gravity wave propagation (ε = 0, α = 0.4): the real part of the complex wave amplitude Re[ψ( of altitude z at time t = 45. The solid curve shows the exact steady-amplitude solution (32) and the dashed curve shows the numerical solution obtained by applying the linear time-dependent boundary condition (18) at the upper boundary.

components remain on the order of a factor of ε smaller than the component corresponding to the forced wavenumber; this is consistent with the expectations of the weakly-nonlinear analysis of Section 2. 5. Concluding remarks In this study we described the derivation and implementation of a time-dependent radiation condition for use in simulating the nonlinear evolution of forced upward-propagating internal gravity waves using a finite computational domain. First the linear time-dependent radiation condition of Campbell and Maslowe [6] was implemented in some linear simulations. This radiation condition is written in terms of a Laplace convolution integral which is nonlocal in time and expensive to compute as it requires the values of the dependent variable at all previous time levels. An approximation was used to make the computations local in time and the results were compared with those obtained using a steady radiation condition or a zero boundary condition. For the special case of the long-wave limit, the numerical results were compared with the exact analytical solution originally derived by Nadon and Campbell [16] and the results obtained with the linear time-dependent radiation condition were found to be in close agreement with the analytical solution. On the other hand, applications of steady and zero boundary conditions both gave inaccurate results. However, when the linear time-dependent radiation condition was implemented in nonlinear simulations, there was a spurious growth in the zero wavenumber component of the solution and this led to the development of instabilities. A nonlinear correction was obtained using weakly-nonlinear perturbation theory and added to the linear time-dependent boundary condition to represent the divergence of the momentum flux at the outflow boundary. This prevented the growth in the zero wavenumber component and allowed the nonlinear simulations to continue further in time without instabilities. We examined a configuration where the buoyancy frequency N and the background mean flow speed u¯ are both constant and there is no critical level so that the forced waves propagate to the outflow boundary without attenuation. The radiation condition could, however, be used in any configuration with more general representations of the background density and with vertical shear, with or without a critical level, provided the buoyancy frequency and the background mean flow speed both vary slowly enough with altitude near the outflow boundary that they can be approximated by constants in the derivation of the radiation condition. The procedure can also be implemented in other similar problems involving simulations of wave propagation in an infinite domain using a finite computational domain; for example, an implementation for Rossby waves on a horizontal plane is suggested by Victor [24]. Appendix A. Steady linear radiation condition A steady linear radiation condition for (10) is obtained by making use of the standard group velocity argument which was first presented by Booker and Bretherton [4] and is summarized here. Considering a solution of the form

ψ(x, z, t ) = φ(z)e ik(x−ct ) + c.c.

(A.1)

for (10) gives the Taylor–Goldstein equation

d2 φ dz2



−

N2

(u¯ − c )2

 + α 2 k 2 φ = 0.

(A.2)

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

89

√

Fig. 8. Nonlinear simulation of gravity wave propagation (ε = 0.05, α = 0.4): contour plots of the wave streamfunction ψ(x, z, t ). Figures (a), (b) and (c) show the numerical solution at times t = 10, 20, 45, respectively, obtained by applying the linear time-dependent radiation condition (18) without the nonlinear correction. Figures (d), (e) and (f) show the numerical solution at times t = 10, 20, 40, respectively, obtained by applying the nonlinear boundary condition (18) and (35) at the upper boundary. The solid contours represent zero and positive values of ψ and the dotted contours represent negative values.

With u¯ constant, φ takes the form e ±imz and (A.2) can be written as



d

∂z

 + im

d

∂z

− im φ( z) = 0,

(A.3)

90

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

√

Fig. 9. Nonlinear simulation of gravity wave propagation (ε = 0.05, α = 0.4): the wave streamfunction ψ(0, z, t ) as a function of time t at altitude z = 15. The dashed curve shows the results obtained by applying the linear time-dependent radiation condition (18) without the nonlinear correction and the solid curve shows the results obtained by applying the nonlinear boundary condition (18) and (35) at the upper boundary.

where

m2 =

N2

(u¯ − c )2

− α 2 k 2 > 0.

(A.4)

This gives a radiation condition

dφ dz

± imφ = 0.

(A.5)

From (A.4) and the gravity wave dispersion relation which relates the frequency

¯ ± ω(k, m) = uk

Nk

(α 2k2 + m2 )1/2

ω to the wavenumbers k and m, (A.6)

,

the vertical component of the group velocity for a solution with vertical wavenumber ±m is

c gz (±m) =

∂ω (u¯ − c )3km . =± ∂(±m) N2

(A.7)

If m is defined to have the same sign as (u¯ − c ), then c gz (m) has the same sign as k, while c gz (−m) has the opposite sign. This means that, for example, if (u¯ − c ) > 0 and k > 0, then the solution e imz has positive (upward) group velocity, while the solution e −imz has negative (downward) group velocity and we must take the negative sign in (A.5) in order to have upward-propagating waves. In general, the appropriate steady radiation condition for an upward-propagating wave is

dφ dz

− sgn(k)i |m|φ = 0,

(A.8)

where m is given by (A.4). In the special case of the long-wave limit,

α = 0 and m is simply given by

m = N /(u¯ − c ).

(A.9)

Appendix B. Horizontal-mean momentum flux divergence The integral of the nonlinear term in (7) over a horizontal wavelength 2π /k is 2π /k

2π /k 2

(ψx ψzzz + α 2 ψx ψzxx − ψz ψzzx − α 2 ψz ψxxx )dx

J (ψ, ∇ ψ)dx = 0

0 2π /k

=

2π /k

(ψx ψzzz − ψz ψzzx )dx = 0

(ψx ψzzz + ψxz ψzz )dx 0

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

91

2π /k

=

(ψx ψzzz − ψz ψzzx + ψz ψzzx + ψzx ψzz )dx 0 2π /k

∂2 (ψx ψzzz + 2ψzz ψzx + ψz ψzzx )dx = 2 ∂z

= 0

2π /k

ψx ψz dx.

(B.1)

0

√

Fig. 10. Nonlinear simulation of gravity wave propagation (ε = 0.05, α = 0.4): the Fourier spectrum of the amplitude of the wave streamfunction ˆ κ , z, t )| as a function of wavenumber κ at time t = 40 and at altitudes (a) z = 15, (b) z = 10 and (c) z = 1, obtained by applying the linear time|ψ( dependent radiation condition (18) without the nonlinear correction and the corresponding results obtained by applying the nonlinear boundary condition (18) and (35) at time t = 30 and at altitudes (d) z = 15, (e) z = 10 and (f) z = 1.

92

V. Nijimbere, L.J. Campbell / Applied Numerical Mathematics 110 (2016) 75–92

√

ˆ 0, z, t ) as a function of altitude z Fig. 11. Nonlinear simulation of gravity wave propagation (ε = 0.05, α = 0.4): the wave-induced mean streamfunction ψ( at time t = 30. The dashed curve shows the result obtained by applying the linear time-dependent radiation condition (18) without the nonlinear correction and the solid curve shows the results obtained by applying the nonlinear boundary condition (18) and (35) at the upper boundary.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, Nat. Bur. Stand., 1964, 1046 pp. P. Baines, Topographic Effects in Stratified Flows, Cambridge University Press, 1995. M. Béland, T. Warn, The radiation condition for transient Rossby waves, J. Atmos. Sci. 32 (1975) 1873–1880. J.R. Booker, F.P. Bretherton, The critical layer for gravity waves in a shear flow, J. Fluid Mech. 27 (1967) 513–539. L.J. Campbell, Wave–mean flow interactions in a forced Rossby wave packet critical layer, Stud. Appl. Math. 47 (2004) 39–85. L.J. Campbell, S.A. Maslowe, Nonlinear critical-layer evolution of a forced gravity wave packet, J. Fluid Mech. 29 (2003) 151–179. A. Eliassen, E. Palm, On the transfer of energy in stationary mountain waves, Geophys. Publ. 22 (1961) 1–23. J.E. Geisler, R.E. Dickinson, Numerical study of an interacting Rossby wave and barotropic zonal flow near a critical level, J. Atmos. Sci. 31 (1974) 946–955. D. Givoli, J.B. Keller, Nonreflecting boundary conditions for elastic waves, Wave Motion 12 (1990) 261–279. S. Goldstein, On the stability of superposed streams of fluids of different densities, Proc. R. Soc. Lond. A 132 (1931) 524–548. M.J. Grote, J.B. Keller, Exact nonreflecting boundary conditions for the time dependent wave equation, SIAM J. Appl. Math. 55 (1995) 280–297. M.J. Grote, J.B. Keller, Nonreflecting boundary conditions for Maxwell’s equations, J. Comput. Phys. 139 (1998) 327–342. H.D. Han, X.N. Wu, Approximation of infinite boundary condition and its application to finite element methods, J. Comput. Math. 3 (1985) 179–192. H.D. Han, X.N. Wu, A survey on artificial boundary methods, Sci. China Math. 56 (2013) 2439–2488. D.A. McLaren, L.J. Campbell, R. Vaillancourt, A sequential, implicit, wavelet-based solver for multi-scale time-dependent partial differential equations, Axioms 2 (2013) 142–181 (Special Issue on Wavelets and Applications). M. Nadon, L.J. Campbell, An exact expression for transient forced internal gravity waves in a Boussinesq fluid, Wave Motion 44 (2007) 340–345. C.J. Nappo, An Introduction to Atmospheric Gravity Waves, Academic Press, New York, 2002. V. Nijimbere, Ionospheric gravity wave interactions and their representation in terms of stochastic partial differential equations, Ph.D. thesis, Carleton University, 2014. V. Nijimbere, L.J. Campbell, Exact expressions for transient forced internal gravity waves and spatially-localized wave packets in a Boussinesq fluid, Wave Motion 58 (2015) 117–128. R.A. Pearson, Consistent boundary conditions for numerical models of systems that admit dispersive waves, J. Atmos. Sci. 9 (1974) 1481–1489. E.A. Spiegel, G. Veronis, On the Boussinesq approximation for compressible fluid, Astrophys. J. 6 (1960) 442–447. B.R. Sutherland, Internal Gravity Waves, Cambridge University Press, 2010, 377 pp. G.I. Taylor, Effect of variation in density on the stability of superposed streams of fluids, Proc. R. Soc. Lond. A 201 (1931) 499–523. N. Victor, Time-dependent non-reflecting boundary conditions for simulations of wave propagation in geophysical fluid flows, M.Sc. thesis, Carleton University, 2010.