Geophysicae Reconstruction of the gravity wave field from convective plumes via ray tracing

We implement gravity wave (GW) phases into our convective plume and anelastic ray trace models. This allows us to successfully reconstruct the GW velocity, temperature, and density perturbation amplitudes and phases in the Mesosphere-Lower-Thermosphere (MLT) via ray tracing (in real space) those GWs that are excited from a deep convective plume. We find that the ray trace solutions agree very well with the exact, isothermal, zero-wind, Fourier-Laplace solutions in the Boussinesq limit. This comparison also allows us to determine the normalization factor which converts the GW spectral amplitudes to real-space amplitudes in the ray trace model. This normalization factor can then be used for ray tracing GWs through varying temperature and wind profiles. We show that by adding GW reflection off the Earth’s surface, the resulting GW spectrum has more power at larger vertical and horizontal wavelengths. We determine the form of the momentum flux and velocity spectra which allows for easy calculation of GW amplitudes in the MLT and thermosphere. Finally, we find that the reconstructed (ray traced) solution for a deep, convective plume with a duration much shorter than the buoyancy period does not equal the Fourier-Laplace Boussinesq solution; this is likely due to errors in the Boussinesq dispersion relation for very high frequency GWs.


Introduction
The Spread F Experiment (SpreadFEx) was performed in Brazil in 2005, as summarized by Fritts et al.(2008a).Key to attainment of SpreadFEx goals were 1) measuring GW Correspondence to: S. L. Vadas (vasha@cora.nwra.com)scales and amplitudes in MLT, thermosphere and ionosphere (TI), 2) linking these measured GWs to convective sources at lower altitudes, and 3) assessing the potential contributions of these GWs to plasma instability dynamics at the bottomside F layer.A central need underlying the connection from tropical convection to GW scales and amplitudes at higher altitudes is an ability to link convection scales and intensities directly to observed GW scales and amplitudes at some altitude where they can be assessed directly.Previous and current studies by Vadas and Fritts (2004) (hereafter VF2004) and Vadas and Fritts (2006) (hereafter VF2006) anticipated specific GW responses to individual convective plumes and/or convective complexes with some assumptions.However, no studies to date have quantified this link with sufficient confidence to extrapolate GW amplitudes to high altitudes at spatial scales not directly measured in the MLT.Given our current measurement capabilities, two steps are needed.The first step is to relate GW amplitudes at the largest observed spatial scales to the characteristics of the convective plumes from which they arose.This is accomplished by reverse ray tracing these medium-scale GWs observed in MLT to specific convective cells identified and quantified from satellite imagery (Vadas et al., 2009a).This constrains as fully as possible the convective spatial and temporal scales important during SpreadFEx.
The second step in linking convection to GWs at the bottomside of the F layer is to calculate the GW amplitudes in the MLT expected to arise for these same convective sources, but for spatial scales not observed easily or directly in the MLT because of steep propagation angles and obscuration by anvil clouds (Vadas et al., 2009a).This is one of the goals of this paper, to quantify the GW source spectra from convection in a form that allows for easy calculation of GW amplitudes at larger spatial scales in the MLT and above.In order to accomplish this goal, it is necessary to incorporate GW phases into the convective source and ray trace models in order to reconstruct the GW fields in the MLT and TI. S. L. Vadas and D. C. Fritts: Reconstruction of the gravity wave field We then compare the ray trace solutions with the Fourier-Laplace (FL) Boussinesq solutions for zero wind, isothermal conditions.This inter-comparison allows for the normalization of the GW spectral amplitudes, so that the amplitudes of GWs can be determined via ray tracing.It also allows for studies involving ray tracing through varying temperature and wind profiles (Vadas et al., 2009b).Note that Fritts and Vadas (2008) employ the GW amplitudes and spectra determined in this paper to anticipate general GW propagation, filtering, and amplitude growth accompanying various thermospheric temperature and wind profiles.Fritts et al. (2008b) also employ these results to infer specific neutral and plasma responses at the bottomside F layer.
Our paper is structured as follows.Section 2 describes the model which calculates the amplitudes, scales and phases of GWs excited from a convective plume with ground reflection, and the model which ray traces these GWs into the MLT and TI.A description of the calculation of GW momentum fluxes and phases is also included.Reconstruction of the GW field via ray tracing, and an inter-comparison of the FL and ray trace solutions in a windless, isothermal atmosphere for a deep, convective plume is presented in Sect.3. We also determine the normalization factor for the spectral wave amplitudes.In Sect.4, we use the results from ray tracing a few, individual GWs to compute a simple, spectral form for the GW momentum flux and horizontal wind spectra at OH airglow altitudes.This allows for easy calculation of GW amplitudes at higher altitudes.Section 5 provides an intercomparison of the FL solutions with the ray trace solution for a deep, convective plume when its duration is much shorter than the buoyancy period.Our conclusions are provided in Sect.6.

Models utilized
Our convective plume model is primarily utilized to calculate the spectral amplitudes and phases of the GWs excited from a deep, convective plume.However, if Fourier-transformed to real space, this model can also be used to calculate the Boussinesq, real-space wave field at any altitude and time, assuming isothermal, zero-wind conditions.This latter feature is essential for normalizing the GW spectral amplitudes used for ray tracing.Our ray trace model inputs the spectral amplitudes and phases from the convective plume model, and ray traces these GWs into the MLT and TI through varying temperatures and winds.We now describe these models.

Convective plume model
Convective sources of GWs can be described equivalently as heating or momentum sources, as these sources are coupled through the vertical momentum equation.One such heating model which describes the spectrum of GWs from a deep, convective plume is given by Walterscheid et al. (2001).Our convective plume model is based on the use of a vertical body force(s) (Vadas and Fritts, 2001) (hereafter VF2001) to describe the uplifting of air which occurs in a convective plume (VF2004; VF2006).This model is linear, Boussinesq, neglects moisture processes, and assumes that the air above the tropopause is stationary within the frame of the mean horizontal wind at the tropopause (U trop ), until a convective plume overshoots the tropopause and pushes the stratospheric air upwards.It solves the linear solutions in a locally unsheared environment with a constant buoyancy frequency.Observations and simulations show that there are typically many small updrafts within the "envelope" of a convectively unstable region, which give rise to a GW spectrum concentrated at small-scales of ∼5-10 km (e.g., Larsen et al., 1982;Alexander et al., 1995).Our model neglects the individual updrafts which generate these small-scale GWs, as these GWs are not likely to propagate to the upper atmosphere and thermosphere (due to wave breaking, critical level absorption, and reflection in the stratosphere).Instead, our model calculates the spectrum of larger-scale GWs excited by the "envelope" of the upward motion of air within a convective plume.(These are the larger-scale GWs, which are more important in the mesosphere and thermosphere.)Each plume is described by its geometry (horizontal width and vertical depth), its duration, its maximum vertical updraft velocity, and its latitude, longitude, and time of occurrence.
This model calculates the amplitudes, scales, and phases of the excited GWs.One of the limitations of this linear model is that the wind is assumed constant in the convective plume region in order to obtain an analytic solution.Although wind shears typically exist in convective plume regions, Beres (2004) showed via comparison between linear and nonlinear models that the linear model gives reasonably good results with a mean background wind in the convective plume region.Also, Lane et al. (2003) found, using a numerical 3-D cloud-resolving model, that the excited GW spectrum in the intrinsic frame of reference is reasonably symmetric, thereby showing that the shear effect is of smaller importance than the overall Doppler effect of the moving wind frame at the tropopause.Therefore, we calculate the excited GW spectrum in the frame of reference of a constant mean wind in the region of the convective plume, then ray trace the GWs out of that region, as described in Sect.2.2.

Gravity wave solutions for vertical body forces
VF2001 solved the 3-D linear, Boussinesq, incompressible, f-plane equations with zero background winds for arbitrary spatial configurations of body forcings and heatings that are spatially and temporally localized.For these solutions, the temporal variation of the forcings and heatings are where the total duration is σ t , the forcing frequency is â=2π n/σ t , and the number of forcing cycles is n.These solutions are calculated in spectral (i.e., Fourier) space; for example, VF2001 solved for the spectral horizontal velocity perturbation u, where the GW zonal perturbation velocity is written as the tilde denotes a Fourier transform, x, y, and z are the geographic zonal, meridional and vertical coordinates, respectively, and k, l, and m are the zonal, meridional, and vertical wavenumbers, respectively.The solutions during and after the zonal (F x (x)), meridional (F y (x)), and vertical (F z (x)) body forcings and heatings (J (x)) are given by Eqs.(3.13-3.17) in VF2001.Because the equations are linear, the solutions in spectral space after the forcings/heatings are finished consist of a mean portion (constant in time) plus a GW portion (oscillates in time).
Here, we use these solutions to model a deep, convective plume.The air in a deep, convective plume moves upwards and "overshoots" the tropopause by 1-3 km into stably-stratified air; thereafter, this air collapses downwards onto the anvil and ceases strong vertical motions (Lane et al., 2001(Lane et al., , 2003)).Since a vertical body force produces a strong, localized updraft and downdraft of air, we model a single convective plume as a vertical body force(s) with a single oscillatory cycle (i.e., n=1).Substituting in F x =F y =J =0, the spectral space solutions after the vertical body force is finished (i.e., when t≥σ t ) from Eqs. (3.13-3.17) in VF2001 are where d=1/[σ t ω 2 ( â2 −ω 2 )] and Here, N is the buoyancy frequency, u, v, and w are the zonal, meridional, and vertical velocities, , ω is the intrinsic frequency, f =2 sin ξ , is the Earth's rotation rate, ξ is the latitude, and overbars and primes denote mean and perturbation quantities, respectively.Here, the dissipativeless, Boussinesq GW dispersion relation is Table 1 lists some useful symbols used in this paper.Note that there is no mean response to a vertical body force.Since we are modeling the excitation of GWs from convective overshoot, which produces very high-frequency GWs with 5-20 min periods (e.g., Larsen et al., 1982;Alexander et al., 1995), we can neglect the Earth's rotation by setting f =0.
The solution after the vertical body force is finished, Eqs.(3-7), then become Here, we have added the subscripts "GW" to emphasize that these solutions are the GW perturbation solutions.For horizontal body forcings and heatings, a mean response is also generated.The FL real-space solutions in an isothermal, zero-wind atmosphere is determined by taking the inverse Fourier transform of Eqs.(10)(11)(12)(13)(14).Note that the GW response depends sensitively on the temporal variability of the force; if the forcing is very slow in time (e.g., σ t →∞), there is no GW response.Since deep convective plumes have large vertical extents (e.g., D z ∼(1/2−2)H), where H is the density scale height, if we approximate these plumes to move upwards in time much faster then the buoyancy period 2π/N 5 min (e.g., an impulsive forcing), then the FL solutions do not accurately represent the excited GW spectrum because of the inaccuracy of the Boussinesq dispersion relation for high frequencies (see Sect. 5).
For this model, we set N x =N y =128 and N z =256 grid points in the x, y, and z directions, respectively.Additionally, the background temperature is T =240 K, N=0.02 rad/s, and f =0.We allow the equally-spaced x, y, and z grid spacings, x, y, and z, respectively, to vary, in order to optimize the reconstructed wave field from ray tracing (see Sect. 2.2).Defining the x, y, and z domain lengths to be L x , L y , and L z , respectively, the grid spacings are Here, we choose y= x, and z= x/2.horizontal velocity perturbation along GW direction of propagation u H 0 3-D FFT of u H (x, y, z, t) after multiplying by exp(−z/2H) w 0 ,... 3-D FFT of w(x, y, z, t) after multiplying by exp(−z/2H)...
Fourier Transform (FFT) used to calculate the spectral solutions are then respectively.Note that smaller grid spacings in real space result in larger grid spacings in spectral space.The largest k, l, and m values are kN x /2, lN y /2, and mN z /2, respectively.Therefore, the minimum zonal, meridional, and vertical wavelengths which are represented in an excited GW spectrum are respectively.

Gravity wave spectral amplitudes and phases
To determine the spectral momentum flux amplitude of a GW with wavenumber (k, l, m) needed for ray tracing, we average the zonal and meridional momentum fluxes, u w * and v w * , respectively, over a wave period, τ =2π/ω.Here, the * denotes the complex conjugate.As is clear from Eqs. (10-14), these spectral amplitudes do not change in time apart from the oscillations of C and S after the force is finished; therefore, we compute these amplitudes when the force finishes, i.e. at t=σ t .Because C and S depend only linearly on sin ωt and cos ωt, we utilize the following expressions: Thus each integral can be calculated simply by evaluating the integrand at t =t and at t =t+τ/4, adding them together, and dividing by two.Therefore, we calculate the time-averaged momentum flux amplitude of a GW as: We multiply by two to obtain the maximum momentum fluxes.Then the maximum dimensionalized zonal and meridional fluxes of vertical momentum at t=σ t are respectively.Here, k l m is included for proper dimensionalization of the spectral amplitudes, because Parseval's theorem states that for any function u, Note however that Parsevals theorem cannot be strictly applied to Eqs. (23-24) because of cross terms of u, v and w.
Eqs. (23)(24) are the GW momentum flux amplitudes we utilize in our ray trace code.
We also wish to calculate the phases of the excited GWs.It seems most sensible to calculate these phases at t=σ t /2 when the body force amplitude is maximum; however there are additional forced oscillations at the frequency â which cause the GW phases to be ambiguous at that time.Because the solution at t≥σ t is composed only of GWs, we calculate each GW's phase at t=σ t , which is the same time we calculated its amplitude (see .Because the GW's amplitude in spectral space is complex, we represent the GW's initial vertical perturbation velocity at t=σ t , for example, as where a is the maximum wave amplitude, φ FL is the FL Boussinesq wave phase, and the subscript "FL" stands for "Fourier-Laplace", as before.Both a and φ FL are real.Note that the initial wave phase from Eq. ( 2) is: Using Eq. ( 26), we then compute the phase of the vertical velocity at t=σ t as: where "real" and "imag" are the real and imaginary parts, respectively, of the complex number.The phase φ FL lies between 0 and 2π.Note from Eqs. (10-12) that u GW , v GW , and w GW are offset by either 0 or 180 o in phase, since they all oscillate as ±S.

Convective plume geometries and ground reflection
We model a single convective plume centered at (x, y, z)=(x 0 , y 0 , z 0 ) as a vertical body force(s) with where g(z) is the vertical distribution of the body force.Additionally, because the total vertical body force is F z F, where F has units of s −1 (see Eq. 1), F 0 is the "force amplitude" with units of m s −1 .In previous papers (e.g.VF2004, VF2006), g(z) was assumed to be a Gaussian in z (see Eq. 31 below).However, because the ground was not included, the air's vertical velocity was not zero at the Earth's surface, nor did the excited downward-propagating GWs reflect upwards from the Earth's surface.
The assumption that the Earth's surface is unimportant for the excitation and propagation of GWs from convective plumes, however, is unrealistic.Bristow et al. (1996a) observed GWs excited from the aurora with horizontal wavelengths of λ H 200−450 km propagating more than 1000 km from their source region which was 2000 km away.Many of these GWs are believed to be Earth-reflected (Samson et al., 1989(Samson et al., , 1990;;Bristow et al., 1994Bristow et al., , 1996b)).Additionally, several of the medium-scale GWs observed during the SpreadFEx experiment were reverse ray traced to convective sources positioned at the reflected source locations (Vadas et al., 2009a).Because the solutions, Eqs.(10-14), are Fourier decompositions in all space, including reflection off the ground is equivalent to inserting a boundary condition at z=0 such that the vertical velocity equals zero there, since air cannot penetrate below the Earth's surface.
In order to include ground reflection, we model a single convective plume as an upward-moving vertical body force centered at (x, y, z)=(x 0 , y 0 , z 0 ) plus a downward-moving "image" vertical body force below the Earth's surface at (x, y, z)=(x 0 , y 0 , −z 0 ).This image body force pushes air downwards at the same time as the body force above the ground pushes air upwards.Additionally, this image body force has the same horizontal and vertical dimensions as the body force above the ground.This causes the GWs excited by the above-ground body force to cancel those from the image body force at z=0, thus causing the vertical velocity to automatically be zero at the ground.Additionally, because the downward-propagating GWs from the aboveground body force reach z=0 at the same time that the upward-propagating GWs from the image body force reach z=0, it appears that those downward-propagating GWs "reflect" upwards at the Earth's surface.Thus the reflected GWs are actually those upward-propagating GWs from the image body force.These vertical distributions are represented mathematically as Here, Eq. ( 31) does not include GW reflection off the Earth's surface, while Eq. ( 32) does include GW reflection off the Earth's surface.We note that the configuration of body forces which results in the reflected spectrum, given by Eq. ( 32), does not take into account wind shear, topography, or boundary layer effects.The full zonal, meridional and vertical widths and depth of each body force are D x ≡4.5σ x , D y ≡4.5σ y , and D z ≡4.5σ z , respectively.We choose equal zonal and meridional widths, i.e.D≡D x =D y .For all of the convective plumes modeled in this paper, we set D=20 km and D z =10 km.Additionally, all of the convective plumes in Sects.2-4 have a duration of σ t =12 min.
We also model a small cluster of convective plumes which appear close together and are nearly simultaneous on satellite images (see Vadas et al., 2009a).If body forces are nearly simultaneous and are separated by less than a few diameters horizontally, the excited GW spectrum is not simply the sum of the excited GW spectrum from each body force individually, because GWs with larger horizontal scales are excited in proportion to the increased size of the cluster (Vadas et al., 2003).We define a small convective cluster to consist of 3 convective plumes which move the air in the troposphere upwards/downwards simultaneously.We situate each of these convective plumes at the corner of an equilateral triangle with sides of length D. If the center of the cluster is (x 0 , y 0 ), the locations of the 3 convective plume are A given "force amplitude" F 0 in Eq. ( 30) results in a maximum vertical velocity of a convective plume/cluster of w max .
If we wish to model a convective plume/cluster with a maximum updraft velocity of w pl instead, then we simply multiply the calculated wave amplitudes by w pl /w max because the solutions are linear.Additionally, we multiply the wave amplitudes by 1/2 because the amplitudes may be overestimated by a factor of ∼2 in linear models as compared to non linear models of convection (Song et al., 2003).Therefore, we multiply the GW amplitudes by w pl /2w max prior to ray tracing, e.g.,

Gravity wave excitation from convective plumes
Figure 1a shows the vertical body force used to model the upward surge of air within a single convective plume with no ground reflection, Eqs.(30)(31).This body force evolves as sin 2 (π t/σ t ) in time from Eq. ( 1); it begins pushing air upwards at t=0, reaches its maximum force at t=σ t /2, and ends at t=σ t .Since we do not include the Earth's surface, this force is the same as in previous plume models (VF2004; VF2006).The GW zonal and vertical velocity perturbations are shown at t=20 min in Fig. 1b and c, respectively.These GWs, excited by the rapid, upward movement of air, are clearly visible propagating symmetrically away from the center of the body force at x=y=0 and z=7 km.Because reflection is not included here, downward propagating GWs continue to propagate downwards below the Earth's surface.
Figure 1d shows the vertical body force we use to model the same convective plume, but with ground reflection via Eqs.( 30) and ( 32).The resulting GW zonal and vertical velocity perturbations are shown at t=20 min in Fig. 1e and f, respectively.It is easy to see the upward movement of air associated with the body force centered at z=7 km and the downward movement of air associated with the image body force centered at z=−7 km in Fig. 1f.GWs are clearly seen propagating away from both body forces at x=y=0, and at z=7 km for the force above the Earth's surface, and at z=−7 km for the force below the Earth's surface.Figure 1f also shows that the vertical velocity perturbation is zero at z=0.With this body force configuration, GW reflection off the Earth's surface is clearly visible; because those upward-propagating GWs excited by the image body force at z<0 reach the ground at the same time as the downwardpropagating GWs from the body force at z>0, it appears as if these latter downward-propagating GWs reflect upwards at the Earth's surface.Because of ground reflection, all of the GWs above the Earth's surface trace back to either z=7 km or to z=−7 km.This is particularly noticeable in Fig. 1e.
Although the non-reflected and reflected GW fields appear similar in Fig. 1, the excited GW spectra exhibit clear differences.Figure 2 shows 2-D spectra of the vertical flux of zonal momentum for GWs with l=0 in the functional form u w * multiplied by k and m: where the overline denotes an average over a wave period.Note that the factor of 2 converts an average to a maximum value.We use this functional form because a GW on any contour in Fig. 2 has the same value of the momentum flux (per unit mass), uw, in real space at any given altitude z, barring temperature and wind effects.(We show that this is the case in Sect.4).Note that although GWs along the same contour line have the same momentum flux amplitudes at a given altitude, they reach that altitude at different times, because they have different vertical group velocities (pink dash-dot lines in Fig. 2).
In Fig. 2, the upper row shows the GW momentum flux spectrum when there is no ground reflection, while the lower row shows the GW spectrum when reflection from the Earth's surface is included.In Fig. 2a and c, Eq. ( 35) is plotted as a function of k and m, while in Fig. 2b and   d, Eq. ( 35) is plotted as a function of λ x and λ z .The GW spectrum shown in the upper row of Fig. 2 is similar to the spectrum shown in Fig. 7 from VF2004, with a single peak at λ H ∼40 km and λ z ∼20 km.The GW spectrum shown in the lower row of Fig. 2, however, has two peaks at λ z ∼10 km and λ H ∼25−30 km and at λ z ∼25 km and λ H ∼40−60 km.This double peak occurs because there is partial destructive interference of the waves at vertical scales of λ z ∼12−18 km.This vertical scale of ∼14 km is the approximate distance between the forces.The larger GW amplitudes at larger λ z occur because there is a cohesiveness of the combined forcings over a larger vertical depth.We note that a double peak was also seen in the intrinsic k and ω I r spectrum for GWs excited by convection using a small-scale, highly-resolved convection model (Fig. 20 of Lane et al., 2003).
The maximum of the reflected spectrum is ∼3.5 times larger than that from the unreflected spectrum.This is partly because there are twice as many GWs in the reflected as in the unreflected spectrum, and partly because the reflected spectrum contains more power in larger vertical scales.Because the grid spacing is x=6.7 km, the minimum horizontal wavelength is 13 km from Eq. ( 17).This cutoff is apparent in the vertical group velocity contours.However, we note that the GW amplitudes are less than 10% for our plume model at λ H ≤13 km.This is because our plume model is designed to calculate the larger-scale GW amplitudes, rather than the smaller-scale GW amplitudes with λ H ∼5−10 km that dominate the convective spectra because of small-scale updrafts (Lane and Sharman, 2006).These smaller-scale GWs, however, are unimportant in the MLT and TI due to wave breaking, critical level absorption, and reflection in the stratosphere (e.g., Lane et al., 2003).
Figure 3a and d shows vertical slices of the body forces employed for the single convective plume and the small convective cluster, respectively.Both include ground reflection.The image vertical body forces can be clearly seen.Figure 3b and e shows horizontal slices of these forces.We note that the plumes are close together for the small convective cluster.Figure 3c and 3f shows the zonal wind velocity perturbations at t=20 min.GWs can clearly be seen propagating away from the center of the vertical body forces.The upward reflection at the ground of those initially downwardpropagating GWs can also be seen.
Estimating w * ∼− u * k/m from the Boussinesq continuity equation with l=0, the square root of Eq. ( 35) yields the horizontal velocity spectrum: (36) Figure 4 shows the GW momentum flux spectrum, Eq. ( 35), in the left column, and the zonal velocity spectrum, Eq. ( 36), in the right column, when ground reflection is included.The upper (lower) row shows the spectra for a single convective plume (small convective cluster).For this figure, we calculate u u * using the anelastic expressions, Eqs. ( 57) and ( 55) with l=0 and ν=α=0.Note that the zonal velocity spectra are broader in λ H and λ z than the momentum flux spectra because of the square root.We see that the small convective cluster excites GWs with larger horizontal and vertical scales than the single convective plume, because the forcing is cohesive over a larger horizontal extent.
We now calculate the real-space momentum fluxes of the GWs excited from a single convective plume in a windless, isothermal atmosphere using the FL solutions; we do this by taking the inverse Fourier transform of Eqs.(10)(11)(12)(13)(14).This results in u(x, y, z, t) and w(x, y, z, t), for example.Since the FL solutions assume that the background density is constant with altitude, the resulting GW amplitudes do not grow with altitude as they should.These solutions do, however, account for wave dispersion.GW amplitudes, such as u, w, and ρ /ρ, grow as 1/ √ ρ (Hines, 1960).In an isothermal atmosphere, T is constant and ρ decreases exponentially with altitude, ρ=ρ( 0 air is unstable below the tropopause during active convection, and therefore cannot support the excitation and propagation of GWs.Therefore, we multiply the FL solutions by exp((z−z trop )/2H) to account for the growth of a GW's amplitude with altitude above the tropopause, where z trop is the altitude of the tropopause.Additionally, we multiply a GW's amplitude by w pl /2w max from Eq. ( 34).Therefore, we multiply the FL solutions u, v, w, , and P by Here, we set z trop =15 km and w pl =40 m/s.Note that we mistakenly left out the factor of 2 in calculating the FL solutions in VF2006.
Figure 5 shows the unaveraged GW zonal momentum fluxes, uw, created from a single convection plume with ground reflection at z=50 km (left column), z=70 km (middle column), and z=90 km (right column).The upper to lower rows show these winds at t=25, 35, 45, and 55 min, respectively.The vertical body forces which model this plume begin at t=0.We see that the GWs propagate upwards from the plume as concentric rings.Because the momentum flux is proportional to sin 2 (k H r), where r= (x−x 0 ) 2 +(y−y 0 ) 2 is the radius from the center of the force, the horizontal wavelength is approximately twice the distance between two white contours.At z=90 km and r∼100 km, the GW horizontal wavelengths are seen to decrease with time, from λ H ∼80 km at t=25 min to λ H ∼40 km at t=55 min.We note that these FL solutions are not accurate at times greater than 60 min for  the numerical box we use here, because some of the fast GWs reflect off the box walls, causing constructive and destructive interference with the other GWs in the box (VF2006).Some of this interference is noticeable at t=55 min.Since these solutions assume constant temperatures and zero winds, they are of limited utility; however, they are of great importance for normalizing the GW spectral amplitudes in ray tracing, as we will see in Sect.3.
Figure 6 shows the corresponding GW zonal, meridional and vertical velocity perturbations in the left, middle, and right columns, respectively, at z=70 km for the same convective plume as in Fig. 5.The upper and lower rows show t=25, 35, 45, and 55 min.As before, the concentric GW rings are clearly visible.Since the perturbation velocities are proportional to sin(k H r), the horizontal wavelength of the GW is approximately the distance between two white or two black contours.Note that the zonal velocity is asymmetric about x=x 0 and the meridional velocity is asymmetric about y=y 0 .The vertical velocity, however, is symmetric in x and y, as is the temperature perturbation (not shown).Note that concentric rings from convective storms, like those in Fig. 6, have been occasionally observed in the OH airglow layer from deep convection (Taylor and Hapgood, 1988;Dewan et al., 1998;Sentman et al., 2003;Suzuki et al., 2007;Yue et al., 2009).

Methodology
Ray tracing has been used for decades for geophysical problems of interest (e.g.,Jones, 1969;Marks and Eckermann, 1995;Cowling et al., 1971;Waldock andJones, 1984, 1987;Hung and Kuo, 1978;Hung and Smith, 1978;Lighthill, 1978;Gerrard et al., 2004;Hecht et al.2004;Lin and Zhang, 2008).Our ray trace model is based on the formalism developed by Lighthill (1978), and allows the wind, density, and other background parameters to change slowly with altitude, horizontal location, and time.If a wave packet is propagating in a background wind, V(x)=(V 1 , V 2 , V 3 )=(U, V , W ), then its evolution in space and time is described by the following equations: ) where the components of the vector group velocity, c g , are c g i =∂ω I r /∂k i , ω I r is the real part of the intrinsic frequency of the GW: and ω r is the real part of the ground-based frequency.Here, the indices i, j =1, 2, 3 indicate the components of the vector quantities x, V , k, and c g , and repeated indices imply a summation.The subscript "r" here denotes the real component of the frequency (Vadas andFritts, 2005, hereafter VF2005), not the relative frequency; the relative frequency (i.e., the frequency in the intrinsic frame of reference moving with the fluid) is denoted by the subscript "I ", where "I " stands for "intrinsic" (see Table 1 for definitions).We do not allow the ground-based frequency of a GW to vary in time in our model.Fourier analysis has shown that energy transfer cannot occur between waves of different frequencies for a set of linear equations when the coefficients do not depend explicitly on time (Lighthill, 1978).For all of the simulations in this paper, in Fritts and Vadas (2008), and in Vadas et al. (2009b), the background parameters are constant in time; therefore, ω r is constant in time for each GW.However, for the reverse and forward ray trace simulations performed for this SpreadFEx campaign in Vadas et al. (2009a), we assume that the background temperatures and winds vary slowlyenough to approximate that ω r is constant in time.We note that several formulations do take into account changing ω r with time as background conditions change with time (e.g.,Jones, 1969).
The GW perturbation quantities, such as u and w, are assumed to be approximately sinusoidal, varying as where is the wave phase as defined in Eq. ( 18) of VF2005.Note that the initial phase for this ray trace model is negative that of the initial phase from the FL convective plume model, φ(t=0)=−φ FL (compare Eqs. 27 and 42 at t=0).From Eq. (42), The change in a GW's phase in time as calculated along the ray path (i.e., at the GW's location as represented by a point particle propagating at speed dx i /dt) is Here, we have used Eqs.( 38), (40), and (43) and the fact that k i c g i =k•c g =0 for a GW, since a GW's phase velocity is perpendicular to its group velocity (e.g., Kundu, 1990).Therefore, a GW's phase along its ray path is

Dissipative, anelastic GW dispersion relation
The GW dispersion relation we use here includes the primary damping mechanisms for high-frequency GWs with large vertical wavelengths: kinematic viscosity and thermal diffusivity.It is nonhydrostatic and compressible, but excludes acoustic waves, similar to Marks and Eckermann (1995).It also neglects the Earth's rotation.Because the Prandtl number is Pr ≈0.7 throughout the mesosphere and lower thermosphere (e.g., Kundu, 1990), kinematic viscosity and thermal diffusivity have comparable damping effects on a GW at approximately the same altitude.This anelastic GW dispersion relation can be rewritten, using Eq. ( 26) from VF2005, as where ν=µ/ρ is the kinematic viscosity, µ is the viscosity coefficient, δ=νm/Hω I r , and δ + =δ(1+ Pr −1 ).This dissipative dispersion relation yields the usual high-frequency GW anelastic dispersion relation when dissipation is negligible, obtained by setting ν=δ=δ + =0 in Eq. ( 47): (48) (Gossard and Hooke, 1975;Marks and Eckermann, 1995).Equation ( 47) differs from the anelastic dispersion relation in Marks and Eckermann (1995) (i.e., Eq.48 when f =0) in that thermospheric dissipation is taken into account for high-frequency GWs which propagate above the turbopause (at z∼110 km).Eq. ( 47) reduces to the Boussinesq dispersion relation (Eq.9) with f =0, when H→∞.Note that the dispersion relation we use here neglects other forms of dissipation such as ion drag and wave-induced diffusion.Additionally, δ is negative for an upward-propagating GW, since m<0.
The inverse decay rate in time for a dissipating GW (Eq. 25 from VF2005) is Since ω I i changes in space and time, we integrate in time along ray paths.Therefore, including GW phases, a GW's horizontal spectral momentum flux (per unit mass) when launched from z=z i and t=t i , divided by k l m, is where we put an absolute value around ω I i to ensure that a GW dissipates even if k 2 <1/4H 2 .Here, u H is the horizontal velocity perturbation along the direction of propagation of the GW.Because Eq. ( 50) does not include the decrease of a GW's amplitude as it propagates into a larger area, we smooth the momentum fluxes horizontally prior to reconstructing the GW field.An important assumption used to derive the expressions for the GW dissipative dispersion relation and GW amplitude decay in time is that the background wind shears are not too large (VF2006).If U H is the background wind in the direction of GW propagation, then In this paper, U H =0, so Eq. ( 51) is automatically satisfied.

Reconstruction of the gravity wave field
In this subsection, we describe how we reconstruct the GW zonal, meridional, and vertical velocity fields, and the temperature and density perturbation fields, using the binned and horizontally-smoothed momentum fluxes from ray tracing.Crucial here is the inclusion of wave phases.We note that mountain waves have been reconstructed to high accuracy using Fourier-ray methods.These methods involve ray tracing in the spectral domain, then inverse Fourier-transforming to the spatial domain (Broutman et al., 2004(Broutman et al., , 2006)).Our method, in contrast, involves calculating the wave amplitudes in the spectral domain, ray tracing the spectral amplitudes and phases in the spatial domain, smoothing horizontally to account for wave dispersion, then multiplying by a normalization factor to convert the spectral amplitudes to real-space amplitudes.
The anelastic, GW, polarization relation between the wind along the direction of GW propagation and the vertical wind is (Eqs.B4-B5 from VF2005): where α≡−k 2 +1/4H 2 +im/H, u H 0 is the spectral scaled horizontal velocity perturbation, w 0 is the spectral scaled horizontal velocity perturbation, ω I =ω I r +iω I i is the complex intrinsic frequency, and Here, the scaled perturbation quantities have subscripts "0", and denote taking the Fourier transforms of the GW perturbation quantities after they have been multiplied by exp(−z/2H ) (see Eq. 17 from VF2005).The complex anelastic dispersion relation is given by Eq. ( 23) from VF2005:  54) into Eq.( 52), we get Equation ( 55) reduces to the correct equation in the nondissipative, Boussinesq limit; setting H →∞ and ν=0, Eq. ( 55) becomes w 0 −(k H /m) u H 0 , which is simply the Fourier transform of the 2-D Boussinesq continuity equation (e.g., Kundu, 1990).When dissipation is unimportant and 4π λ z H, w 0 and u H 0 are in phase (but may have opposite signs) for an upward-propagating GW, because m<0 but k>0 (k<0) for eastward (westward) propagating GWs.However, since Eq. ( 55) is complex in general, the horizontal and vertical velocity components may not be in phase for very large λ z and during strong dissipation.The vertical velocity squared is determined via multiplying Eq. ( 55) by the spectral momentum flux: The zonal velocity squared is determined via: The zonal and meridional components are then determined from Eq. ( 57) via The temperature and density perturbations are then determined from the GW dissipative, anelastic, polarization relations (Eqs.B6-B7 from VF2005): respectively.

Inputs, outputs, and specifications
The convective plume model outputs the average amplitudes, scales, intrinsic frequencies, and phases of the GWs excited from single or multiple convective plumes at t=σ t .The ray trace model then inputs these quantities.However, this inputted intrinsic frequency was calculated from the Boussinesq dispersion relation, Eq. ( 9), with f =0.We therefore assume that the amplitude calculated from the Boussinesq model is correct for a given GW with wavenumber vector (k, l, m), and recalculate the GW's intrinsic frequency using Eq. ( 48) and the ground-based frequency using Eq. ( 40).By doing this, we ensure that ω I r is compatible with the dispersion relation in the ray trace model.Note, however, that the anelastic dispersion relation does not allow GWs to have frequencies near N, as is possible with the Boussinesq dispersion relation.This is because the cutoff frequency limits the highest frequency obtainable; this is a large effect for GWs with large λ H (Marks and Eckermann, 1995).
After recalculating ω I r , we ray trace each GW through varying winds and temperatures, and add its horizontal momentum flux amplitude from Eq. ( 50), wavenumbers, intrinsic frequencies, and vertical velocity phases into large 4-D arrays.Here, the wavenumbers, intrinsic frequencies, and phases of each GW are weighted at (x, t) by the GW's total horizontal momentum flux.Note that we only need to keep track of the phase of the vertical velocity, because all of the other quantities are related to w 0 through the polarization relations (see Eqs. 52, 59, and 60).These 4-D arrays are binned in (x, y, z, t), with bin sizes of 4 km, 4 km, 2 km, and 2 min in x, y, z, and t, respectively, for x=[−240, 240] km, y=[−240, 240] km, z=[20, 100] km, and t=[0, 2] h.After all GWs have been ray traced, each bin contains the total momentum flux, average horizontal and vertical wavenumbers, average intrinsic frequency, and average phase of the vertical velocity perturbation.
Because we are comparing our ray trace solutions with the FL solutions in this paper, we chose an isothermal temperature profile, zero winds, and ν=0.Thus, the GWs we ray trace here experience no critical level filtering, no dissipation, and no thermal or Doppler ducting.Other ray trace studies utilize realistic, variable winds and temperatures (Vadas et al., 2009a, b).Here, we set T =240 K, N=0.02 rad/s, H=7 km, µ=0, Pr =0.7, γ =1.4 and X MW =29.We ray trace each GW from x=y=0 and z=z 0 =7 km; these were previously shown to be good assumptions (see Fig. 1e-f).We also ray trace each GW from the time the vertical body force is maximum, at t=σ t /2=6 min, even though the GW phases and amplitudes are calculated at t=σ t (when the force is finished).We do this because the excitation and radiation of GWs is strongest at t=σ t /2; thus t=σ t /2 best represents when convective overshoot occurs.From a technical point of view, if we instead ray trace the GWs from t=σ t , the ray trace solution is exactly the same, but is shifted forwards in time by σ t /2=6 min.We will show in the next section that ray tracing the GWs from the middle of the body force at the peak time of t=σ t /2 agrees very well with the FL solutions in space and time for σ t =12 min.Finally, to achieve accurate numerical solutions, a 4th-order Runge Kutta routine (Press et al., 1992)

Reconstruction of the wave field in real space from ray tracing
In this section, we ray trace the excited GW spectrum from a single convective plume through an isothermal, windless atmosphere up through the OH airglow layer.We calculate the resulting momentum flux and velocity perturbations as a function of altitude, radius, and time, up to a constant normalization factor.We then compare these (anelastic) solutions with the FL (Boussinesq) solutions at the same times and altitudes.This allows us to determine the normalization factor in the limit that the two solutions are equal, which is very important because it allows for the conversion of the ray traced spectral wave amplitudes to real-space wave amplitudes for generalized wind and temperature profiles (e.g., Vadas et al., 2009a, b).

Grid space, spectral resolution, and concentric rings
First, we demonstrate how changing the grid spacing of the convective plume model affects the distribution of wavelengths in the excited GW spectrum.In Fig. 7a-d, we show the reflected GW horizontal velocity spectra (solid lines) as a function of λ H and λ z for a single convective plume with ground reflection and with horizontal resolutions of x=4.4,6.7, 8.9, and 13.3 km, respectively.We also overplot each GW in the excited spectrum by a dot.Although the grid spacings increase in real space from Fig. 7a to d, they decrease in spectral space from 7a to 7d using Eq. ( 16).As x increases, the spectrum is increasingly represented by GWs with larger horizontal wavelengths.The spectrum in Fig. 7a is focused on smaller scales with λ H ∼2−50 km, while the spectrum in Fig. 7d is focused on larger scales with λ H ∼25−100 km.Most importantly, there are far more GWs representing the largest horizontal scales of λ H =100−400 km in Fig. 7d as compared to Fig. 7a.Note that the wave amplitudes for the GWs with very large λ H and very large λ z in Fig. 7a are much larger than in Fig. 7d, because each of these GWs in Fig. 7a carries a larger amount of momentum flux in order to represent the same excited spectrum of GWs.The minimum horizontal wavelengths for the spectra shown in Fig. 7a, b, c, and d are λ H ≥8.9 km, λ H ≥13.3 km, λ H ≥17.8 km, and λ H ≥26.6 km, respectively, from Eq. ( 17).
We show these minima wavelengths in Fig. 2a and c as vertical dotted lines.Although the spectra with minima horizontal wavelengths smaller than 17.8 km (or a grid resolution of x≤9 km) has a negligible effect on the excited GW spectrum, the spectrum with a minima horizontal wavelengths of 26.6 km (or a grid resolution of x=13.3 km) eliminates some of the excited GWs.(Remember, however, that our model does not include small-scale updrafts, as discussed previously.)Clearly, the resolution one chooses should be based on the GW scales important for the problem at hand.If one is ray tracing GWs into the mesosphere and thermosphere where scales of λ H ∼30−400 km are important, then one might choose x=8−13 km, whereas if one is ray tracing GWs into the stratosphere where small-scale waves are important instead (e.g., Lane and Sharman, 2006), one might choose x=1−6 km instead, and include smaller-scale dynamics.
Figure 8 shows the zonal flux of vertical momentum of the GWs in the airglow layer at z=90 km at t=45 min for the same convective plume as in Fig. 5, but obtained instead from ray tracing the GW spectra.In Fig. 8a-d, the GW spectra used are shown in Fig. 7a-d, respectively.Since these are the reflected spectra, we only ray trace those GWs which propagate upwards from z=7 km.Note that the GW amplitudes are scaled to reflect this convective plume's updraft velocity of w pl =40 m/s from Eq. ( 34) prior to ray tracing, and that they are additionally multiplied by where the tropopause altitude is z trop =15 km, so that the wave amplitudes are correct at the tropopause (i.e., not increasing with altitude below the tropopause).We see that circular, concentric rings result from ray tracing the GWs, especially in Fig. 8b-d, and that these rings are located at similar radii in Fig. 8b-d.Comparing with Fig. 5i, we see that the circular rings in Fig. 8b-d agree well qualitatively with the FL solutions.We will perform a more detailed comparison of the FL and ray trace solutions momentarily.Because it is difficult to distinguish one concentric ring in Fig. 8a from another, it is clear that the spectrum in Fig. 7a is completely inadequate for reconstructing the wave solution in the OH airglow layer and at higher altitudes.This is because this spectrum does not have adequate resolution for GWs with large vertical wavelengths of λ z >30 km and with intermediate frequencies of ω I r /N<0.2−0.6.The spectrum used in for Fig. 8b is better, although there are quite a few locations at lower wave frequencies (i.e., larger radii) where not enough GWs are present to represent the wave field adequately (e.g., at r∼140−220 km).The concentric rings in Fig. 8c-d are well resolved.However, because Fig. 8c (8d) includes GWs with λ H ≥17 km (λ H ≥27 km), the spectrum in Fig. 7c is likely the best choice for reconstructing the wave field at z∼90 km with a GW spectrum represented by N x =N y =128 and N z =256 grid points.If a larger number of grid points were employed instead, a different spectral grid spacing might be optimal.Note that for ray tracing convectively-generated GWs into the thermosphere, the spectrum in Fig. 7d may instead be the best choice.There-fore, we only consider GW spectra excited from convective plumes with 8.9≤ x≤13.3 km (and 4.4≤ z≤6.7 km) here.

Reconstruction and normalization of the wave field in real space
There are 2 possible methods by which to ray trace the GWs excited from our modeled convective plume.First, we can ray trace only those upward-propagating GWs using the amplitudes from the reflected spectrum (Fig. 2c-d) from the center of the body force at z=7 km.This is the approach we used in Fig. 8. Second, we can ray trace those upward and downward-propagating GWs using the amplitudes from the unreflected spectrum (Fig. 2a-b) from the center of the body force at z=7 km, and allow the downward-propagating GWs to reflect upwards from the ground.We employ the former method first.
Figure 9 shows the zonal flux of vertical momentum in real space from ray tracing the GWs excited from the same convective plume as in Fig. 5.These solutions are shown at the same altitudes and times as in Fig. 5. Here, we use the GW spectrum which includes reflection off the Earth's surface, and therefore only ray trace those upward-propagating GWs from z=7 km.Note that Fig. 9c is blank because no GWs have propagated to z=90 km at 25−6=19 min after convective overshoot.(Remember that the GWs are ray traced from t=6 min, at the maximum time of the body force, which we equate to be the physical time of convective overshoot).Comparing with Fig. 5, we see that the ray trace model wave scales (as a function of radius and time) agree well with those from the FL solution.
In Fig. 10a-c, we show the maximum values of the smoothed momentum fluxes from the FL solutions, multiplied by exp((z−z trop )/H), as diamonds at z=50, 70, and 90 km, respectively.We see that the momentum fluxes increase then decrease in time at each altitude.Because the density decreases with altitude, the momentum fluxes are much larger at z=90 km than at z=70 km.Here we use the reflected spectrum in Fig. 2c-d, and only ray trace those upward-propagating GWs.This plot is unsmoothed in x, y, z, and t.The resolution for the spectra is x=8.9 km.Fig. 10.Row 1: Average momentum fluxes from the smoothed ray traced solutions at z=50, 70, and 90 km in columns 1, 2, and 3, respectively.Solid (dash) lines show the ray trace results for the reflected (unreflected) GW spectrum normalized by Eq. (63) (Eq.62).Diamonds shows the average momentum fluxes from the smoothed FL solutions shown in Fig. 5 and at t=15 and 60 min.(Note that the FL solutions have been multiplied by exp((z−z trop )/H) to take into account the density decrease with altitude.)Row 2: Same as in (a), but multiplied by (z−z 0 ) 2 exp(−(z−z trop )/H).Row 3: Average GW horizontal wavelengths as a function of radius at z=90 km and y=0.
Results are shown at t=35, 45, and 55 min in columns 1, 2, and 3, respectively.Solid (dash) lines show the ray trace values using the reflected (unreflected) GW spectrum.Diamonds shows the results from the FL solutions.
Because the momentum fluxes from the ray trace model must equal the momentum fluxes from the FL model in this isothermal, windless limit, we divide the binned, spectral momentum fluxes, u w * k l m (which have units of momentum flux times a volume), by a volume factor ξ .For the unreflected GW spectrum, we employ ξ = x y z t/(50 s), and for the reflected GW spectrum, we employ The GW momentum fluxes (per unit mass) in real space are then computed as: where the sum is over all of the GWs that enter a (x, y, z, t) bin.The normalization factor, 1/ξ , is smaller for the reflected spectrum because all of the GWs arrive at a given altitude at the same time; for the unreflected spectrum, half of the GWs arrive somewhat later because they first reflect off the Earth's surface.In Fig. 10a-c, we overplot the maximum values of the smoothed momentum fluxes from ray tracing using Eq. ( 64).The momentum fluxes using the reflected spectrum from Fig. 2c-d with Equation ( 63) is shown as solid lines, and the momentum fluxes using the unreflected spectrum from Fig. 2a-b with Eq. ( 62) is shown as dash lines.Here, we have smoothed the horizontal slices of the momentum fluxes as we did for the FL solutions, so that we are comparing average values in all cases.This smoothing consists of a 7 point running average (in both x and y), then a 3 point running average 3 additional times.
In the middle row of Fig. 10, we show the same momentum fluxes at the same altitudes as in the upper row of Fig. 10, but multiplied by to approximately account for the increase in GW amplitudes from the density decrease with altitude above the tropopause, and the decrease in GW amplitudes as the GWs disperse into increasingly larger areas for z>z 0 (Fritts and Vadas, 2008).Note that the curves in Fig. 10d-f have approximately the same value; therefore, the ray trace model (with horizontal smoothing) is properly taking into account wave dispersion and wave amplitude growth with altitude.We see that at all three altitudes, the momentum fluxes increase and decrease with time.At higher altitudes, the peak occurs at later times; this is expected, because it takes longer for the GWs to propagate to higher altitudes.The GWs from the reflected spectrum result in momentum fluxes which agree well with the FL solutions.The GWs from the unreflected spectrum, however, cause the momentum fluxes to decay too slowly with time after the peak time, as compared to the FL solutions.This may be due to incomplete (and incorrect) cancellation of wave amplitudes for GWs with λ H ∼12−18 km.Note that if we had instead ray traced both GW spectra from t=σ t rather than from t=σ t /2, each curve in Fig. 10a-f would have shifted by +6 min (later in time), which would have agreed less well with the FL results.Thus, ray tracing the reflected GW spectrum from t=σ t /2 (i.e., when the body force is maximum) is the optimal choice.Figure 10g-i shows the average horizontal wavelengths as a function of radius from the center of the convective plume at t=35, 45, and 55 min, respectively.We calculate the average horizontal and vertical wavelength in each bin in the ray trace model, λ H =2π/k H and λ z =2π/m, respectively, using Here, the sum is over all of the GWs which enter this bin.We see that the ray traced values agree well with the horizontal and vertical wavelengths of the FL solutions.Additionally, note that λ H increases rapidly with radius at a fixed time, but λ H decreases in time at a fixed radius.For example, at z∼200 km, λ H ∼110, 90, and 80 km at t=35, 45, and 55 min, respectively.
Figure 11 shows the reconstructed GW perturbation velocities in real space, u, v and w, at z=70 km and t=55 min.These wave fields are determined from the momentum fluxes using Eqs.(55-58), and by normalizing by Eq. ( 64).Note that we smoothed the binned momentum fluxes prior to reconstructing the wave fields in order to ensure that the wave amplitudes decrease as the waves spread out into a larger area from the source.In the upper row, we show the reconstructed fields using the unreflected spectrum from Fig. 2a-b with upward and downward-propagating GWs.In the lower row, we show the reconstructed fields using the reflected spectrum from Fig. 2c-d with upward-propagating GWs only.We see that the solutions are similar.However, because the reflected spectrum agrees better with the FL solutions in time (see Fig. 10c-f), we only use the reflected GW spectra and the normalization factor given by Eq. ( 63) to convert the ray traced GW spectral to real space amplitudes for the rest of this paper.
Figure 12 shows the reconstructed GW velocity fields in real space, u, v and w, at z=70 km for the same plume from Fig. 9.These solutions are shown at the same times as the FL solutions in Fig. 6.Although noisier than Fig. 6, especially at larger radii (i.e., lower frequencies) where there are proportionately less GWs represented in the spectrum (see Fig. 7), we see that the reconstructed velocity fields agree very well with the FL solutions.
We now compare the amplitudes and phases of the FL solutions with the ray trace solutions.In Fig. 13, we show the zonal velocity, vertical velocity, temperature, and density perturbations from left to right, respectively, at z=70 km and y=0.The upper to lower rows show the solutions at t=35, 45, and 55 min, respectively.The solid and dash lines show the ray trace and FL solutions, respectively.Overall, the agreement is excellent.In particular, the amplitudes agree very well.Additionally, the phases are in very good agreement for x<125 km for all times and for x<200 km for times greater than 50 min.There are several reasons that the phases do not agree as well for x>125 km for early times.First, the FL Boussinesq solutions assume that H→∞.Yet, the earliest GWs which arrive at this altitude have very large λ z >50 km for which the anelastic corrections (of order λ z /2πH) cannot be neglected.Second, the large-λ z GWs are not as well represented in the GW spectrum (see Fig. 7).Therefore, there may not be enough GWs in the spectrum to reconstruct the GW perturbation fields accurately at early times and large radii.At later times, when smaller λ H and λ z GWs reach this altitude, λ z /2πH can be neglected and the excited spectrum is better represented with more GWs (see Fig. 7).Upper row shows the results using the unreflected GW spectrum (Fig. 2a-b), while the lower row shows the results using the reflected GW spectrum (Fig. 2c-d).From left to right, the maximum values of u, v, and w are: 1st row: 17, 17 and 11 m s −1 .2nd row: 14, 14 and 9 m s −1 .The resolution is x=8.9 km.
The convective plume model, ray trace model, and reconstruction methods described in this paper have been recently applied to the modeling of an observed deep, convective plume near Fort Collins, Colorado, which produced nearly concentric rings in the OH airglow layer for ∼1.5 h (Yue et al., 2009).Using radar measurements and satellite imagery, we estimated that this plume had a width of D∼15 km, a depth of D z ∼10 km, a duration of σ t ∼10 min, and an updraft velocity of w pl =35 m s −1 .The GW spectrum excited from this plume was calculated using the formalism described in Sect.2.1, including ground reflection.The reflected spectrum was then ray traced through vertically-varying winds using the formalism described in Sect.3, with the calculated normalization factor given by Eq. ( 63).It was found that the concentric rings became "squashed" or arc-like if the intervening winds were strong, and that the apparent center of the concentric rings shifted with respect to the convective plume if the average horizontal wind between the tropopause and the OH layer was large (Vadas et al., 2009b).It was also found that the model results agreed reasonably well with the amplitudes, phases, horizontal wavelengths, and periods of the observed waves, and that a disappearance of part of the concentric GW rings after an hour could be explained by the reflection of high-frequency GWs propagating opposite to the direction of the background, horizontal wind.

Spectral amplitudes at airglow altitudes
In this section, we verify that the functional forms of the GW spectra used to plot Figs. 2 and 4, i.e.Eqs. ( 35) and (36), yield the same wave amplitudes at higher altitudes along any constant contour line, barring temperature and wind effects.We then calculate the normalized GW horizontal wind spectrum at the OH airglow altitude for a single convective plume and a small convective cluster.
We choose 2 GWs with values that are 10% and 20%, each, of the maximum value of the uw * spectrum in Fig. 4a.The parameters for these 4 GWs are shown in Table 2.We ray trace these individual GWs with l=0 from x=y=0, z=7 km, and t=6 min, and calculate their wavelengths as a function of altitude and time.We show λ x and λ z in Fig. 14a and b, respectively.The GWs from left to right in Table 2 are shown as solid, dot, dash, and dash-dot lines, respectively.As expected, λ x and λ z for these individual GWs are unchanged with altitude because the atmosphere is isothermal and windless.
Using these ray trace results, we calculate the locations and times that these individual GWs reach a given altitude.Using the ray trace results from Figs. 9 and 12, we calculate the nearest (x, y, z, t) bin to each of these locations and times for each GW.We overlay the average λ x and λ z values in  each of these bins determined from Eq. (66) in Fig. 14a and b, respectively.The GWs from left to right in Table 2 are shown as diamonds, triangles, squares, and x's, respectively.We see that the average wavelengths are reasonably accurate, yielding zonal and vertical wavelengths that are similar to the individual ray trace results.Differences are somewhat larger at lower altitudes because the source is closer and the bin sizes are too large to accommodate a single wave packet.Instead, a few wave packets likely enter each (x, y, z, t) bin, thereby skewing the average somewhat.However, Fig. 14a-b shows that the 4-D array is binned reasonably finely-enough for z≥50 km for our purposes here.
In Fig. 14c, we show the binned, smoothed momentum fluxes multiplied by (z−z 0 ) 2 exp(−(z−z trop )/H) at these same locations and times for the GWs from Fig. 14a-b.For all altitudes, we see that the momentum fluxes are approximately double for those bins closest to the GWs with 20% of the maximum spectral amplitude, as compared to those bins closest to the GWs with 10% of the maximum spectral amplitude.Variations occur partly because each bin contains several wave packets, as discussed previously, and partly because these GWs are in a region of the spectrum in Fig. 4a where the momentum flux changes rapidly with λ H and λ z ; therefore the binned momentum fluxes vary rapidly with radius.This latter problem might be lessened if we could choose large λ H and λ z GWs in the slowly-varying region of the GW spectrum; however, we cannot, because they exit our numerical box prior to reaching z=100 km.In conclusion,  we have shown that contours of Eq. ( 35) yield GWs with the same values of uw at higher altitudes in a windless, isothermal atmosphere.
We also choose 2 GWs with contours that are 30% and 60%, each, of the maximum value of the u H spectrum in Fig. 4b.The parameters for these 4 GWs are shown in Table 3.We ray trace these individual GWs with l=0 from x=y=0, z=7 km, and t=6 min.We show the results in Fig. 15 in the same format as in Fig. 14.As before, λ x and λ z are constant with altitude for these individual GWs (solid, dot, dash, and dash-dot lines for the GWs from left to right in Table 3, respectively).We see that the average zonal wavelengths (Fig. 15a) and average vertical wavelengths (Fig. 15b) for the bins nearest to the individual GWs (shown as diamonds, triangles, squares, and x's for the GWs from left to right in Table 3, respectively) are similar to the zonal and vertical wavelengths of the individual GWs.This implies, as before, that the bins are small enough to admit reasonably localized wave packets above z>50 km.In Fig. 15c, we plot the smoothed values of |u H | in the bins nearest to these four individual GWs.We see again that the zonal velocities are approximately double for those bins closest to the GWs with 60% of the maximum spectral amplitude, as compared to those bins closest to the GWs with 30% of the maximum spectral amplitude.Therefore, we have shown that contours of Eq. ( 36) yield GWs with the same values of u H at higher altitudes.
Figure 16a shows the GW zonal velocity spectra for the same convective plume as in Figs. 5 and 9 with ground reflection (solid lines).This plume has an updraft velocity of w pl =40 m/s. Figure 16b shows the GW zonal velocity spectra for a small convective cluster with the same plume parameters as in Fig. 16a.Both spectra were shown previously in Fig. 4b and d, respectively.Here, however, we show the values of the GW horizontal velocity amplitudes at z=87 km if the atmosphere is windless and isothermal, and if no wavebreaking occurs.Here, we plot the normalized zonal velocity spectra via combining Eqs. ( 36) and (64).Given λ H and λ z for a GW, these plots show the maximum zonal velocity that this GW has in the OH airglow layer.Because the zonal velocity oscillates in time and a wave packet amplitude increases and decreases in time, the amplitudes shown in Fig. 16 represent an upper bound on the zonal velocities.The maximum values of u H in Fig. 16 are 60 and 100 m s −1 for the single plume and the small convective cluster, respectively.Typically, however, the larger-amplitude GWs break at lower altitudes where their nondimensional amplitudes, u/(c−U ), reach one (Lindzen, 1981;Fritts and Alexander, 2003).This effect is not included in Fig. 16.
The GW spectra shown in Fig. 16 were used in Vadas et al. (2009a) to compare theoretical estimates with the observed amplitudes of the six medium-scale waves observed by Taylor et al. (2008) in the OH layer, in order to constrain this theory.We first converted the measured intensity amplitudes to horizontal velocity amplitudes.We then determined the vertical wavelength of each GW at the tropopause via reverse ray tracing through realistic winds and temperatures.Here, the value of λ z at the tropopause is important because the GWs are excited there via convective overshoot.For each of the GWs, we then plotted the observed λ H and λ z (at the tropopause) onto a spectra similar to Fig. 16, except for somewhat different convective plume parameters.We compared the modeled and observed amplitudes, and found that the measured wave amplitudes agreed reasonably well with theory in the cases where λ z >10 km and wave ducting did not likely occur.
The spectra shown in Fig. 16 were also used to infer GW amplitudes at higher altitudes in the TI for this SpreadFEx campaign (Fritts and Vadas, 2008;Fritts et al., 2008b).GWs with λ z > ∼ 100 km are those we expect to penetrate to the greatest altitudes, based on the viscous dispersion relation developed by VF2005 and the ray tracing studies by Vadas (2007) and Fritts and Vadas (2008).Fritts et al. (2008b) employed the results of this study to infer maximum horizontal velocities near z∼80 km of ∼1−2 ms −1 for GWs with λ z ∼150 km and λ H ∼200−400 km arising from a single convective plume, with amplitudes ∼2 times larger in response to a small convective cluster.Assuming that these GWs penetrate to higher altitudes, they then inferred corresponding GW horizontal velocities as large as ∼100 ms −1 or larger at altitudes of ∼250 to 300 km for GWs propagating against the tidal winds and refracting to higher intrinsic frequencies.Corresponding predicted electron density perturbations of ∼10 to 20% appear to agree with observed variations measured by digisondes and seen to occur in TIMED/GUVI tomographic reconstructions of plasma densities.This suggests general agreement between GW amplitudes and scales predicted by our present theoretical description of GWs arising from deep convection and direct measurements of plasma density fluctuations and inferred velocities at the bottomside F layer.

GW spectra for nearly impulsive vertical forcings
We showed in Sects.2-3 through comparison with the FL solutions that we can accurately reconstruct the GW field for a typical convective plume with σ t =12 min via ray tracing the excited GWs using an anelastic dispersion relation.In this section, we compare the FL solutions for a nearly impulsive vertical body force and its image to that from the ray traced solution.
Figure 17 shows the unaveraged GW zonal momentum fluxes, uw, for an impulsive, vertical body force and its image at z=50, 70, and 90 km from the left to right columns, respectively.These are the FL (Boussinesq) solutions obtained from Fourier-transforming Eqs. ( 10) and ( 12) to real space.The upper and lower rows show these winds at t=35 and 45 min, respectively.The vertical body forces which model this plume begin at t=0.Here, the vertical body force has a full width of D=20 km, a full depth D z =10 km, a duration of σ t =0.01 min, and an updraft velocity of w pl =40 m/s.As before, concentric rings of GWs are excited from these forces.Additionally, there are GWs with very high frequencies present in the spectrum, because there are concentric rings very close to the center of the forcing at z=90 km.If a GW propagates upwards with angle ψ to the vertical, then cos ψ=ω r /N in the Boussinesq approximation where H→∞ (e.g., Kundu, 1990).Since tan ψ=r/(z−z 0 ), all GWs at radius r=40 km and at z=90 km propagate at an angle ψ=26 o from the vertical in a windless, isothermal atmosphere, and have frequencies of ω r =0.9N.For the buoyancy period we assume here, this equates to a wave period of τ r =5.8 min.From Fig. 17c, this GW has λ H ∼35 km at t=35 min.From Eq. ( 9) with f =0, we calculate a very large vertical wavelength of λ z ∼72 km for this GW.Comparing with Fig. 5f, we see that for σ t =12 min, this high-frequency GW is not excited with a large amplitude because there are seemingly no waves at r∼40 km and z=90 km at this time.This difference in wave amplitude between the two Boussinesq solutions occurs because σ t ∼12 min is somewhat larger than the characteristic time scale of the forcing, which is (VF2001) For this forcing, τ c =9 min.Therefore, setting σ t =12 min results in the highest frequency portion of the GW spectrum to be substantially reduced in amplitude (VF2001).
In Fig. 18, we show the ray traced, real space, momentum fluxes for the GWs excited from the same impulsive body force in Fig. 17, and at the same altitudes and times as in Fig. 17.We see that there is a large region for small radii where GWs are excluded for the ray trace (anelastic) solutions.This is caused by the cutoff frequency in the anelastic GW dispersion relation (Marks and Eckermann, 1995); GWs are not allowed with frequencies larger than this cutoff frequency, which depends on λ H .At z=90 km, this region is r<70−80 km at t=35 min and r<50−60 km at t=45 min.This region decreases with time as the slower GWs with smaller λ H and λ z propagate to z=90 km, since GWs with smaller λ H have larger cutoff frequencies.For the GW in Fig. 17c with λ H ∼35 km and λ z ∼72 km, we recalculate ω r /N∼0.847 using Eq. ( 48), which is only slightly smaller than the Boussinesq value.The angle this GW makes with respect to the vertical, however, is: for a zonally-propagating GW in an isothermal, windless environment.Since for l=0 (using Eq. 48), Eq. ( 69) becomes Note that the anelastic correction term, 1/4m 2 H 2 , is of order one when λ z ∼2πH.For this GW, since the "correction" term is 0.67, the angle this GW makes with the vertical is ψ=39 o , leading to a radius of r∼70 km at z=90 km.This explanation corresponds well with the smallest concentric rings seen in Fig. 18c.
In Fig. 19, we directly compare the zonal and vertical velocity perturbations at z=70 km and y=0 at t=35 min (upper row) and 45 min (lower row).The solid and dash lines show the ray trace anelastic and FL Boussinesq solutions, respectively.We see that the amplitudes and phases are in excellent agreement for x>50 km.However, for smaller x (or r), which corresponds to frequencies near N, the FL model vastly overestimates the wave amplitudes.Again, this is because the FL model assumes that H→∞, which is not a good approximation for GWs with λ z ∼2πH.Therefore, the Boussinesq solutions do not accurately represent the solutions to deep vertical forcings when the forcings are nearly impulsive (with a duration smaller than the buoyancy period), because they result in large-amplitude GWs with large λ z and with frequencies near N.This artificial component of the GW Boussinesq spectrum results in large-amplitude concentric rings that are too close to the center of the forcing at airglow altitudes.

Summary and conclusions
In this paper, we incorporated GW phases into our convective plume and ray trace models in order to reconstruct the Fig. 19.Ray traced GW amplitudes as compared to FL amplitudes at z=70 km and y=0.The first and second columns show the GW perturbation fields u and w, respectively.The upper and lower rows show the results at t=35 and 45 min, respectively.Solid lines show the ray trace results, while dash lines show the FL Boussinesq results.
GW perturbation fields at higher altitudes for deep convective plumes.We represented a single convective plume as an upward-moving vertical body force plus a downwardmoving image vertical body force in order to take into account the Earth's surface and reflection of the GWs off the ground.We found that for this single convective plume, the reconstructed wave fields from ray tracing result in concentric GW rings at the OH airglow layer in an isothermal, windless atmosphere.We found that ray tracing the upwardpropagating GWs from the reflected GW spectrum, rather than ray tracing the upward and downward-propagating GWs from the unreflected GW spectrum, yields the best results, and agrees very well with the exact FL solutions, which are the Boussinesq solutions to vertical body forces in an isothermal, windless atmosphere.Via intercomparision between the ray trace and FL solutions, we determined the normalization constant needed to convert the spectral, ray traced momentum fluxes to the real-space momentum fluxes.Because the ray trace model is generalizable to non-isothermal temperatures and non-zero winds, the formalism developed in this paper allows for the determination of GW effects in the MLT and TI for more realistic atmospheric environments (see Vadas et al., 2009b).
Additionally, via ray tracing a few, individual GWs and comparing their locations and times with the ray trace solution for the entire convective plume, we determined the functional form of the GW spectra needed to easily calculate the GW amplitudes in the MLT and TI from a deep convective plume.We then calculated the zonal velocity amplitudes of GWs excited from a single convective plume and a small convective cluster at the OH layer as a function of λ H and λ z .These spectra were used to infer GW amplitudes at higher altitudes in the TI during this SpreadFEx campaign (Fritts and Vadas, 2008;Fritts et al., 2008b).They were also used to constrain the amplitudes of medium-scale GWs detected in the OH airglow layer and reverse ray traced to deep convective plumes (Vadas et al., 2009a).
Finally, we found that the Boussinesq solutions are not accurate for deep, impulsive vertical body forcings, because the solutions result in large-amplitude GWs at unrealistically small radii at much higher altitudes.This is not unexpected, considering the inaccuracy of the Boussinesq www.ann-geophys.net/27/147/2009/Ann.Geophys., 27,2009 dispersion relation for GWs with frequencies near the buoyancy frequency (Marks and Eckermann, 1995).

Fig. 1 .
Fig. 1.(a) Vertical body force, F z , representing a single convective plume with no reflection off the Earth's surface.(b) GW zonal velocity perturbation from the body force in (a) at t=20 min and y=0.(c) GW vertical velocity perturbation from the body force in (a) at t=20 min and y=0.(d-f) Same as (a-c), but for a single convective plume with reflection off the the Earth's surface.Contours for F z are in intervals of 10% of the maximum value.Contours for the perturbation velocities are 10, 30, 50, 70, and 90% of the maximum value.Solid lines indicate positive values, and dash lines indicate negative values.The maximum values in (b), (c), (e),and (f) (unscaled, with constant background density) are 0.5, 0.5, 0.9, 0.4 m/s, respectively.The resolution is x=4.4 km.

Fig. 2 .
Fig. 2. Vertical flux of zonal momentum in the form given by Eq. (35) as a function of horizontal and vertical wavenumbers (left column) and wavelengths (right column) for upward and eastward-propagating GWs (solid lines) excited from a single convective plume.Upper row: No ground reflection.Lower row: Ground reflection included.Contours show 10% intervals of the maximum value.Pink dash-dot lines indicate c gz in intervals of 15 m/s.Blue dash lines indicate c H in intervals of 50 m/s.Dotted lines in (a) and (c) show λ H =13.3, 17.8 km, and 26.6 km.The resolution is x=6.7 km.

Fig. 3 .
Fig. 3. Row 1: Single convective plume represented by a vertical body force and its image.Row 2: Small convective cluster represented by 3 vertical body forces and their images.Vertical slices of these body forces are shown at y=0 in (a) and at y=9 km in (d).Horizontal slices of these body forces at z=11 km are shown in (b) and (e).GW zonal velocity perturbation contours are shown at t=20 min and y=0 in (c) and (f).Body force contours are in intervals of 10% of the maximum value.Contours for the perturbation velocities are shown at 10, 30, 50, 70, and 90% of the maximum value.Solid lines indicate positive values, and dash lines indicate negative values.The maximum values in (c) and (f) (unscaled, with a constant background density) are 0.9 and 1.2 m/s, respectively.The resolution is x=4.4 km.

Fig. 5 .
Fig. 5. Horizontal slices of the vertical flux of zonal momentum for a single convective plume with reflection using the FL solutions Fourier transformed to real space.The left, middle, and right columns show the GW response at z=50, 70, and 90 km, respectively.The upper to lower rows show the images at t=25, 35, 45, and 55 min, respectively.Maximum values (scaled for the decreasing background density) for each image are approximately twice those values shown in Fig. 10a-c as diamonds.Maximum positive values are white, while maximum negative values are black.The resolution is x=6.7 km.

Fig. 6 .
Fig. 6.Horizontal slices of u, v, and w in columns 1, 2, and 3, respectively at z=70 km for the same convective plume in Fig. 5. Rows 1-4 show the images at times t=25, 35, 45, and 55 min, respectively.From left to right, the maximum values (scaled for the decreasing background density) for the u and w images are: 1st row: 8 and 6 m s −1 .2nd row: 15 and 12 m s −1 .3rd row: 18 and 14 m s −1 .4th row: 17 and 11 m s −1 .The maxima for the v images are identical to those of the u images.

Fig. 8 .
Fig. 8. Zonal flux of vertical momentum at z=90 km and t=45 min from ray tracing the excited GW spectra from a single convective plume with reflection.The GW spectra used in (a), (b), (c), and (d) are those shown in Fig. 7a, b, c, and d, respectively.

Fig. 9 .
Fig.9.The vertical flux of zonal momentum for the GWs ray traced from the same convective plume as in Fig.5.The times and altitudes are the same as in Fig.5, as labeled.Here we use the reflected spectrum in Fig.2c-d, and only ray trace those upward-propagating GWs.This plot is unsmoothed in x, y, z, and t.The resolution for the spectra is x=8.9 km.

Fig. 11 .
Fig. 11.Horizontal slices of the reconstructed u, v and w fields from ray tracing at z=70 km and t=55 min.Upper row shows the results using the unreflected GW spectrum (Fig.2a-b), while the lower row shows the results using the reflected GW spectrum (Fig.2c-d).From left to right, the maximum values of u, v, and w are: 1st row: 17, 17 and 11 m s −1 .2nd row: 14, 14 and 9 m s −1 .The resolution is x=8.9 km.

Fig. 12 .
Fig. 12. Horizontal slices of the reconstructed u, v and w fields in columns 1, 2, and 3, respectively, from the solutions shown in Fig. 9 at z=70 km.Here, we show the results at t=25, 35, 45, and 55 min in the upper to lower rows, respectively.From left to right, the maximum values of u and w are: 1st row: 3 and 2 m s −1 .2nd row: 15 and 11 m s −1 .3rd row: 17 and 12 m s −1 .4th row: 14 and 9 m s −1 .The maxima for the v images are identical to those of the u images.

Fig. 13 .
Fig.13.Ray traced GW amplitudes as compared to FL amplitudes at z=70 km and y=0 for the same convective plume.The first, second, third, and forth columns show the GW perturbation fields u, w, 100 T /T , and 100 ρ /ρ, respectively.Upper to lower rows show the results at t=35, 45, and 55 min, respectively.Solid lines show the ray trace results, while dash lines show the FL Boussinesq results.

Fig. 14 .
Fig. 14.(a) Solid, dot, dash, and dash-dot lines show λ x from ray tracing each individual GW from left to right in Table 2, respectively.Diamonds, triangles, squares, and x's show the average values of λ x in the ray trace bins closest to the individual GWs from left to right in Table 2, respectively.(b) Same as in (a), but for λ z .(c) Average momentum fluxes times (z−z 0 ) 2 exp(−(z−z trop )/H) for the bins closest to the individual GWs.Symbols are the same as in (a).

Fig. 15 .
Fig. 15.(a-b) Same as in Fig. 14a-b, but for the GWs in Table 3. (c) Values of |u|(z−z 0 ) exp(−(z−z trop )/2H) for the ray trace bins closest to the individual GWs in a-b.Symbols are the same as in (a).

Fig. 16 .
Fig. 16.GW horizontal velocity amplitudes at z=87 km in intervals of 5 m/s (light solid lines).These are the maximum amplitudes at this altitude.(a) Single convective plume with ground reflection.(b) Small convective cluster with ground reflection.Pink dash-dot lines indicate c gz in 15 ms −1 intervals.Blue dash lines indicate c H in 50 ms −1 intervals.

Fig. 17 .
Fig. 17.Horizontal slices of the vertical flux of zonal momentum for an impulsive vertical body force and its image using the FL solutions Fourier-transformed to real space.For this forcing, D=20 km, D z =10 km, and σ t =0.01 min.The left, middle, and right columns show the GW response at z=50, 70, and 90 km, respectively.The upper and lower rows show the solutions at t=35 and 45 min, respectively.The resolution is x=6.7 km.

Fig. 18 .
Fig. 18.Horizontal slices of the vertical flux of zonal momentum for the GWs ray traced from the same convective plume as in Fig.17.Times and altitudes are the same as in Fig.17.This plot is unsmoothed in x, y, z, and t.The resolution for the spectra is x=8.9 km.

Table 1 .
Symbols and notation., ... Fourier transform of u, v, F z , ... u GW , v GW , ... The GW portion of the solution, u, v, ..., after the force is finished φ phase of wave in ray trace model (x 0 , y 0 , z 0 ) center of the convective plume envelope σ x , σ y , σ z half widths at half-max of the vertical body force in x, y, and z D=4.5σ x , D z =4.5σ z full widths and depths of the body force is employed to advance the ray equations and Eq.(45) in time.

Table 2 .
Individual GW parameters for verifying uw * spectrum.

Table 3 .
Individual GW parameters for verifying u H spectrum.