Numerical modeling study of the momentum deposition of small amplitude gravity waves in the thermosphere

We study the momentum deposition in the thermosphere from the dissipation of small amplitude gravity waves (GWs) within a wave packet using a fully nonlinear two-dimensional compressible numerical model. The model solves the nonlinear propagation and dissipation of a GW packet from the stratosphere into the thermosphere with realistic molecular viscosity and thermal diffusivity for various Prandtl numbers. The numerical simulations are performed for GW packets with initial vertical wavelengths ( λz) ranging from 5 to 50 km. We show that λz decreases in time as a GW packet dissipates in the thermosphere, in agreement with the ray trace results of Vadas and Fritts (2005) (VF05). We also find good agreement for the peak height of the momentum flux ( zdiss) between our simulations and VF05 for GWs with initialλz ≤ 2πH in an isothermal, windless background, whereH is the density scale height. We also confirm that zdiss increases with increasing Prandtl number. We include eddy diffusion in the model, and find that the momentum deposition occurs at lower altitudes and has two separate peaks for GW packets with small initial λz. We also simulate GW packets in a non-isothermal atmosphere. The netλz profile is a competition between its decrease from viscosity and its increase from the increasing background temperature. We find that the wave packet disperses more in the non-isothermal atmosphere, and causes changes to the momentum flux andλz spectra at both early and late times for GW packets with initial λz ≥ 10 km. These effects are caused by the increase inT in the thermosphere, and the decrease in T near the mesopause.


Introduction
The propagation of gravity waves (GWs) in the Earth's atmosphere was first investigated by Hines (1960).Lower atmospheric disturbances such as deep convective events, tsunamis, typhoons, tropical cyclones, airflow over orography, cold fronts, shear and geostrophic adjustment are important sources of lower atmospheric GWs (Fritts and Alexander, 2003).
If a GW can survive critical level filtering, wave-breaking, and evanescence in the lower atmosphere, it can propagate into the thermosphere (Vadas andFritts, 2005, 2006, hereafter VF05 andVF06, respectively).Indeed, GWs have been observed in the thermosphere for decades, and are considered "ubiquitous" there (Oliver et al., 1997;Djuth et al., 2004).Because of the exponential decrease in the neutral density, a GW's amplitude increases exponentially with altitude.Therefore, GWs with larger vertical wavelengths (λ z ) and initially small amplitudes can have large amplitudes in the thermosphere.Only GWs with large λ z can propagate deep into the thermosphere (Vadas and Liu, 2009).Eventually, every GW dissipates in the thermosphere (Pitteway and Hines, 1963;Hines, 1973;Richmond, 1978;Hickey and Cole, 1988;Zhang and Yi, 2002;VF05;Walterscheid and Hickey, 2011;Vadas and Nicolls, 2012;Nicolls et al., 2012).
Molecular viscosity and thermal diffusivity have been incorporated into ray trace models for temporally and spatially variable GW sources (VF05; VF06; Vadas and Liu, 2009;Vadas and Crowley, 2010).If the spectral characteristics of the GW source and the background state are known, ray tracing can be a useful tool for studying the propagation of GWs and their influence on the background atmosphere (Jones, 1969;Schoeberl, 1985;Marks and Eckermann, 1995;Ding et al., 2003;VF05;VF06;Broutman and Eckermann, 2012;Vadas et al., 2009Vadas et al., , 2012)).Ray tracing is applicable for nonbreaking, small amplitude GWs which propagate in slowly varying background atmospheres (Einaudi and Hines, 1970;VF05).Its can also be applied to the large-amplitude GWs after properly including a scheme to parameterize the effects of GW saturation and/or breaking (Marks and Eckermann, 1995;Vadas and Crowley, 2010;Vadas et al., 2012).It is important to note that the theoretical basis for ray tracing is the WKB (Wentzel-Kramers-Brillouin) method, which requires the variation of the background atmosphere to be slow relative to the GW vertical wavelength (Pitteway and Hines, 1963;Einaudi and Hines, 1970;Nappo, 2002).For example, in a constant-wind (dissipative) thermosphere, the slowly varying background assumption requires λ z ≤ 4πH (VF05; Vadas, 2007).Here, H is the density scale height.
Using analytic solutions and numerical ray tracing, VF05 showed that λ z decreases as a GW dissipates in an isothermal atmosphere.For a Prandtl number of Pr = 1, analytic solutions showed that λ z decreased by a factor of 2 during dissipation.Vadas (2007) generalized this result for a nonisothermal atmosphere, and found that λ z can decrease, remain approximately constant, or slightly increase during dissipation, depending on the altitude and the background temperature profile.However, it was also found that the WKB approximation may not be valid H to 2H above the peak height of the momentum flux (z diss ) because the residue (Res) (defined by Einaudi and Hines, 1970) becomes > 1 there.It is therefore not known what the solution is at and above the altitude where Res > 1.However, VF05 and Vadas (2007) argued that this may not be important, because the GW momentum fluxes are typically negligible where Res > 1.
In order to determine if the VF05 solutions are valid during GW dissipation, a fully nonlinear numerical model solving the exact primitive fluid equations is needed in order to (1) determine the behavior of λ z as a GW (within a wave packet) dissipates, (2) compare the exact solutions with the Vf05 solutions, and (3) determine under what conditions the WKB assumptions (used in the VF05 solutions) are valid.
Fully nonlinear numerical simulations of GW packets with quasi-monochromatic spectra were performed by Zhang and Yi (2002) to study the propagation of temporally and spatially localized GW packets in a dissipative thermosphere.They showed that λ z decreases as the GW packet dissipates, because of the vertical inhomogeneity of molecular viscosity.Additionally, they showed that their results differed from the results of a ray trace model.However, their comparison was inadequate to rule out the soundness of the WKB ray tracing method in the thermosphere, because their ray trace model was non-dissipative.Additionally, only GW packets with initial λ z < 16 km were examined; thus, their results are not generalizable to larger initial λ z .
The main objective of this paper is to determine the temporal evolution of the momentum flux profiles of small amplitude, dissipating GWs within a wave packet using a twodimensional, fully nonlinear, compressible, numerical model with molecular viscosity and thermal diffusivity (Liu et al., 2008(Liu et al., , 2009)).Our results include various Prandtl numbers, eddy diffusion, and isothermal and non-isothermal backgrounds.We then compare our results to the VF05 solutions, in order to better understand the validity of the WKB ray trace approach.
Our paper is arranged as follows.Our numerical model is described in Sect. 2. In Sect.3, we simulate the GW momentum deposition in isothermal and non-isothermal atmospheres which include molecular viscosity and eddy diffusion.We also compare our results with the VF05 ray tracing model solutions.Our conclusions are contained in Sect. 4.

Governing equations and calculations
The governing equations solved in our model are the twodimensional compressible, non-linear Navier-Stokes equations and the equation of state for an ideal gas (Batchelor, 1967), Where x and z are the horizontal and vertical coordinates, g is the acceleration of gravity, ρ and p are the density and pressure of atmosphere, u and w are the horizontal and vertical velocities, respectively, T is the temperature, c v = 718 J kg −1 K −1 is the specific heat at constant volume, γ = 1.4 is the ratio of specific heat, and R = 287 J kg −1 K −1 is the gas constant for dry air.In the equations, u , w and T are the time-dependent fluctuations associated with the GWs, and are calculated by subtracting the average values over one horizontal wavelength from the total values (Xu et al., 2003).
The Rayleigh friction α(z) is adapted from Xu et al. (2003), and is utilized to reduce or eliminate spurious reflection from the lower and upper boundaries.Note that the above equations are slightly different from those used in Liu et al. (2008Liu et al. ( , 2009) ) with respect to the diffusion terms in Eqs.
The ratio between the molecular viscosity ν and thermal diffusivity κ is defined as the Prandtl number Pr.The kinematic molecular viscosity ν is adopted from Banks and Kockarts (1973), ν(z) = 3.5 × 10 −7 T (z) 0.69 /ρ(z). (6) It has almost the same profile as that used by Yu et al. (2009) and Liu et al. (2010) in their local or global models.Its profile  (min) 11.72 11.76 11.91 12.16 12.51 12.95 is shown as the solid line in Fig. 1a.The thermal diffusivity κ(z) is calculated from Pr and ν(z).To directly compare with the results in VF05, we choose Pr = 0.7, 1.0 and ∞ in this study.
The vertical extent of our model is 300 km including a sponge layer of ∼ 40 km at the upper boundary.The horizontal domain employs periodical boundary conditions.Here, we set the horizontal domain to equal one GW horizontal wavelength.These values are listed in Table 1.The grid sizes in the vertical and horizontal directions are 0.2 km and λ x /20 km, respectively.Due to the large molecular viscosity in the thermosphere, the time step is set to 5 × 10 −4 s.

Background atmosphere and initial GW perturbations
Two atmospheric background conditions are investigated in this study.First a windless, isothermal atmosphere with a constant temperature of 239 K (shown as a dash line in Fig. 2) is used in order to compare directly to the results in VF05.Thus, the scale height H = RT/g is ∼ 7 km and the buoyancy frequency N = {(g/T ) × [(dT /dz) + g/c p ]} 1/2 is ∼ 0.02 rad s −1 .These values are consistent with those used in VF05.The second background state is a non-isothermal but windless atmosphere calculated from NRLMSIS-00 (Picone et al., 2002) at 45 • N on day 173 (shown as a solid line in Fig. 2).This is designed to study the effects of a nonisothermal background temperature on GW momentum deposition.
The initial horizontal velocity perturbation is similar to that used by Xu et al. (2003) and Liu et al. (2008), To exclude nonlinear processes such as wave-wave interactions, wave-flow interactions, wave saturation/breaking, etc., the initial amplitude is chosen to be small in all cases (A = 1.0 × 10 −3 m s −1 ).The horizontal and vertical wave numbers are k x = 2π/λ x and k z = −2π/λ z , respectively.We set z 0 = 60 km, which is the peak height of the wave packet energy; therefore, the wave packet is launched from ∼ 60 km.
The width of the wave packet in the vertical direction is proportional to the vertical wavelength, σ z = λ z .The other initial perturbations are derived from u (x, z, t = 0) in Eq. ( 6) using the polarization equations for a linear GW in a nondissipative atmosphere (Fritts and Alexander, 2003), since molecular viscosity is unimportant at z 0 = 60 km.The ini-Fig. 3. Snapshots of GW B with initial λ z = 10 km (see Table 1) propagating in the isothermal atmosphere with Pr = 1.We show the square root of the density-weighted horizontal wind disturbances (m s −1 ) at times of 0, 37.5, 75, 112.5, and 150 min.
tial GW vertical and horizontal wavelengths and wave periods are listed in Table 1.The wavelengths of GWs A, B, C, and F are purposely chosen to be the same as those used in VF05.It is important to note that VF05 used monochromatic GWs in their ray trace study, while GW packets having different spectral parts are used here.In an isothermal atmosphere, there is only a small difference in the results for a wave packet and a monochromatic GW because of the small dispersion of the wave packet (see Sect. 3.3).In a nonisothermal atmosphere, the dispersion of a wave packet and the variations of λ z are much larger (see Sect. 3.5).

GW propagation in the thermosphere
In this subsection, we utilize GW B as an example to illustrate the propagation and dissipation of GW packets within an isothermal atmosphere.Here, the Prandtl number is set to 1. Figure 3 displays snapshots of the GWs in square root density-weighted horizontal wind disturbances u (x, z) = [ρ 0 (z)/ρ 0 (z r )] 1/2 u (x, z) at times of 0, 37.5, 75, 112.5, and 150 min, where ρ 0 (z r ) is the background density at a reference height of 60 km (Huang et al., 2010).Figure 3 shows that u (x, z) does not change significantly with altitude and time in the absence of dissipation below the height of 100 km and before 75 min.This is because wave energy is conserved before the GW packet enters the dissipative thermosphere.After the GW packet propagates upward into the lower thermosphere (120-140 km), the amplitude of u (x, z) decreases rapidly.At t = 150 min, the GW packet is almost completely dissipated by molecular viscosity and thermal diffusion.
It is important to note that during this dissipative process, the GW phase lines do not become vertical; in fact, the slope of the phase lines are similar before and during dissipation.This shows that λ z does not change significantly during dissipation since the slope of a phase line equals −k z /k x = λ x /λ z .Therefore, qualitatively, the predictions of VF05 are correct for this wave packet.
The GW momentum flux and its corresponding body force are important quantities to illustrate the effect of the GWs on the mean flow.They are defined as (Andrews et al., 1987), and The time-height cross-section of the momentum flux is shown in Fig. 4a.One can see that most of the momentum flux carried by the GW packet is concentrated in a layer at z ∼ 120-140 km, with a peak at z ∼ 128 km (shown as a white dash line) and at a time interval of t ∼ 75-130 min.Before ∼ 75 min, the wave packet is located far below z ∼ 128 km, and has an exponentially smaller momentum flux because the neutral density is larger there.After t ∼ 130 min, the GW packet has been dissipated by the large molecular viscosity and thermal diffusivity in the thermosphere.All of this momentum has been transferred into the background fluid.
To obtain the vertical wavelength power spectra, a discrete Fourier transform is first applied at each time step to the vertical profile of u (x mid , z), where x mid denotes the mid-point of the horizontal computational domain.Then the power spectrum of each vertical wavelength is normalized by the maximum power obtained from u (x mid , z); we refer to this spectrum as the normalized vertical wavelength power spectrum.Note that a limited spatial coverage of a signal leads to an uncertainty of the wave number power spectrum.For example, if the uncertainty in resolving λ z is zero, then the spatial coverage of the signal must be infinite (Folland and Sitaram, 1997).Thus, for the fixed computational vertical domain used in our model (300 km), the uncertainty in resolving λ z increases with the increasing λ z .This will be further analyzed in Sect.3.3.Figure 4b shows the temporal evolution of the normalized vertical wavelength power spectrum.We find that the vertical wavelength decreases from λ z ∼ 10 km at t ∼ 70 min to about λ z ∼ 8.5 km at t ∼ 150 min.
Comparing Fig. 4a and b, we see that λ z begins to decrease when the GW packet reaches the altitude where it begins to dissipate (i.e., at z ∼ 128 km).Additionally, λ z continues to decrease during and after the wave packet dissipates.Thus, λ z decreases during dissipation because of the increasing molecular viscosity and thermal diffusivity in the thermosphere.During dissipation, the momentum flux within the GW packet is deposited into the background fluid.
For the special case of Pr = 1, with the help of Eq. ( 62) of VF05, the variation of k z with time can be written as, For an upward propagating GW with negative k z , Eq. ( 10) shows theoretically that k z becomes increasingly negative because of the gradient of the viscosity in the thermosphere.As a result, λ z decreases after the GW encounters large viscosity in the thermosphere.If a monochromatic GW is replaced by a wave packet in the VF05 ray trace model, λ z also decreases in altitude because of the increasing viscosity.In order to explore how λ z changes during dissipation for monochromatic GWs, we changed the spectral width of our initial GW source by changing the width of the wave packet; indeed, we found the same results that λ z decreases as the wave packet dissipates (not shown here).Thus, we find good agreements between our numerical model results and the VF05 ray trace solutions.

GW momentum fluxes and body forces for differing GW packets
We now simulate the propagation and dissipation of GWs C, D, and F from Table 1 in the isothermal atmosphere (dash line in Fig. 2) with Pr = 1.These GWs have initial vertical wavelengths of 20, 30, and 50 km.corresponding momentum flux profiles, and the normalized vertical wavelength power spectra.This figure is similar to Fig. 4 in that (1) the momentum fluxes are concentrated in altitude and time, and (2) λ z decreases while each GW packet dissipates.This decrease is larger for the GWs with smaller λ z ; it is ∼ 25 % for the GW packet with initial λ z = 20 km, and is ∼ 15 % for the GW packet with initial λ z = 30 km.We see that there is numerical noise, especially for the GW packet with initial λ z = 50 km.This noise adds additional uncertainty to the peak values of λ z inferred for these GWs.Additionally, we note from Fig. 5 that the momentum fluxes are concentrated at increasingly higher altitudes as the initial λ z increases.
For a small amplitude GW in a windless environment, changes in its momentum flux profile are caused by molecular viscosity and thermal diffusivity.As in VF05, we define z diss to be the altitude where the momentum flux per unit mass, MF(z), is maximum.For our wave packet simulation, we define MF(z) as, where t 0 and t 1 are the initial and final times for the simulations, respectively.In order to determine an accurate mo-mentum flux profile, the wave must dissipate prior to t 1 .Figure 6 shows the MF(z) profiles for the GW packets with initial λ z ranging from 5 to 50 km (solid lines).All GW packets have the same initial amplitudes.We see that GW packets with larger initial λ z dissipate at higher altitudes and achieve larger amplitudes at z diss .This agrees qualitatively with the results of VF05.The corresponding body force profiles are also shown in Fig. 6 (dash lines).For all of the GW packets shown here, the maximum body force occurs ∼ 4.5-5.0km (i.e., ∼ 0.6-0.7H ) abovez diss ; this result quantitatively agrees with VF05, since the peak of the body force is higher than z diss by ln(2)H (Eqs.54 and 55 of VF05).However, the body force profiles created by the larger λ z GW packets are not as narrow as that predicted by the ray tracing model of VF05.This discrepancy may be attributable, in part, to the difference between the monochromatic GW ray traced in VF05, and the Gaussianshaped wave packets with different spectra simulated here.

GW momentum fluxes for differing Prandtl numbers
To study the effects of different Prandtl numbers on the GW momentum fluxes, we perform simulations for GW packets with different initial λ z in an isothermal atmosphere with Pr = 0.7 and ∞.We combine these results with those from The Prandtl numbers for Fig. 6a, b and c are 0.7, 1.0 and ∞, respectively.We also overlap the VF05 results (taken from Table 1 of VF05) (square, ).
Sect. 3.1.The values of z diss for Pr = 0.7, 1.0 and ∞ are shown in Fig. 7a, b, and c, respectively, as crosses (×).In addition, the values of z diss listed in Table 1 of VF05 are also plotted in Fig. 7 (as squares, ) for comparison.Figure 7 indicates that: 1.The dissipation heights calculated from our simulations agree well with the ray trace results from VF05, especially for GW packets with initial λ z = 5 km, 10 km and 20 km.VF05 did not calculate z diss for initial λ z = 30 and 40 km.For GW packet with initial λ z = 50 km, the dissipation height calculated from our simulation is slightly higher than that from VF05, independent of Prandtl number.We discuss this difference further in the next section.
2. Our simulation results and those of VF05 show that GW packets in a Pr = 1.0 environment dissipate at slightly higher altitudes than those in a Pr = 0.7 environment (by ∼ 1 km).This is because the thermal diffusivity is slightly larger if Pr is smaller.More generally, our simulation results and those of VF05 show that the dissipation heights of GW packets increase as Pr increases.This is because the Prandtl number (ν/κ) is inversely proportional to the thermal diffusivity.For fixed ν, larger Pr corresponds to smaller thermal diffusivity and smaller loss of wave energy at a given altitude, thereby allowing GW packets to propagate to higher altitudes before dissipating.

Discussion of the comparison between our results and the VF05 ray trace
The good agreement obtained between our simulation results and VF05 in Sects. of GW packets from temporally and spatially varying GW sources.The slight discrepancies that exist between our simulation results and VF05 for GW packets with λ z = 50 km can be attributed both to the uncertainty in resolving λ z in a limited vertical domain, numerical noise, and the small change of λ z in time (at a fixed altitude in the thermosphere) due to the dispersion of different spectral parts of the GW packet.Equation ( 73) of VF05 provides a reasonable way to estimate the uncertainty of the dissipation height.Taking the derivative of this equation with respect to k z , we obtain As previously mentioned, since the vertical extent of our model is 300 km, the uncertainty of resolving k z increases with increasing initial λ z (e.g., Folland and Sitaram, 1997).
For the GW packet with initial λ z = 50 km (5 km), there are ∼ 6 (60) wave cycles in the vertical domain.If we estimate the maximum error in resolving k z to be 1 wave cycle, then the maximum relative error of k z /|k z | is ∼ 1/6 (1/60).Then, Eq. ( 12) implies that the uncertainty in z diss is ∼ 3.5 km and 0.35 km for the GWs with initial λ z = 50 km and 5 km, respectively.Additionally, this theoretical estimation can be proved from Fig. 5. From Fig. 5f, we found that the error in determining the mean λ z is larger for the GW with initial λ z = 50 km.While for the GWs with initial λ z = 20 and 30 km (Fig. 5b and d), these errors are small.Finally, wave dispersion causes slight differences in the propagation times and dissipation altitudes of the GWs within the packet.This occurs because the wave packet has different spectral parts with slightly different vertical wavelengths and group velocities.As a result, the vertical wavelength and dissipation altitude may change with time.This effect causes an error in calculating the mean λ z and z diss .We assume that the initial wave packet has a central vertical wavelength of λ z and a spectral width of 1/λ z = |k z |/(2π ).For the GW packet with initial λ z = 50 km, if we assume that one half of the wave energy is attributable to the central wave number k z and one fourth to the "side" wave numbers then there is an uncertainty of in resolving the dissipation height of the GW packet.Combining these effects, the total uncertainty in calculating z diss for the GW packet with initial λ z = 50 km is as large as ∼ 5-6 km.Since the difference between our simulation results and VF05 is less than 5 km for the GW with initial λ z = 50 km in an atmosphere with Pr = 0.7, 1.0 and ∞. (see Fig. 7), we conclude that this difference is not significant within the error bars.Therefore, taking into account uncertainties in calculating z diss , we conclude that our simulation results agree well with those of VF05.
In order to derive an analytic form for the GW dissipative dispersion relation, VF05 assumed that ν was locally constant during GW dissipation.This assumption was shown to be approximately equivalent to λ z ≤ 4π H at the altitude where the GW dissipates (Vadas, 2007).Our numerical results show that this assumption holds for GWs with λ z ≤ 2π H at the dissipation altitude (i.e., λ z ≤ 40-50 km, for H = 7 km).

Effects of eddy diffusion
In the mesopause region, turbulence can be quite large, and can be confined to a relatively shallow layer mainly because of GWs breaking (Lindzen, 1981;Balsley et al., 1983;Hocking, 1987;Lübken, 1997).Here a simple Gaussian function is chosen to represent eddy diffusion (unit in m 2 s −1 ) (Xu et al., 2009): The eddy diffusion coefficient η(z) is added directly to the molecular viscosity ν(z) in Eq. ( 2) and to the thermal diffusivity κ(z) in Eq. ( 3).Here, we set Pr = 1 and η max = 100, 150, 200 m 2 s −1 , which are similar to previous study values of 10-1000 m 2 s −1 (Hocking, 1987;Xu et al., 2009).The total diffusion which includes both molecular diffusion and eddy diffusion is shown in Fig. 1.Including eddy diffusion results in a sharp increase of the total diffusion at z ∼ 75-100 km, and in a sharp decrease at z ∼ 100-110 km.Above 110 km, molecular viscosity dominates, and increases approximately exponentially with altitude.
Because GWs with smaller λ z are more sensitive to eddy diffusion than those with larger λ z due to the scale dependence of eddy diffusion (Xu et al., 2009;Walterscheid and Hickey, 2011), we simulate those GW packets with λ z = 5 km and λ z = 10 km here, in order to illustrate the effects of the total diffusion on the GW momentum flux profiles.Figure 8a-d shows the temporal evolution of the momentum flux profiles for the GW packet with initial λ z = 5 km propagating in the isothermal atmosphere with Pr = 1 and η max = 0, 100, 150, 200 m 2 s −1 , respectively.Whereas the momentum flux is localized in a single, relatively thin height layer when η max = 0 m 2 s −1 , eddy diffusion causes the momentum flux to be spread out in altitude.In particular, there are two peaks for η max = 150 m 2 s −1 .One is located at z ∼ 114 km, consistent with z diss obtained for the case without eddy diffusion.The other is located at z ∼ 95 km, and is caused by eddy diffusion since the magnitude of the total diffusion at z ∼ 95 km is equal to the molecular viscosity at z ∼ 114 km.For η max = 100 m 2 s −1 (200 m 2 s −1 ), most of the GW's momentum flux is concentrated in a layer centered at z ∼ 114 km (95 km).These results indicate that the momentum flux profile, and the deposition of momentum, caused by a GW packet with initial λ z = 5 km is sensitive to eddy diffusion in addition to molecular viscosity and thermal diffusivity.
As for the case of the GW packet with initial λ z = 10 km (figures not shown here), although the peak of the maximum   and (b, d, f) normalized vertical wavelength power spectra for the cases of (left column) GW A, (middle column) GW B and (right column) GW C propagating in the non-isothermal atmosphere.The white dash lines show the heights of z diss .The scaling factors are shown at the right of the color bar.momentum flux decreases as eddy diffusion increases, z diss does not change for this GW packet.Thus the momentum deposition of this GW packet is less sensitive to the additional eddy diffusion than for the GW packet with initial λ z = 5 km.This is because the impact of eddy diffusion on the momentum flux is scaled by k 2 z , and is less important for GWs with larger λ z (Eq.5d of Xu et al., 2009;Eq. 13 of Walterscheid and Hickey, 2011).

Effects of a non-isothermal background temperature
Figure 9 shows the temporal evolution of the momentum flux (upper row) and normalized vertical wavelength power spectra (lower row) for GW cases A, B and C (see Table 1) propagating in the non-isothermal (solid line in Fig. 2) and windless atmosphere with Pr = 1.The white lines in the upper row show the heights of z diss .It can be seen that λ z www.ann-geophys.net/31/1/2013/Ann.Geophys., 31, 1-14, 2013 decreases before the GWs enter the lower thermosphere, unlike in the isothermal atmosphere (see Fig. 4).Since viscosity below the mesopause region is negligible, the decrease of λ z in the region z ∼ 80-100 km is mainly due to the decreased background temperature in the height range of ∼ 60-90 km (Vadas, 2007;Zhou and Morton, 2007;Fig. 5 of Liu et al., 2012).In particular, the smaller T creates a larger buoyancy frequency N .Since λ z ∝ 1/N, this decreases the GW's vertical wavelength.
Comparing with Figs.4-5, we see that the GWs with initial λ z = 10 and 20 km dissipate at lower altitudes than the same GW packets propagating in the isothermal atmosphere.At first glance, this appears to be in contrast with Vadas (2007), who finds that GW propagates to higher altitudes when the temperature in the thermosphere is larger.However, that general rule applies to GW with relatively large λ z that can propagate to high-enough altitudes in the thermosphere where λ z is able to increase enough (from increasing T ) to compensate for its decrease near the mesopause.In fact, Vadas (2007) shows that a GW must dissipate at z diss > 175 km in order to fully recover from the decrease in λ z near the mesopause; otherwise, z diss may be lower in the non-isothermal atmosphere (Eq.73 of VF05; Vadas, 2007).
An unusual feature in Fig. 9 is that a second peak occurs in the momentum flux profile and λ z spectra at early times; this is especially noticeable for GW C (with initial λ z = 20 km), although it also occurs to a lesser extent for GW B (with initial λ z = 10 km).For GW C, the 2nd peak occurs at t = 25-30 min, z ∼ 140-150 km, and λ z ∼ 20-30 km.The occurrence of the second peak in the momentum flux profile may occur because some of the fast GWs at the leading edge of the packet are more affected by the increasing thermospheric temperature.To illustrate this, we show the horizontal wind perturbations in Fig. 10 for GWs B and C in the isothermal atmosphere, and in Fig. 11 for GWs B and C in the non-isothermal atmosphere.We see wave dispersion occurring within each wave packet, especially in the non-isothermal atmosphere.This is especially noticeable for GW C. Within each packet, it is clear that the fastest (highest in altitude) GWs have the largest λ z , since the slope of a GW's phase line is proportional to λ z .As these GWs enter the thermosphere, the increasing background temperature causes λ z to increase because of the decreasing buoyancy frequency (VF06).Larger vertical wavelengths result in larger vertical group velocities, and therefore faster propagating times (Fritts and Alexander, 2003).Larger vertical wavelengths also result in these fast GWs achieving higher altitudes prior to dissipating (Vadas, 2007).This effect is likely responsible for the second peak in the vertical wavelength power spectrum in Fig. 9.
We now summarize how the momentum deposited by GWs propagating in the non-isothermal background is different from that in the isothermal background.First, the momentum deposition height is lower than that in the isothermal cases for GWs A, B, and C.This occurs because the increase in λ z in the thermosphere is not enough to offset the decrease in λ z in the mesopause region, because z diss is not high enough (Vadas, 2007;Zhou and Morton, 2007;Liu et al., 2012).Second, the momentum deposition height increases as the initial λ z increases for both the isothermal and non-isothermal atmospheres.This occurs because GWs with larger λ z have larger horizontal phase speeds and group velocities, and can therefore propagate to higher altitudes (where ν is larger) before dissipating.
Moreover, we found that during GW dissipation and momentum deposition, the variation of λ z depends on the temperature profile in the thermosphere.During dissipation in the non-isothermal atmosphere, λ z decreases significantly for GW A, decreases somewhat for GW B, and is relatively constant for GW C. The results are different from the isothermal results, for which λ z decreases significantly during dissipation (see Figs. 4-5 for GWs B and C).From Fig. 10, z diss = 101, 113, and 128 km for GWs A, B, and C, respectively.From Fig. 2, the background temperature is still approximately constant at z = 101 km; this causes GW A to dissipate in an approximately isothermal background, thereby causing λ z to decrease significantly during dissipation.For GWs B and C, the background temperature is increasing substantially where the GWs dissipate.Since larger temperatures yields larger λ z (if the GW is not dissipating), there is therefore a competition between the decrease in λ z because of GW dissipation, and the increase of λ z because of the increasing background temperature.This yields the result that λ z is relatively constant during dissipation for these GWs.A similar result wave obtained using the ray tracing model (see Fig. 4 in Vadas and Nicolls, 2012).

Conclusions
In this paper, we simulated the propagation of spatially and temporally localized GW packets from the upper stratosphere to the thermosphere using a fully nonlinear, twodimensional, compressible, numerical model, which includes eddy diffusion in the mesosphere and molecular viscosity and thermal diffusivity in the thermosphere.These simula-tions were performed for GW packets with small initial amplitudes, in order to ensure that the GWs evolved linearly.Additionally, the GWs had initial vertical wavelengths of λ z ≤ 50 km.The background is windless in order to simplify the interpretation of the results, and included both isothermal and non-isothermal temperature profiles.We also employed varying Prandtl numbers.
We first simulated a GW packet with initial λ z = 10 km in an isothermal atmosphere with only molecular viscosity and thermal diffusivity.We found that after the GW packet entered the thermosphere, their amplitudes decreased significantly in time, and became negligible within an hour.Importantly, the slope of the phase lines did not appear to change substantially during GW packet dissipation.An analysis of the vertical wavelength power spectrum was performed.We found that the peak λ z decreased, not increased, during GW packet dissipation.This result agrees well with the ray tracing solutions of VF05, but disagrees with the steady-state results of Walterscheid and Hickey (2011) (for which λ z increases exponentially during dissipation up to a very large value).
Next, we simulated GW packets with initial λ z varying from 5 to 50 km in the isothermal atmosphere.We found that the altitude at which dissipation occurs depends on the vertical wavelengths of the GW packets, with larger-λ z GWs dissipating at higher altitudes than smaller-λ z GWs.For all of the GW packets, λ z decreased during GW dissipation, although less so for the GW packets with larger values of λ z , this is in qualitative agreement with VF05.We then compared our simulation results quantitatively with the results of the ray tracing model in VF05 for high-frequency GWs with initial λ z ≤ 50 km.We found very good agreement between our results and VF05 for the heights of the maximum momentum fluxes and body forces (created by the dissipation of GW packets), this agreement includes their sensitivity to varying Prandtl numbers.For the GW packets modeled here, then, such agreement validates the use of the VF05 ray tracing model for calculating the dissipation of GW packets within the thermosphere.
We then simulated GW packets using molecular viscosity and thermal diffusivity in the thermosphere and eddy diffusion near the mesopause.For GW packet with initial λ z = 5 km, inclusion of eddy diffusion created significant dissipation below the turbopause.For certain eddy diffusion amplitudes, a second momentum flux layer was created.This phenomenon did not occur for the GW packet with initial λ z = 10 km.These results confirm previous findings that the dissipation of GWs from eddy diffusion is vertical wavelength dependent.This also demonstrates our model's ability to resolve the effects of eddy diffusion on the momentum deposition of GW packets.
Moreover, we simulated the propagation of GW packets through a non-isothermal background temperature profile.We found that λ z decreases significantly as the GW packets propagate through the mesosphere, because the buoyancy frequency increases there.This confirms ray tracing results (Vadas, 2007).We also found that λ z continues to decrease (or stays relatively constant) above the turbopause for all of the GW packets simulated here.This occurs because all of the GW packets simulated in this paper dissipate within a few density scale heights above the altitude where the temperature begins to increase; therefore the λ z profiles are influenced by both the increasing thermospheric temperature (which causes λ z to increase), and dissipation from the rapidly increasing molecular viscosity (which causes λ z to decrease).The competing effects determine the variation of λ z during dissipation.Note that such competing effects have been seen using more-idealized temperature profiles with the VF05 ray tracing model (see Fig. 3 of Vadas, 2007); indeed, they found that λ z can only significantly increase in the thermosphere if λ z is large-enough near the turbopause, and if the thermosphere is hot enough.
Because λ z decreased significantly in the mesosphere, and because the GWs simulated here dissipated near the turbopause, we find that the momentum deposition altitudes in the thermosphere in the non-isothermal cases are lower than in the isothermal (lower-temperature) cases.This result seems to be inconsistent with Vadas (2007), who found that GWs propagate to higher altitudes in a hotter thermosphere.However, Vadas (2007) found that a GW propagates to higher altitudes in a hotter thermosphere only for those GWs with λ z > ∼ 60 km and horizontal phase speeds of c px > ∼ 60 m s −1 (see Fig. 4a, c, e, g in that paper).Such larger-λ z GWs have dissipation altitudes of z diss > 175 km, which are much higher than the dissipation altitudes of the GWs simulated here.Additionally, all of the GWs simulated here have initial λ z ≤ 50 km, which are smaller than those threshold values.
Finally, we showed in detail how the GW packets disperse in altitude and time in both the isothermal and nonisothermal atmospheres for the GWs with initial λ z = 10 and 20 km.We showed that those GW spectra with slightly larger λ z (within a packet) reached the thermosphere earlier and dissipated at higher altitudes than the mean-λ z GWs.
We emphasize that all of the simulations performed in this paper were windless.This was done in order to compare our results directly with VF05 and to obtain robust, simple understandings of the effects of dissipation and temperature on localized GW packets.However, the winds in the thermosphere are large, and significantly influence λ z , GW propagation, and GW dissipation (e.g., Fritts and Vadas, 2008).We will investigate the effects of wind and wind shear on GW wave packets in a future paper.

Fig. 4 .
Fig. 4. Temporal evolutions of (a) momentum flux and (b) normalized vertical wavelength power spectra for the GW shown in Fig. 3.

Fig. 8 .
Fig. 8. Temporal evolutions of the momentum flux profiles for the GW with initial λ z = 5 km propagating through the isothermal atmosphere (shown as the solid line in Fig. 2) with Pr = 1 and η max = 0, 100, 150, 200 m 2 s −1 , respectively.The scaling factors increase from left to right.Dotted lines show the altitude z = 114 km, which is the momentum deposition height for the case without eddy diffusion.

Fig. 9 .
Fig. 9. Temporal evolutions of (a, c, e) momentum flux (unit is m 2 s −2 ) and (b, d, f) normalized vertical wavelength power spectra for the cases of (left column) GW A, (middle column) GW B and (right column) GW C propagating in the non-isothermal atmosphere.The white dash lines show the heights of z diss .The scaling factors are shown at the right of the color bar.

Fig. 10 .Fig. 11 .
Fig. 10.(a) Horizontal wind perturbations for GW with initial λ z = 10 km propagating in the isothermal atmosphere with Pr = 1 at various times (labeled).(b) Same as (a), but for GW with initial λ z = 20 km.The times are different than in (a).

Table 1 .
The wavelengths and periods of the initial GWs.