Numerical effects on vertical wave propagation in deep‐atmosphere models

Ray‐tracing techniques have been used to investigate numerical effects on the propagation of acoustic and gravity waves in a non‐hydrostatic dynamical core discretized using an Arakawa C‐grid horizontal staggering of variables and a Charney–Phillips vertical staggering of variables with a semi‐implicit timestepping scheme. The space discretization places limits on resolvable wavenumbers, and redirects the group velocity and the propagation of wave energy towards the vertical. The time discretization slows the wave propagation while maintaining the group velocity direction. Wave amplitudes grow exponentially with height due to the decrease in the background density, which can cause instabilities in whole‐atmosphere models. Although molecular viscosity effectively damps the exponential growth of waves above about 150 km, additional numerical damping might be needed to prevent instabilities in the lowermost thermosphere. These results are relevant to the Met Office Unified Model, and provide insight into how the stability of the model may be improved as the model's upper boundary is raised into the thermosphere.


Introduction
It is well-known that numerical methods can affect waves in various ways. In this article, ray-tracing techniques are used to compare the vertical propagation of acoustic and gravity waves for the discrete governing equations in an atmospheric model with that for the corresponding continuous governing equations. An important application of this study is to the development of a whole-atmosphere extension of the Met Office Unified Model (UM). The UM is based on the ENDGame dynamical core (Wood et al., 2014), which solves the non-hydrostatic compressible Euler equations, and commonly runs with an upper boundary at around 80 km. Initial attempts to raise the upper boundary gave rise to instabilities (D. R. Jackson, 2015; personal communications); it is believed that these result from a combination of numerical effects and missing physical processes. The present ray-tracing study will allow some of these numerical effects to be quantified, and also for the stabilizing effect of realistic molecular viscosity and diffusion, as well the tunable numerical damping due to the semi-implicit time integration scheme, to be estimated.
The instabilities found in the UM appear to be related to the propagation of fast acoustic and gravity waves from the lower atmosphere into the thermosphere. The velocity perturbations of such waves grow roughly exponentially with height because of the roughly exponential decrease of the background density with height. The numerical methods used can affect the generation of these waves, their upward propagation, and their eventual fate when they reach large amplitudes. The main focus of the present study is on the wave propagation, but it will be important to note also the role of the other effects.
On large scales the lower atmosphere is close to hydrostatic and geostrophic balance. Unbalanced motions such as gravity waves and especially acoustic waves are generally weak. Nevertheless, a variety of mechanisms, including nonlinear spontaneous emission by the balanced flow, can generate unbalanced motions, and this wave generation can be greatly exaggerated in numerical models (e.g. Mohebalhojeh and Dritschel, 2000). We have carried out a standard baroclinic instability test case (Jablonowski and Williamson, 2006) with a semi-implicit semi-Lagrangian dynamical core that uses almost identical numerical methods to ENDGame, and found that unbalanced density and divergence perturbations, i.e. acoustic waves, are generated near the surface with relative amplitude ρ /ρ 0 of order 10 −4 to 10 −3 , where ρ is the density perturbation and ρ 0 is the background density. Although very weak, these are significantly stronger than is realistic. Moreover, because typical model grids are much finer in the vertical than in the horizontal, any marginally resolved waves will tend to be generated with nearly vertical wave vectors.
A single vertical column version of ENDGame (ENDGame1D) has been developed in order to supplement the ray-tracing calculations reported here. The column model is able to represent the purely vertical propagation of acoustic waves. Like the three-dimensional version, it uses a semi-Lagrangian advection scheme with a semi-implicit time integration scheme. When the model is initialized with a relative density perturbation 10 −3 near the bottom boundary, a wave packet propagates vertically, growing exponentially as expected, until it reaches an altitude of around 110 km, at which point the relative amplitude ρ /ρ 0 is close to order 1 and the amplitude is so large that it causes the code to fail. Clearly the semi-implicit semi-Lagrangian numerical method, which is very efficient and accurate for modelling tropospheric dynamics, is not suitable for handling very large amplitude acoustic waves that result from upward propagation.
Ray tracing has been used previously to study the propagation and dissipation of gravity waves in the real atmosphere (e.g. Vadas, 2007, and references therein). Here, the emphasis will be on how the propagation and dissipation of waves is modified in numerical models by the numerical methods used, for both gravity and acoustic waves. Because of the background temperature structure of the atmosphere, waves propagating upwards at certain angles may be refracted back towards the surface. However, it is well-known that the propagation of waves can be captured imperfectly by numerical models, for example through the appearance of computational modes and parasitic waves, distortion of group velocity, and spurious refraction, reflection, and trapping of waves (e.g. Trefethen, 1982;Vichnevetsky, 1987aVichnevetsky, , 1987bLong and Thuburn, 2011). Of particular interest is whether the numerical methods used in ENDGame, in combination with the grid anisotropy, could underestimate the downward refraction of wave packets and hence cause an excessive focussing of wave energy into the thermosphere.
A physical process that becomes very important in the thermosphere, but is currently missing from the UM, is molecular viscosity and diffusion. It might be expected that molecular viscosity and diffusion, if included, would help to control the amplitude of upward propagating waves and hence help to stabilize the model. The ray-tracing calculations here have been extended to include the effects of molecular viscosity and diffusion on wave amplitude. The results suggest that molecular visocity and diffusion are effective at controlling wave amplitudes, but only above about 150 km. This result is confirmed in experiments with ENDGame1D. Thus, waves generated in the troposphere with a typical relative density amplitude of 10 −3 become large enough to cause the model to fail around 110 km, before they reach the molecularly diffused region. In other words, damping due to realistic physical processes appears to be insufficient, on its own, to stabilize the model.
ENDGame's numerical methods contribute a certain amount of numerical damping. In particular, the semi-implicit time integration scheme includes an off-centring parameter α ∈ [0.5, 1.0] that can be varied to change the implicitness of the scheme. Setting α = 0.5 gives a centred time step which has no damping, while higher values for α give more implicit time steps, which are artificially damping to high-frequency waves but can provide more stability. The effects of varying α can be estimated and included in the ray-tracing calculation, enabling an investigation into whether a modest off-centring can stabilize upward propagating waves until they reach the altitude of strong molecular damping.
The formulation used by ENDGame, detailed by Wood et al. (2014), is considered here. However, the results should be applicable to any atmosphere models that use the same space and time discretizations. Note that results relating to acoustic waves are only applicable to non-hydrostatic models because acoustic waves are not represented by hydrostatic models (Monin and Obukhov, 1958). In the horizontal, ENDGame uses the Arakawa C-grid staggering of variables (Arakawa and Lamb, 1981) and in the vertical it uses the Charney-Phillips staggering of variables (Charney and Phillips, 1953). The other most common choice for vertical grid staggering is the Lorenz grid (Lorenz, 1960). The propagation of gravity waves would be affected by the averaging of the buoyancy term needed with the Lorenz grid, but it is not considered here. Discussions of vertical discretizations and their ability to capture wave propagation can be found in Kasahara (1974), Tokioka (1978), Leslie and Purser (1992), Fox-Rabinovitz (1994), Thuburn and Woollings (2004) and Liu (2008).
Section 2 recalls the principal ideas behind the ray-tracing technique, and derives the wave dispersion relations for the continuous and the discrete equations. Results of the ray-tracing calculations are presented in section 3 for both acoustic waves and gravity waves. Both observations (e.g. Nastrom and Gage, 1985) and model simulations (e.g. Skamarock and Klemp, 2008) indicate a broad energy spectrum in the atmosphere; therefore this study considers a wide range of wavenumbers. The results are summarized in section 4.

Ray-tracing equations
The theory of ray tracing (e.g. Lighthill, 1978) assumes that waves are locally sinusoidal in space and time, and that the scales on which the wave vector varies are large compared to the wavelength. Here, attention is restricted to propagation in a two-dimensional vertical cross-section of the atmosphere. Then a nearly monochromatic wave packet with wave vector k = (k, m) and frequency ω approximately satisfies the local dispersion relation ω = (k, m; x, z, t), where the dependence on x = (x, z) and time t can arise, for example, through variations in the background wind, temperature, or stratification. In this work, any dependence on x or t is neglected. The wave packet then propagates along a space-time trajectory (a 'ray') at the group velocity c g and the wave vector also evolves according to The derivative operator D c g /Dt = ∂/∂t + c g · ∇ x represents the rate of change of a variable with time at a position moving with the group velocity c g . It can be found from the ray-tracing equations (1) and (2) that D c g ω/Dt = 0, provided ∂ /∂t = 0, i.e. provided the background is steady. In this study, c g depends on x only through the temperature T, which is treated herein as a function of z only. Lighthill (1978) also introduces the wave energy conservation law that describes the evolution of wave amplitude: Here, W represents the wave energy per unit volume and τ represents the time-scale for any processes that act to damp the waves. Note that, for a stratified atmosphere in which the background density ρ 0 decreases exponentially with height, a constant W corresponds to an exponential growth of wave fields such as relative density perturbations ρ /ρ 0 . In this section, expressions for , ∇ k and ∇ x are obtained for the continuous and discrete cases in order to solve the ray-tracing equations (1) and (2) and the wave energy conservation law (3). The wave propagation and the rate of wave amplitude growth may then be compared when using the analytical governing equations or the spatially and temporally discrete governing equations.

Continuous equation set
The 2D Euler continuity, thermodynamic and momentum equations considered here are obtained from the continuous governing equations described by Wood et al. (2014) by neglecting the Coriolis and source terms, as follows: Dw Dt Here D/Dt is the material derivative, ρ is the density, u = (u, w) is the velocity, is the potential temperature, C p is the specific heat capacity at constant pressure, ≡ (p/p ref ) κ is the Exner pressure, where p is the pressure and p ref is a constant reference pressure, κ ≡ R/C p where R is the gas constant per unit mass, and g is the gravitational acceleration. This equation set is fully non-hydrostatic, although it can be made hydrostatic by neglecting the Dw/Dt term from the vertical momentum equation (7). The equation of state for this system of equations is

Linearized equation set
A wave equation for a single purely oscillatory variable, and hence a dispersion relation, can be obtained by linearizing the governing equations (4)-(8) about a basic state. This is done by separating the current state of the system (u, w, ρ, , ) into the basic state (u 0 , w 0 , ρ 0 , 0 , 0 ) which satisfies the governing equations at rest and in hydrostatic balance, and a small perturbation term (u , w , ρ , , ), so that any variable ϕ can be expressed as ϕ = ϕ 0 + ϕ . Note that ρ 0 = ρ 0 (z), 0 = 0 (z) and 0 = 0 (z), and products of the small perturbation terms may be neglected. The resulting linear system will have solutions proportional to exp(−iωt) for some ω, so the time derivatives can be replaced by a factor −iω, giving Here, N 2 is the buoyancy frequency, given by (Salby, 1996). The equation of state (8) is linearized as

Obtaining an equation for a single variable
Equations (9)-(12) and (14) form a system of equations for (u , w , ρ , , ). Each variable can be eliminated in turn to obtain a single equation for one unknown variable. Begin by eliminating ρ and between (9), (10) and (14), and eliminate between (10) and (12), using the fact that the basic state satisfies hydrostatic balance to simplify the result. Use the resulting two expressions to eliminate w , and substitute this result into ∂(11)/∂x to get the following equation for only: Here H ρ is the density scale height, given by Finally, divide through (15) by C p 0 and note that the speed of sound is given by (Batchelor, 2000), where T 0 ≡ 0 0 is the basic state temperature, to obtain

Obtaining a wave equation
The solutions to (18) for have a strong, nearly exponential, dependence of the amplitude on height. In order to avoid this near-exponential dependence, a change to a more wavelike variable is needed. Consideration of the case with isothermal background, for which N 2 and c 2 s are constant, indicates that the required change of variable is Making this change of variable in (18) (without assuming an isothermal background), using the identity (Daley, 1998;Thuburn et al., 2002), and defining the inverse height scale leads to In the last step a WKB approximation has been made in order to bring the factor 1/(ω 2 − N 2 ) outside a vertical deriative. Equation (22) for q agrees with the wave equation derived by Thuburn et al. (2002) when T 0 is constant.

Dispersion relations and group velocities
In this section, the wave equation (22) is used to obtain expressions for the dispersion relation and group velocity for the continuous wave equation. The changes that are made to these expressions when the wave equation is discretized in space and time are then discussed. For both the continuous and discrete cases it is assumed that q has wavelike solutions that are locally proportional to exp{i(kx + mz − ωt)}. In other words, k, m and ω are the values of wavenumber and frequency that would be obtained by applying a Fourier transform to the solution for q. The objective here is to determine the dispersion relation ω = (k, m; x, z).
If ω = (k, m; x, z) is defined to be the dispersion relation for the continuous system (22), then it turns out that the continuous and the discrete cases can all be written in the form ω = ( k, m; x, z) for particular choices of k(k), m(m), and ω(ω). Here, k, m and ω may be thought of as the effective wavenumbers and frequency as seen by the continuous dispersion relation.

Continuous dispersion relation
In the continuous case, k and ω are simply given by k = k and ω = ω, respectively. Substituting a solution for q proportional to exp{i( kx + mz)} into (22) with ω = ω gives a quadratic for ω 2 : The four solutions for ω from (23) are where If the second '±' in the solution for ω (24) is set to be a '+', then the contents of the outer square root are quite large, giving a large frequency ω, which corresponds to the two high-frequency acoustic wave solutions. Conversely, if the second '±' is set to be a '−', then the two A 2 terms almost cancel. In this case, the contents of the outer square root are quite small, giving a smaller frequency ω, which corresponds to the two lower-frequency gravity wave solutions.
As an aside, the pure acoustic wave dispersion relation can be considered, and is: where K 2 = k 2 + m 2 , which can be obtained by neglecting the buoyancy frequency N 2 and 2 in the full dispersion relation (24). Similarly, letting c s → ∞ and neglecting in the quadratic for ω 2 (23) gives the pure gravity wave dispersion relation: Next, the group velocity c g = ∇ k can be obtained as follows: In the continuous case, ∂ω/∂ ω = 1, ∂ k/∂k = 1, ∂ m/∂m = 1, so only ∇ k is needed. Define D( ω) to be the right-hand side of the quadratic for ω 2 given by (23). Then D is the derivative of D with respect to ω and is: Expressions for ∂ /∂ k and ∂ /∂ m can be found by differentiating the quadratic for ω 2 (23) with respect to k and m and rearranging: Next, ∇ x is found, which is used in the ray-tracing equation that describes the wavenumber along the trajectory given by (2). Since our background state is independent of x, ∂ /∂x = 0. Then ∇ x is as follows: so k is constant along the ray. In the analytical case, ∂ω/∂ ω = 1, so only ∂ /∂z is needed. This can be found by differentiating the quadratic for ω 2 (23): where the subscripts | k, m and | k, m, ω indicate that these values are held constant as the derivatives are taken. Evaluating (33) and rearranging gives To evaluate (34) for ∂ /∂z, analytical expressions for ∂(N 2 )/∂z, ∂(c 2 s )/∂z and ∂ /∂z are required, which are found in Appendix A.

Spatially discrete dispersion relation
In this section, expressions for k and ∂ k/∂k, ∂ m/∂m are found in the spatially discrete case, which appear in the equations for the dispersion relation (24), the group velocity (28) and ∂ /∂z (32). These equations will then describe wave propagation behaviour in the dynamical core discretized in space using an Arakawa C-grid horizontal staggering of variables (Arakawa and Lamb, 1981) and a Charney-Phillips vertical staggering of variables (Charney and Phillips, 1953).
In the spatially discrete case, derivatives are replaced with finite differences. First consider horizontal derivatives. Provided the grid resolution x is uniform (or slowly varying), it is still possible to obtain locally wavelike solutions (Long and Thuburn, 2011). However, the finite-difference approximation to the horizontal derivative of a wavelike variable is rather than ike ikx . Thus, the effective horizontal wavenumber seen by the dispersion relation in the spatially discrete case is The effect of discretization on vertical derivatives is more subtle, since the sinusoidal structure of the perturbation fields will be scaled by a roughly exponential dependence on height. Consequently, the dispersion relation then sees not only a modified effective vertical wavenumber: but also a modified effective N 2 and (Thuburn, 2006, gives details). However, provided the vertical resolution z is finer than the density scale height (16), as is usually the case in numerical models, the effective N 2 and should be close to their true values. Therefore, to simplify the calculations below, N 2 and are taken to have their true values. These expressions for k and m are then used in (24) for ω and in (30), (31) and (34) for ∂ /∂ k, ∂ /∂ m and ∂ /∂z respectively. The group velocity is given by (28), but ∂ k/∂k and ∂ m/∂m are now found by differentiating (36) and (37): Note in this study that x and z are assumed to be uniform in space. Non-uniformity of x or z would introduce an explicit spatial dependence into the expressions for k and m ( (36) and (37)), which would then require additional terms in (32) and (33).

Temporally discrete dispersion relation
In this section, expressions for ω and ∂ω/∂ ω are found in the temporally discrete case, which appears in the equations for the dispersion relation (24), the group velocity (28) and ∂ /∂z (32). These equations will then describe wave propagation behaviour in the dynamical core discretized in time using a semi-implicit timestepping scheme, as described by Wood et al. (2014).
Recall that the linearized equation set (9)- (14) is of the form δ t ϕ + L(ϕ) = 0, with i ω as an eigenvalue for the L operator. Discretizing δ t ϕ between times n and n + 1 for one eigenmode using a centred-in-time semi-implicit scheme gives: For wavelike solutions, it is expected that The use of this expression and the identity gives ω in terms of ω: In the temporally discrete case, the group velocity is given by (28) and ∇ x is given by (32), but ∂ω/∂ ω is now found by differentiating (40) for ω as follows:

Time-scales for decaying processes
In this section, expressions for the time-scales of decaying processes are found, so that they may be used with the wave energy conservation law (3) to describe how they affect the growth of wave amplitudes. The effect of off-centring parameters in the semi-implicit scheme can be considered by finding the resulting decaying time-scale. Discretizing δ t ϕ, for some variable ϕ, between times n and n + 1 using a semi-implicit scheme with off-centring parameter α ∈ [0.5, 1.0] gives the following: Rearranging and expressing ϕ (n+1) in terms of ϕ (n) gives: Rearranging this equation gives an expression for ω: Note also that ω is comprised of a real part ω that describes oscillations, and an imaginary part −1/τ α that describes the decaying time-scale, as follows:  (Hedin 1991), and the corresponding decaying time-scale for molecular viscosity and diffusion τ diff at various altitudes found using (50), with a total wavenumber of K = 10 −3 rad m −1 . The imaginary part of (44) gives the required time-scale for offcentring which can be used with the wave energy conservation law (3) to describe the effect of off-centring on the growth of wave amplitudes: The effect of molecular viscosity and diffusion can also be considered by finding its decaying time-scale. Consider the temperature diffusion equation: where λ = C p μ/P r is the thermal conductivity, with P r ≈ 0.7 the Prandtl number and μ the scalar coefficient of molecular viscosity, which may be approximated as (Vadas 2007). Considering the approximate magnitude of the individual terms in (47) gives the following: Here, τ diff is the decaying time-scale for molecular diffusion and T is a temperature perturbation. The time-scale for molecular diffusion is then just For simplicity, this is used by the wave energy conservation law (3) to describe the effect of molecular diffusion on the growth of wave amplitudes. Values for τ diff at different altitudes are shown in Table 1: a smaller time-scale means that molecular viscosity and diffusion act more quickly. Of interest here is whether these damping mechanisms are strong enough to offset the non-dissipative wave growth due to the background density variation. In the absence of dissipation, the velocity perturbations u and relative density perturbations ρ /ρ 0 grow in proportion to ρ ∝ exp(z/2H ρ ). Hence, the growth timescale following an acoustic wave packet is given by: Whichever of the decaying time-scales is most significant will effectively damp wave growth in the dynamical core if the following is satisfied: Spatially-discrete wave on a uniform grid Altitude (km) Altitude (km) Horizontal distance (km)

Ray-tracing simulations
In this section, acoustic and gravity waves are simulated using the 2D wave equation with a variation in background density (22) by solving the ray-tracing equations (1)-(3). A sufficiently small time step is used such that numerical errors are negligible in the solutions of (1)-(3). Note that this time step is unrelated to the time step t of the semi-implicit scheme introduced in section 2.5.3, and discussed further in sections 3.2 and 3.5. The purpose of these calculations is to explore how wave propagation in the 'real' atmosphere, simulated using the analytical wave equation, differs from wave propagation with a discrete Arakawa C-grid horizontal staggering of variables (Arakawa and Lamb, 1981), a Charney-Phillips vertical staggering of variables (Charney and Phillips, 1953) and a semi-implicit timestepping scheme, using the spatially and temporally discrete wave equations.

Acoustic waves in a spatially discrete grid
Here, acoustic waves are simulated by using the solution for the dispersion relation (24) that corresponds to the frequency for acoustic waves. Figure 1 shows a comparison of ray plots: the trajectories of wave packets, between the analytical and spatially discrete wave equations for rays initiated with a small total wavenumber K = 4 × 10 −5 rad m −1 and a range of wavevector angles from 0 to π/2. The background US Standard Atmosphere (USSA) temperature profile (COESA, 1976) has the effect of refracting waves and, for small wavenumbers up to 4 × 10 −5 rad m −1 , waves are refracted back to the surface before passing the tropopause, as in Figure 1. An examination of rays with varying wavenumbers shows that rays initiated with total wavenumbers up to 7 × 10 −4 rad m −1 propagate to an altitude of no higher than 110 km before turning over.
However, the coarse horizontal resolution in the model's Arakawa C-grid ( x ∼100 km in operational climate simulations) causes waves using the spatially discrete wave equation to be less well-resolved in the horizontal direction than in the vertical direction which has a much finer grid spacing ( z ∼1 km). In the analytical case, the group velocity is approximately parallel to the wavevector; but this is no longer necessarily true in the spatially discrete case, as demonstrated later in Figure 4. This means that even waves initiated with low wavevector angles in the spatially discrete case have high group velocity angles and can reach up to 110 km before turning over.
Total wavenumber K (rad m −1 ) Critical angle for wave turning (rad) Lowest angle for resolvable k Region of WKB validity Acousti canalytical disp. relation Full analytical disp. relation Discrete-in-space Figure 2. The critical wavevector angle (from parallel to the surface, 0, to perpendicular to the surface, π/2) up to which acoustic wave turning occurs, for a range of total wavenumbers K. The dotted black line represents the critical angles for the analytical acoustic wave dispersion relation (26), and the solid black line represents the critical angles for the acoustic wave solution of the full analytical dispersion relation (24). The dashed-dotted black line represents the critical angle for the spatially discrete wave equation with x = 10 km and z = 1 km. In order for the horizontal wavenumber k to be resolved horizontally, it is required that k < π/ x. This restriction in the numerical case is illustrated with the grey area showing where k is not resolved for this choice of x. The area to the right of the grey dashed line indicates wavenumbers that satisfy WKB theory, although it is shown in Appendix B that WKB theory can still be used to describe the behaviour of wave propagation outside this limit.
However, wave turning does not always occur; whether or not waves turn over also depends on the total wavenumber K and the initial angle of the wavevector k = (k, m). This relationship can be seen in Figure 2. This plot contains a lot of information and will be referred to a number of times in this section.
In the most general case, waves turn over at the height where the vertical component of the group velocity c g is no longer positive, i.e. ∂ /∂ m = 0. From the general expression for ∂ /∂ m (31), it follows that acoustic waves turn over at the height where m = 0. For pure acoustic waves described by the pure acoustic wave dispersion relation (26), the height at which waves turn over, where m = 0, corresponds to the height at which c 2 s (z) = ω 2 / k 2 . The dotted black line in Figure 2 represents the critical wavevector angle of acoustic waves for a range of total wavenumbers K, using the pure acoustic wave dispersion relation (26). The critical wavevector angle is the angle measured from the horizontal up to which rays with a particular wavenumber will be refracted and turn over, and above which these rays will propagate upwards into the thermosphere. In reality, these waves would be very strongly dissipated above ∼150 km by molecular viscosity and diffusion (as discussed further in section 3.3), and this would need to be reflected in the dynamical core either with some form of numerical damping or in a more physically sound way with molecular viscosity and diffusion. The critical angle for the pure acoustic wave dispersion relation is a constant ∼0.32π for all wavenumbers, because it can be seen from (26) that its group velocity is independent of K. Note that the wavevector angles and the group velocity angles do not necessarily match: their relationship for acoustic waves is shown later in Figure 4.
For acoustic waves described by the full dispersion relation (24), the height at which waves turn over, where m = 0, is the height at which the following expression is satisfied: where c s , N and are all functions of height z. The solid black line in Figure 2 represents the critical angle of acoustic waves for a range of total wavenumbers K, using the acoustic wave solution to the full dispersion relation (24). Waves with small wavenumbers up to 7 × 10 −5 rad m −1 always turn over, even if they are initiated with vertical wavevectors. This is due to the relatively large buoyancy frequency N 2 and that act to bring the expression in (53) to 0 at a lower height z, causing the wave to turn over sooner. This is illustrated by the ray plots of Figure 1, where the waves always turn over and no wave energy reaches the thermosphere. For large wavenumbers and higher frequency waves, the buoyancy frequency N 2 and 2 become negligible, so the critical angle is the same as for the pure acoustic wave dispersion relation. The unusual feature of the critical angle line for acoustic waves described by the full dispersion relation (24) is the dip that occurs for total wavenumbers just above 7 × 10 −5 rad m −1 . This warrants a closer analysis: a sample of rays of varying wavenumbers near this dip at an angle of 0.3π are shown in Figure 3. For total wavenumbers just to the left of the dip: between 6 × 10 −5 and 7 × 10 −5 rad m −1 , the waves turn over very sharply at about 110 km. It turns out that there is a feature of the USSA temperature profile at 110 km that makes N 2 particularly large at this altitude (this is shown and described in Appendix A). This acts to bring the expression in (53) to 0, and hence, m to 0 at about 110 km, causing the wave to turn over.
For wavenumbers slightly larger than 7 × 10 −5 rad m −1 with an initial wavevector angle of 0.3π , the peak of N 2 is no longer sufficient to bring the expression in (53) to 0 at 110 km. Hence, these waves do not turn over and instead propagate to the top of the thermosphere, as can be seen by the dashed-dotted line in Figure 3, and is indicated in Figure 2 by the dip in the solid black line. As the wavenumber increases further, the waves become more like pure acoustic waves as the large wavenumbers k and m dominate the gravity terms in the full dispersion relation (24). Therefore, the waves are more easily refracted by the higher temperatures in the thermosphere, and turn over more quickly for increasingly high wavenumbers, as illustrated by the grey line in Figure 3.
The dashed-dotted line on Figure 2 shows the critical angle for the spatially discrete wave equation. However, the region marked by the grey area denotes wavevectors for acoustic waves that are unable to be resolved in a discrete grid with x = 10 km. For the horizontal component of waves to be resolved, the horizontal wavenumber must satisfy the condition that 0 < k < π/ x.
Where the grey area of Figure 2 puts a lower limit on the wavevector angles that can be used with the spatially discrete wave equation, an interesting phenomenon occurs. Figure 4 compares the initial angles of the wavevector and group velocity for analytical and spatially discrete acoustic waves. In the analytical case, the initial group velocity is almost parallel to the initial wavevector. However, as the discretized case is only resolved for 0 < k < π/ x, only a small range of wavevector angles, Initial angle of wave vector (rad) Initial angle of group velocity (rad) K = 10 −4 , Analytical K = 10 −4 , Discrete-in-space Where k = π/ x, the angle of the wavevector is as low as possible in order for the wave to still be resolved. The analytical group velocity c g(an) goes in the same direction as the wavevector, but the numerical group velocity c g(num) is still vertical.
between ∼0.4π and 0.5π , are able to be considered in this particular case, so that the horizontal wavenumber k is not too big. For the higher allowable values of k, and hence lower initial wavevector angles, the initial group velocity angle diverges and becomes much higher than the initial wavevector angle. This phenomenon can be explained by the group velocity equation (28). In the spatially discrete case, the ∂ k/∂k factor of the horizontal component of the group velocity is given by (38). This term cos(k x/2), approaches 0 as k approaches its upper limit π/ x. This means that, for waves initiated with the lowest resolvable wavevector angle, the horizontal component of the group velocity is 0: so only the vertical component of the group velocity has a positive value, and hence the wave propagates vertically. This may cause excessive amounts of wave energy to be channelled upwards into the thermosphere in the dynamical core, when in the real atmosphere more wave energy would be channelled horizontally. This effect is illustrated in Figure 5.
Referring back to the critical angle plot (Figure 2), the thick dashed grey line indicates that WKB theory is valid to the right of the line. For WKB theory to be valid, it is ideal for the wavelength to be significantly shorter than the length-scale L at which the background varies (L ∼10 km for z) (Howison, 2005), or equivalently: where θ is the angle of the wavevector. Note that while this line provides a rough indication of the conditions of the validity of WKB theory, it should still be indicative of the behaviour of wave propagation outside the WKB limit. Justification for this statement is given in Appendix B, and is supported by Dingle (1973) and Vadas (2007). In the critical angle plot (Figure 2), it can be seen that there is a small range of wavevectors where analytical waves will turn over but spatially discrete waves will propagate to the top of the thermosphere: some of the rays in this region are shown in Figure 6. The rays in the spatially discrete case are not able to propagate as far in the horizontal direction because the waves are not well resolved horizontally. Another feature is that the initial group velocity is at a lower angle for higher wavenumbers, but this changes at ∼280 km where the rays appear to cross over, and the group velocity for the higher wavenumbers becomes more vertical. This again is due to the cos(k x/2) factor which is applied to the horizontal component of the group velocity, influencing the angle of wave propagation.
Next, the analytical and spatially discrete propagation of waves for higher total wavenumbers are compared. Figure 7 shows another comparison of analytical and spatially discrete wave propagation, as in Figure 1, but for a larger total wavenumber of 1 × 10 −4 rad m −1 for a small range of wavevector angles from ∼0.4π to 0.5π , for which the spatially discrete waves are resolved in the horizontal direction. For this wavenumber, there are no circumstances in which a wave in the spatially discrete case can be refracted and turn over. In the spatially discrete case, waves are also unable to propagate as far in the horizontal direction as in

Altitude (km)
Horizontal distance (km) Analytical wave in USSA temperature profile Spatially-discrete wave in USSA temperature profile Figure 8. Ray plot of acoustic waves with a total wavenumber of 1 × 10 −3 rad m −1 initiated with a wavevector angle of 0.3π using either the analytical wave equation or the temporally discrete wave equation with a centred-in-time (α = 0.5) time step of 3 s, and using the USSA temperature profile. The dots represent 1 min intervals, and the simulation was run for 30 min. the analytical case. This is due to the effect of the cos(k x/2) term on the horizontal component of the group velocity preventing many waves from being resolved in the horizontal direction. Figure 4 shows that each initial group velocity angle corresponds to two potential initial wavevectors, so each red curve in Figure 7 corresponds to two wavevector angles.

Acoustic waves with a semi-implicit scheme
Here, acoustic wave propagation is simulated with the semiimplicit time discretization. Figure 8 demonstrates the slowing effect that the temporally discrete wave equation has on acoustic wave propagation. Here, a simulation has been run for 30 min, and the time step for the temporally discrete case is just 3 s: an unrealistically short time step chosen for illustrative purposes. The rays both have the same trajectory, but it can be seen that the wave in the temporally discrete case travels at a slower speed. This phenomenon can be explained by the equation for group velocity (28). In the temporally discrete case, the ∂ω/∂ ω factor of the group velocity is given by (41), which is constant along the ray, as D c g ω/Dt = 0. This acts to slow the horizontal and vertical components of the group velocity by the same constant factor while maintaining the angle of propagation.
The larger the value of ω t, the greater the slowing effect the time discretization has on the wave. For a larger more realistic time step, an acoustic wave propagates upwards extremely slowly,

Altitude (km)
Wave amplitude growth factor Analytical in isothermal atmosphere Analytical/temporally-discrete in USSA Temporally-discrete with molecular diffusion Figure 9. The wave amplitude growth factor at each height for a vertically propagating acoustic wave with a total wavenumber of 1 × 10 −4 rad m −1 under different conditions: 1. using the analytical wave equation with an isothermal atmosphere of 250 K, 2. using the analytical wave equation with the USSA temperature profile, 3. using the temporally discrete wave equation with centred-in-time (α = 0.5) 5 min time steps and the USSA temperature profile (curves 2 and 3 are coincident; see section 3.3) and 4. using the temporally discrete wave equation with centred-in-time (α = 0.5) 5 min time steps and molecular viscosity with the USSA temperature profile. but for waves with smaller frequencies ω, the time discretization has less of a slowing effect. For this reason, gravity waves are less slowed by the time discretization, as shown in section 3.5.
The experiments performed in this section are not discretein-space. However, when including the effects of both space and time discretizations, they interact as expected.

Wave amplitude growth of acoustic waves
Here, the evolution of the wave amplitude growth factor W for vertically propagating acoustic waves is considered. W represents the wave energy per unit volume, or the wave amplitude growth factor, if W = 1 is taken at the launch of the acoustic wave (Lighthill, 1978). In a stratified atmosphere, a constant value for W corresponds to an exponential growth of the relative density ρ /ρ 0 , while a strong drop-off in W is required for the ρ /ρ to remain small. Note from the wave energy conservation law (3) that if dissipation is negligible and the background is independent of x, then c z g W is constant along the ray (where c z g denotes the vertical component of c g ), and so the vertical profile of W reflects the vertical profile of c g (k, T), which in turn reflects that of the background temperature profile. This is why the plots of W resemble the mirror image of the USSA temperature profile. Figure 9 shows how W changes with height for vertically propagating acoustic waves with a total wavenumber of K = 1 × 10 −4 rad m −1 with the analytical wave equation and the temporally discrete wave equation with centred-in-time (α = 0.5) time steps, with and without the decaying time-scale for molecular viscosity from (50). The values for W from the analytical and temporally discrete wave equations are both given by the grey line in Figure 9. This is because ω and ω are both constant along the ray, so the wave amplitude growth factors in each of these cases follow the same line, but the waves propagate upwards at different speeds.
The wave amplitude growth factor from the temporally discrete wave equation with the decaying timescale for molecular viscosity τ diff (50) is given by the black line in Figure 9. Below 130km, W is indistinguishable between the cases with and without molecular viscosity. Between 130km and 150 km, the decaying effect of molecular viscosity makes W noticeably smaller, and above 200 km, molecular viscosity makes W decrease very rapidly, which can be seen more clearly in the log plot for W later in Figure  18, and effectively prevents the growth of waves above this altitude. Including molecular viscosity in the dynamical core would not only improve its accuracy by incorporating a new physical process Altitude (km) Wave amplitude growth factor that has a significant effect in the thermosphere, but it would likely also improve its stability by preventing excessive wave growth in the upper atmosphere of the model.
The features that can be observed in Figure 9 (such as the sharp increase in W due to the peak in N 2 at 110 km, and the altitude at which molecular viscosity begins to have a significant damping effect) are dependent on the total wavenumber K. Figure 10 shows how sensitive W is to different values of K. The reasons for this dependence are twofold: first, diffusive damping varies like the diffusion timescale τ diff (50), which is proportional to K 2 . Second, the time discretisation also has the effect of slowing the group velocity of the acoustic wave by the factor ∂ω/∂ ω (41). Using the approximate acoustic wave dispersion relation (26), the propagation time-scale τ prop can be given by τ prop ∝ ∂ ω/∂ω ≈ (c s K t/2) 2 + 1.
For smaller wavenumbers/bigger wavelengths, W gets very large and diffusion does not have a significant damping effect until the wave reaches an altitude above 150 km. Fortunately, for wavenumbers smaller than 8 × 10 −5 rad m −1 , the waves turn over before the wave amplitudes get too large. For larger wavenumbers/smaller wavelengths, W does not grow much, and diffusion has a more significant damping effect at lower altitudes.
The decaying time-scale for the off-centring parameter from the semi-implicit timestepping scheme (46) can be altered to account for different amounts of off-centring, as shown in Figure 11. It can be seen that increasing levels of off-centring artificially damps the wave amplitude growth throughout the whole atmosphere, regardless of whether the time-scale for molecular viscosity is also included. These results also vary depending on the total wavenumber K that is used. It must be considered that the damping also has to compensate for the exponential growth of the wave, which means that the rapid decrease in W with the increase in off-centring shown in Figure 11 is not as severe as it might appear.
In Figure 11, a small time step of t = 1 minute is used for illustrative purposes. For larger t, the damping effect of the higher off-centring parameters is even greater. This is due to the contribution of the effects of the decaying time-scale for off-centring (46) which depends on ω t and α, and of the propagation timescale τ prop which depends on t. The effect of changing t on the time taken for the wave to propagate up to 600 km is shown in Table 2.

Gravity waves in a spatially discrete grid
Here, gravity waves are simulated by using the solution for the dispersion relation (24) that corresponds to the frequency for gravity waves. Gravity waves are characterised by their smaller frequencies ω and slower group velocities c g that propagate Altitude (km) Wave amplitude growth factor α = 0.5 α = 0.501 α = 0.505 α = 0.51 α = 0.52 α = 0.53 Figure 11. The wave amplitude growth factor at each height for a vertically propagating acoustic wave with a total wavenumber of 1 × 10 −4 rad m −1 using the discretized-in-time wave equation with 1 min time steps, with no molecular viscosity but varying amounts of off-centring, and the USSA temperature profile.
[Colour figure can be viewed at wileyonlinelibrary.com]. approximately perpendicularly to the wavevector. Note that gravity waves initiated with a wavevector angle of π/2 will not propagate: from the pure gravity wave dispersion relation (27), it can be seen that at this angle, there is no gravity wave mechanism to produce gravity waves. Gravity waves initiated with small total wavenumbers K = 2 × 10 −5 rad m −1 (chosen to look at cases where the rays always turn over) at a range of wavevector angles from π/2 to π using the analytical and spatially discrete wave equations are shown in Figure 12. It can be seen in Figure 12 that the propagation of gravity waves with a small total wavenumber is not affected much by the space discretisation, because the gravity wave is still well resolved in the grid for this wavenumber.
As with acoustic waves, whether or not gravity waves turn over due to the background USSA temperature profile depends on the total wavenumber K and the initial angle of the wavevector. This relationship can be seen in Figure 13. For pure gravity waves described by the pure gravity wave dispersion relation (27), the height at which waves turn over, where m = 0, corresponds to the height at which N → ω + . In addition, the group velocity obtained from the pure gravity wave dispersion relation is given by The dotted black line in Figure 13 represents the critical angle of gravity waves for a range of K using the pure gravity wave dispersion relation (27). The critical wavevector angle here is the angle measured from π/2 to π , up to which the wave reaches an altitude of 600 km, and above which the wave is refracted and turns over. The critical angle for the pure gravity wave dispersion relation is a constant ∼0.78π for all wavenumbers, because its group velocity is independent of K. In the analytical gravity wave case, group velocity is almost perpendicular to the wavevector, but this is no longer true in the spatially discrete case, as demonstrated later in Figure 15. Horizontal distance (km) Figure 12. A comparison of ray plots of gravity waves with a total wavenumber of 2 × 10 −5 rad m −1 initiated with a range of wavevector angles measured from the horizontal, from π/2 to π, using the USSA temperature profile. Total wavenumber K (rad m −1 ) Critical angle for wave turning (rad) Discrete-in-space k Full analytical disp. relation Gravity wave analytical disp. relation Lowest angle for resolvable k Region of WKB stability Figure 13. The critical wavevector angle (from perpendicular to the surface, π/2, to parallel to the surface, π) above which wave turning occurs, for a range of total wavenumbers K. The dotted black line represents the critical angles for the analytical gravity wave dispersion relation (27), and the solid black line represents the critical angles for the gravity wave solution of the full analytical dispersion relation (24). The dashed-dotted line represents the critical angle for the spatially discrete wave equation with x = 10 km and z = 1 km. In order for the horizontal wavenumber k to be resolved horizontally, it is required that k < π/ x. This restriction in the numerical case is illustrated by the grey area showing where k is not resolved for this choice of x. The area to the right of the grey dashed line indicates wavenumbers that satisfy WKB theory, although it is shown in Appendix B that WKB theory can still be used to describe the behaviour of wave propagation outside this limit.
For gravity waves described by the full dispersion relation (24), the height at which waves turn over, where m = 0, is the height at which the expression in (53) is satisfied. The solid black line in Figure 13 represents the critical angle of gravity waves for a range of total wavenumbers K, using the gravity wave solution to the full dispersion relation (24). For large wavenumbers and higher frequency waves, the buoyancy frequency N 2 and become negligible, so the critical angle is the same as for the pure gravity wave dispersion relation.
The unusual feature of this critical angle line is the peak that occurs between the total wavenumbers 2 × 10 −5 and 1 × 10 −4 rad m −1 . This also warrants a closer analysis: a sample of rays of varying wavenumbers near the peak at a wavevector angle of 0.9π are shown in Figure 14. For total wavenumbers Horizontal distance (km) Altitude (km) K = 4 × 10 −5 K = 5 × 10 −5 K = 7 × 10 −5 Figure 14. Ray plot of gravity waves using the analytical wave equation with an initial wavevector angle of 0.9π from the horizontal, for a range of total wavenumbers near the peak on the critical angle plot in Figure 13.
to the left of the peak in Figure 13, gravity waves turn over very quickly, never propagating above the tropopause at 10 km before turning over, as in Figure 12. On the right side of the peak, if the wave does turn over, it still reaches a high altitude close to 600 km. These waves slow down greatly before turning over, giving the appearance on the ray plot that they turn over very sharply.
The slowing of gravity waves as they reach the upper atmosphere is due to N 2 becoming small at high altitudes due to the USSA temperature profile. The small N 2 feeds into the group velocity, an estimate for which is given by (55). From this, it can be seen that as m → 0, both c x g → 0 and c z g → 0, but c x g → 0 at a faster rate, so the horizontal group velocity becomes smaller more quickly, resulting in the cusp of the rays that turn over as seen in Figures 14, 16 and 17.
As the wavenumber increases further, the large k and m terms dominate the gravity terms in the full dispersion relation (24). Therefore, the gravity waves are more easily refracted by the higher temperatures in the thermosphere and turn over more quickly for increasingly high wavenumbers.
The dashed-dotted line on Figure 13 shows the critical angle for the spatially discrete wave equation. The region marked by the grey area denotes wavevectors for gravity waves that are unable to be resolved in a discrete grid with x = 10 km. Where this puts an upper limit on the wavevector angles that can be used with the spatially discrete wave equation, since the condition that 0 < k < π/ x must be satisfied, the same phenomenon occurs as for spatially discrete acoustic waves in Figure 4. Figure 15 compares the initial angles of the wavevector and group velocity for analytical and spatially discrete gravity waves. The pure gravity wave dispersion relation (27) suggests that the initial group velocity is perpendicular to the initial wavevector angle. When using the analytical gravity wave solution to the full dispersion relation (24), the initial group velocity is similar to the pure gravity wave case for bigger total wavenumbers K. As the spatially discrete case is only resolved for 0 < k < π/ x, only a small range of wavevector angles between 0.5π and ∼0.6π can be considered in this particular case for x = 100 km. For the higher allowable values of k, and hence higher initial wavevector angles, the initial group velocity angle diverges from the angle perpendicular to the initial wavevector angle.
This phenomenon is equivalent to what happens to spatially discrete acoustic waves: the cos(k x/2) term in the horizontal component of the group velocity (28) approaches 0 as k reaches its upper limit π/ x, meaning that only the vertical component of group velocity has a positive value, which forces the wave straight upwards.
In the critical angle plot: Figure 13, it can be seen that there is a small range of wavevectors where analytical waves will turn over, but spatially discrete waves will propagate to the top of the thermosphere: some of the rays in this region are shown in Figure 16. For higher wavenumbers (Figure 16), which are less well resolved in the grid, the space discretisation significantly affects Initial angle of wave vector (rad) Initial angle of group velocity (rad) K = 10 −4 , Analytical gravity wave K = 10 −4 , Spatially-discrete gravity wave K = 10 −4 , Pure gravity wave dispersion relation Figure 15. Initial group velocity angles for a range of initial wavevector angles for gravity waves with a total wavenumber of 1 × 10 −4 rad m −1 , either using the pure gravity wave dispersion relation, the gravity wave solution to the full dispersion relation or the spatially discrete wave equation with x = 100 km, and using the USSA temperature profile. Horizontal distance (km) Figure 16. A comparison of ray plots of gravity waves with a total wavenumber of 2 × 10 −4 rad m −1 initiated with a range of wavevector angles measured from the horizontal from 0.77π to 0.83π, using the USSA temperature profile. These plots illustrate the difference in ray propagation where the lines for the critical angle of the analytical and spatially discrete wave equations diverge in Figure 13. the extent of horizontal and vertical gravity wave propagation compared to well resolved small wavenumbers ( Figure 12).

Gravity waves with a semi-implicit scheme
Here, gravity wave propagation is simulated with the semiimplicit time discretisation. Figure 17 demonstrates the slowing effect that the temporally discrete wave equation has on gravity wave propagation. Here, a simulation has been run for 8 h, and the time step for the temporally discrete case is 5 min. The rays both have the same trajectory, but it can be seen that the wave in the temporally discrete case travels at a much slower speed. This is explained in the same way as for temporally discrete acoustic waves: the ∂ω/∂ ω factor of the group velocity (28) takes a value less than 1 in the temporally discrete case, which acts to slow the horizontal and vertical components of the group velocity by the same constant factor while maintaining the angle of propagation. The wave can be seen to slow down as it reaches the cusp, before speeding up again after it turns over.
The time discretization has a lesser slowing effect on gravity waves than on acoustic waves: a 3 s time step, as used in Figure 8 for the temporally discrete acoustic wave, has a negligible effect on gravity waves, whereas it severely slowed the propagation of an acoustic wave. A 5 min time step would have an even greater Altitude (km) Horizontal distance (km) Analytical wave in USSA temp. profile Spatially-discrete wave in USSA temp. profile Figure 17. Ray plot of gravity waves with a total wavenumber of 1 × 10 −3 rad m −1 initiated with a wavevector angle of 0.85π using either the analytical wave equation or the temporally discrete wave equation with a centred-in-time (α = 0.5) time step of 5 min, and using the USSA temperature profile. The dots represent 1 h intervals, and the simulation was run for 8 h.

Altitude (km)
Wave amplitude growth factor Acoustic wave with molecular diffusion Gravity wave with molecular diffusion Figure 18. Log plot of the wave amplitude growth factor of acoustic waves and gravity waves propagating with an initial wavevector angle of 0.7π and a total wavenumber of 1 × 10 −4 rad m −1 using the temporally discrete wave equation with centred-in-time (α = 0.5) 5 min time steps and molecular viscosity with the USSA temperature profile. slowing effect on acoustic waves. The amount of slowing caused by the time discretisation is again governed by the size of ω t: as ω is much smaller for gravity waves, ω t is also much smaller, so gravity waves are less slowed by the time discretization.

Wave amplitude growth of gravity waves
Here, the evolution of the wave amplitude growth factor W for gravity waves is considered. In general, it behaves in a similar way to acoustic waves, but there are a few differences: in Figure 18, it can be seen that that the molecular viscosity begins to have a significant effect on gravity waves at a lower altitude than it does for acoustic waves. This is because gravity waves propagate more slowly, so there is more time for molecular viscosity to act on them, and so it can act to damp gravity waves more effectively at slightly lower altitudes. The altitude at which molecular viscosity begins to have a significant damping effect on gravity waves is also dependent on the total wavenumber K, similarly to acoustic waves in Figure 10.
The effect of the decaying time-scale for the off-centring parameter (46) from the semi-implicit timestepping scheme on the wave amplitude growth factor is shown in Figure 19 for different amounts of off-centring. As for acoustic waves, increasing the amount of off-centring artificially damps the wave amplitude growth throughout the whole atmosphere. It can be seen by comparing with Figure 11 that off-centring parameters have less of a damping effect on low-frequency gravity waves than on high-frequency acoustic waves.

Summary
In this study, ray-tracing experiments were performed to study the effects of numerical methods on wave propagation in numerical models.

Altitude (km)
Wave amplitude growth factor α = 0.5 α = 0.501 α = 0.505 α = 0.51 α = 0.52 α = 0.53 α = 0.54 α = 0.55 Figure 19. The wave amplitude growth factor at each height for a gravity wave propagating at an angle of 0.7π with a total wavenumber of 1 × 10 −4 rad m −1 using the temporally discrete wave equation with 1 min time steps, with varying amounts of off-centring, and the USSA temperature profile. [Colour figure can be viewed at wileyonlinelibrary.com].
A coarse horizontal resolution in numerical models prevents waves from being well-resolved in the horizontal direction. The limit on resolvable horizontal wavenumbers has a significant effect on wave propagation, redirecting waves vertically that would usually propagate at much lower angles, causing excessive, unphysical amounts of wave energy to be transferred into the upper atmosphere.
An important effect that has not been considered in this ray-tracing analysis is that of varying grid model spacing. A coarser vertical resolution at higher altitudes would trap short wavelengths at low altitudes, preventing them from propagating upwards and transferring excessive amounts of wave energy into the upper atmosphere. It would therefore likely have a positive effect on dynamical core stability, but its effect on wave propagation has not been studied here.
The semi-implicit timestepping scheme is severely slowing to high-frequency acoustic waves, but less so to gravity waves. On its own, it has no effect on the direction of wave propagation.
The inclusion of the timescale for molecular viscosity and diffusion confirms that molecular viscosity has a significant damping effect on wave amplitudes in the thermosphere. This suggests that the inclusion of molecular viscosity and diffusion in numerical models should be beneficial for the model's stability when simulating the thermosphere by preventing excessive wave growth at high altitudes. Simulating acoustic waves in non-hydrostatic models seems potentially more problematic than simulating gravity waves because they are less damped by molecular viscosity. On the other hand, they are more strongly damped by the off-centring of the semi-implicit timestep.
The improvement of the model's stability with the inclusion of molecular viscosity and diffusion has been confirmed by experiments with ENDGame1D. A small amount of off-centring is required in order to extend the top boundary up from 110 km to the molecularly diffused region above 150 km, whereupon the implementation of molecular viscosity acts to regulate wave amplitude growth, keeping the model stable up to altitudes of 600 km. These results also provide guidance for how to proceed with improving the stability of the UM and the full 3D version of the ENDGame dynamical core as they are extended into whole-atmosphere models.

B. The WKB approximation
In section 3, a number of experiments are performed that are outside the limit of validity of WKB theory. However, the assertion is made that WKB theory can still be used to describe the behaviour of wave propagation outside this limit. In this appendix, this assertion is justified. Dingle (1973) discusses the WKB approximation, and states that the accuracy of solutions obtained through this approximation are high when either (i) the wavevector k is slowly varying, or (ii) k is large in magnitude. In section 2.4, it is assumed that k varies slowly with height z, but k is not large enough in magnitude for the wavelength to be significantly shorter than the background length-scale. To confirm whether changes in the wavenumber and wave amplitude agree with the WKB theory, experiments have been performed to compare the wavenumbers and wave amplitudes produced by the ray-tracing scheme with those obtained from ENDGame1D.
Vertically propagating waves are simulated in ENDGame1D by forcing an oscillation at the bottom boundary. Once it reaches a steady state, the wave amplitudes at each height level i of ENDGame1D from (19), accounting for temperature variations, can be recorded as where 0 indicates the linearly interpolated value and is the perturbation from the initial . The local wave amplitude q 0i , wavenumber m i and phase φ i can be estimated by fitting a function q i = q 0i cos(m i z i + φ i ) ( B 2 ) to the data over a set of (2n + 1) grid points near level i (n = 2 was found to be sufficient). The fits for q 0i , m i and φ i are found by minimizing the mean square error (MSE i ) given by using a minimization scheme such as the downhill simplex method (Nelder and Mead, 1965). The wavenumbers m i and wave amplitudes q 0i can then be compared with outputs from the raytracing scheme, to determine whether changes in the wavenumber and wave amplitudes indeed agree with WKB theory. The results of this experiment are shown in Figures B1 and B2. Figure B1 shows a comparison of the wave amplitudes q i from Altitude (km) Wave amplitude growth factor ENDGame1D Ray tracing Figure B1. A comparison of the relative wave amplitudes q 0 from ENDGame1D and the wave amplitude growth factors W from the ray-tracing scheme. In ENDGame1D, there is a continuous oscillation to at the bottom boundary with a frequency of ω = 0.2 from which the wave amplitude is measured. In both cases, a uniform grid is used with z = 500 m and a centred time step is used with t = 1 s.

Altitude (km)
Wavenumber (× 10 −4 rad m −2 ) ENDGame1D Ray tracing Figure B2. A comparison of the wavenumbers m from ENDGame1D and the raytracing scheme at each height. In ENDGame1D, there is a continuous oscillation to at the bottom boundary with a frequency of ω = 0.2 from which the wavenumber is measured. In both cases, a uniform grid is used with z = 500 m and a centred time step is used with t = 1 s.
ENDGame1D to the wave amplitude growth factors W from the ray-tracing scheme, for an acoustic wave in cases where the wavenumbers and frequency being used are outside the limit of validity of WKB theory according to Figure 2. The two measures of the wave amplitude match up quite well in this plot, although the molecular viscosity appears to take slightly longer to reduce the wave amplitude in the thermosphere in ENDGame1D. Figure B2 shows a comparison of the wavenumbers m in ENDGame1D and the ray-tracing scheme. Again, the wavenumbers produced by both models are very similar, up to an altitude of just above 200 km. At this point, the wave amplitudes in ENDGame become smaller than the drift in the background and the wavenumber quickly falls to 0. However, the lines in Figures B1 and B2 are very similar, and it can be concluded that WKB theory can still be used to describe the behaviour of wave propagation for the small wavenumbers that are dealt with here.