## Abstract

In recent years, aircraft data from mountain waves have been primarily analyzed using velocity and temperature power spectrum and momentum flux estimation. Herein it is argued that energy flux wavelets (i.e., pressure–velocity wavelet cross-spectra) provide a more powerful tool for locating and classifying different types of mountain waves. In the first part of the paper, pressure–velocity cross-spectra using various linear mountain-wave solutions are shown to be capable of disentangling collocated waves with different propagation directions and wavelength. A field of group velocity vectors can also be determined.

In the second part, the energy flux wavelet technique is applied to five cases of mountain waves entering the stratosphere from the Terrain-Induced Rotor Experiment (T-REX) in 2006. Perturbation pressure along the flight track is determined using aircraft static pressure corrected hydrostatically with GPS altitude. In four of the cases, collocated long up-propagating and short down-propagating waves are seen in the stratosphere. These waves have strong, but opposite, *p*′*w*′ cospectra. In one of these cases, a patch of turbulence is collocated with the up and down waves. In two other cases, trapped waves riding on the tropopause inversion layer (TIL) are seen. These trapped waves have *p*′*w*′ quadrature spectra that reverse sign across the tropopause. These newly discovered wave types may arise from secondary wave generation (i.e., a nonlinear transfer of energy from the long vertically propagating waves to shorter modes).

## 1. Introduction

While balloon soundings and remote sensing have contributed to mountain-wave observation, the most detailed observations derive from horizontal flight legs of research aircraft. Aircraft mountain-wave projects have been carried out in many regions of the world including the Alps, the Pyrenees, the Rockies, and the Sierra Nevada. The three major innovations in measurement technology in aid of these measurements were 1) gliders with a recording variometer, 2) aircraft with temperature and pressure instruments, and 3) aircraft with a gust probe and inertial platform. The mathematical analysis of aircraft data has advanced also, from using potential temperature to chart streamline deflection to using longitudinal and vertical velocity to compute momentum flux and Fourier cross-spectra to determine dominant wave scales and phase relationships. This history probably begins with Kuettner’s study of stationary updrafts and downdrafts in mountain waves using a sailplane variometer over the Riesengebirge in Sudeten (Kuettner 1938).

Eliassen and Palm’s (1960, hereafter EP60) theoretical description of energy and momentum fluxes within stationary mountain waves laid a foundation for modern observations of momentum fluxes. Eliassen and Palm showed that the vertical flux of horizontal momentum (**MF**) is constant with height regardless of atmospheric conditions so long as the wave is steady, linear, propagating, and nondissipative. The vertical energy flux (EF) varies with height proportional to the wind speed and the two fluxes are related by the formula EF = −*U* · **MF**. Bretherton (1969) and others showed how waves might be absorbed in turbulent patches and critical layers.

Lilly and Kennedy (1973, henceforth LK73) first measured wave momentum fluxes over the Rockies, using aircraft leg altitudes all the way up to 60 hPa. While LK73’s measurements were generally consistent with EP60, data quality issues limited the interpretation of their measurements. Gust probes and inertial data systems were still in their infancy and were unavailable on most of the research aircraft. Instead, wind vertical velocity was computed from the potential temperature field that was computed using the static pressure and temperature measurements. LK73 calculated momentum fluxes using three different methods and found a generally constant negative momentum flux that was seemingly absorbed in a turbulent stratospheric layer associate with wave breaking. LK73 also made the first attempt to use Fourier analysis on mountain-wave data and had limited success. LK73’s momentum flux results were soon confirmed in another study looking at a different flight within the same field campaign (Lilly and Lester 1974). Brown (1983) further expanded upon LK73 by examining momentum fluxes related to trapped waves.

Georgelin and Lott (2001, henceforth GL01), following Bougeault et al. (1990), examined momentum fluxes in one case of trapped lee waves over the Pyrenees. Power spectra in GL01 showed more than one wave packet at an altitude of about 4 km with the shorter-wavelength packets evanescing with height so that the dominant wavelength changed from 8.5 km at an altitude of 4 km to a dominant wavelength of 15.5 km at an altitude of 6 km. GL01 interpreted their results to indicate the superposition of two wave packets at low levels. By applying bandpass filters, GL01 showed that trapped waves can contribute significantly to wave momentum fluxes.

Gravity waves over the Alps appear to be more complicated than over the Rockies or Pyrenees (Smith et al. 2002; Smith and Broad 2003; Smith et al. 2007). The complex terrain, directional wind shear, and low-lying critical levels there often absorb waves in the middle troposphere. Boundary layer absorption prevents resonant trapped waves; momentum flux distributions are patchy with occasional reversals.

The most recent contribution to the observational literature on mountain waves was the study of waves entering the stratosphere over the Sierras in California during the Terrain-Induced Rotor Experiment (T-REX; Smith et al. 2008, henceforth S08). Momentum fluxes were mostly negative, with one striking exception in research flight 10 (RF10). Using GPS-altitude-corrected static pressure, mountain-wave energy fluxes were observed for the first time. All cases, with negative or positive fluxes, satisfied the Eliassen–Palm relation between energy and momentum flux. Using Fourier pressure–velocity cross-spectra, three types of waves were identified: long upward propagating waves, short downward propagating waves, and short trapped waves near the tropopause. Finding these short waves in one event raises a question about the frequency of such waves.

Another aspect of the T-REX wave studies (S08) was the attention paid to background layering. If the background wind profile is layered or sheared and that nonuniform flow is then displaced vertically by a mountain wave, an aircraft flying at a constant altitude through that wave field would sample horizontal velocity perturbations resulting from the background layering rather than wave-induced perturbations. Pressure perturbations do not suffer the same drawback and can be much more cleanly measured through a wave field.

The major limitation in S08 was the lack of an analysis method to distinguish wave packets with different locations, wavelengths, and propagation directions. A shortcoming of their Fourier spectral analysis was the lack of spatial information. Wavelet transforms can be used to simultaneously present location and wavelength information. Through the use of complex wavelets, the continuous wavelet transforms of two variables, such as *p*′ and *w*′, can be combined to create a cross-spectrum, in this case representing vertical energy flux. Wavelet energy flux calculations have not previously been employed in aircraft mountain-wave studies, although complex wavelet analysis of satellite brightness temperatures (Alexander and Barnet 2007; Alexander et al. 2009), two-dimensional numerical simulations (Alexander 1996), and surface pressure (Hauf et al. 1996; Fritts et al. 2003) have been used to diagnose gravity wave propagation.

In this study, we develop a wavelet energy flux analysis method and apply it to an extended set of Gulfstream V (GV) observations from T-REX and idealized model results (Table 1). There are GV data available from 13 T-REX flights during March and April 2006. In addition to the race track flight legs over the Sierras, we include long ferry legs between Colorado and California that also contain mountain-wave events.

## 2. Methodology

For our purposes, a continuous wavelet transform and nonorthogonal complex wavelet function are needed to calculate vertical fluxes using wavelet cross spectra. The most commonly used complex mother wavelet, and the one used in this study, is the Morlet wavelet, which is a plane wave whose amplitude is modulated by a Gaussian. It resembles the pattern expected in a resonant mountain wave. The Morlet wavelet has the form

where *ω*_{0} is the wavenumber and *ζ* is the nondimensional spatial coordinate. To satisfy the admissibility condition, which requires that the wavelet have a 0 average, we set *ω*_{0} = 6 to ensure that otherwise necessary correction terms are of the same order as typical numerical rounding errors (Farge 1992). The continuous wavelet transform at some position index *n* and scale *s* is defined as the convolution of the series of interest *y* with the complex conjugate of the wavelet, which is scaled to have unit variance by

where Δ*x* is the horizontal spacing of the data. The choice of a wavelet with a number of oscillations within the wave packet provides superior frequency resolution at the expense of spatial resolution. With such a spatially broad wavelet, we expect to see wavelet transform signatures that are more spread out than the signal being analyzed.

Because we are dealing with a noncyclic spatial dataset of finite length, boundary influences will induce errors near the edge of the domain that will propagate further into the interior of the domain at longer scales. Transforms and boundary influences were processed using scripts modified from a package written by Torrence and Compo (1998, henceforth TC98) and distributed with the National Center for Atmospheric Research (NCAR) Command Language (NCL; see acknowledgments section).

Aircraft data were linearly interpolated from the 1-Hz data to similar spatial transects with 200-m resolution. Aircraft transects were required to have minimal heading or altitude changes, which was nearly universal due to the experiment design. No smoothing was performed on the data itself. Detrending removed the series mean and least squares trend for Fourier and wavelet analysis.

The detrended series which subjected to a 10% split-cosine-bell taper before Fourier spectra were calculated. The discrete Fourier transform of the detrended and tapered series *y*′ is a convolution with sinusoidal functions:

The periodogram was subjected to a three-point modified Daniell smoother to account for spectral leakage. The smoothed transform, *Y*′* _{k}*, was then scaled to conserve the variance of the detrended and tapered series by

Energy and momentum fluxes were calculated explicitly using the horizontal velocity and pressure perturbations. The equations for calculating the leg-integrated momentum and energy fluxes explicitly from velocity (*w*′) and pressure (*p*′) perturbations and average air density (*ρ*) are

The units of (5) and (6) are Newtons per meter and watts per meter, respectively. To calculate the energy flux in Eq. (6), the static pressure (*p*) from aircraft data must be corrected for aircraft altitude fluctuations (*z*) detected from GPS, using the formula

where the reference altitude *Z*_{REF} is taken to be the average altitude of the flight leg. To analyze wave flux we expressed the *p*′*w*′ cross spectrum using real and imaginary parts. The energy flux, showing wave propagation, was calculated in both wavelet and Fourier space by

where *P̃ _{n}*(

*s*) represents the wavelet transform of the pressure perturbation and

_{j}*W̃**

_{n}(

*s*) is the complex conjugate of the transform of the vertical velocity perturbation each at spatial index

_{j}*n*, with

*s*being the wavelet scale at wavenumber

_{j}*j*. Conversely, the

*p*′

*w*′ quadrature, indicating wave trapping, is {

*P̃*(

_{n}*s*)

_{j}*W̃**

_{n}(

*s*)} and {

_{j}***

_{k}Ŵ_{k}}. Significance was calculated for the 95% confidence interval (see appendix A). The single-variable power spectra for vertical velocity are |

*W̃*(

_{n}*s*)|

_{j}^{2}and |

*Ŵ*|

_{k}^{2}in wavelet and Fourier space.

Total energy flux can also be calculated from either the Fourier or wavelet representation using Parseval’s theorem:

where Δ*j* is the wavenumber resolution used in the wavelet calculations, and *C _{δ}* = 0.776 is a reconstruction factor that is constant and unique for each mother wavelet as explained in TC98. The Morlet wavelet scales,

*s*= Δ

_{j}*x*2

^{jΔj+1}, were adjusted to match Fourier wavelengths,

*λ*= 4

*πs*/(

_{j}*ω*

_{0}+ ), in all discussions and displays of wavelet scale and wavelength. The EF given in (6), (10), and (11) is given in Watts per meter of ridge length and agree to within about 1%–2%. In each case momentum fluxes were calculated analogously by swapping the horizontal velocity perturbation for the pressure perturbation and scaling by the density.

## 3. Queney solution

The simplest case of a vertically propagating mountain wave is the Queney solution for a hydrostatic nondispersive, linear wave in a single-layer atmosphere using the Boussinesq approximation (Queney 1948). The lower boundary in the Queney solution is a Witch of Agnesi ridge defined by

Parcel displacement is given by

where *h _{m}* is the mountain height,

*a*is the mountain half-width, and

*l*= |

*N*/

*U*| is the Scorer parameter, which is a ratio of the buoyancy frequency to the wind speed (Scorer 1949). The wind and pressure fields can then be solved for using the kinematic relation

*w*′ =

*U*(∂

*η*/∂

*x*), continuity

*u*′ = −

*U*(∂

*η*/∂

*z*), and the Bernoulli relation

*p*′ = −

*ρ*

_{0}

*Uu*′, yielding

where *ρ*_{0} = 1.292 kg m^{−3} is the reference density. The Queney solution for an atmosphere with *U* = 10 m s^{−1}, *N* = 10^{−2} s^{−1}, *a* = 7 km, and *h _{m}* = 1000 m is shown in Fig. 1a. The momentum and energy fluxes for the Queney solution are

**MF**= −(

*π*/4)

*ρ*

_{0}

*NUh*

_{m}^{2}and EF = (

*π*/4)

*ρ*

_{0}

*NU*

^{2}

*h*

_{m}^{2}. These analytic solutions can be evaluated on a discrete grid and analyzed using our wavelet methods. Wavelet analysis performed at numerous altitudes indicates that wavelet power and cospectra are independent of altitude with errors on the order of 10

^{−4}for single variable spectra and 10

^{−5}or less for cross-spectra. The error present in the Queney solution provides an upper bound of the accuracy that can be expected in experimental analysis.

An example of the vertical velocity and pressure transforms is given in Fig. 2, with the associated cross-spectrum shown in Fig. 3. As can be seen from Eqs. (14)–(16), the vertical velocity function is spatially more compact and therefore shifted to shorter scales than the pressure or horizontal velocity spectra. The cross-spectrum, in turn, has a maximum power at a scale in between the vertical velocity and pressure spectra maxima. Only one altitude is arbitrarily shown for the Queney solution since the wavelet power and cross-spectra signatures are independent of altitude. As the amplitude of the Queney solution scales linearly with the height of the mountain, the magnitude of the wavelet power scales with the square of the height of the mountain. While vertically propagating waves have no quadrature spectrum in Fourier space, the same is not true in wavelet space where the waves show offsetting areas of weak positive and negative quadrature over the mountain slopes (see appendix B).

## 4. Linear model

### a. Linear case 1: Two isolated hills

To test the effectiveness of our wavelet techniques on idealized flow, we used a three-layer, two-dimensional, dispersive linear model based on the equations from Smith (2002) and Smith et al. (2002) to generate theoretical mountain-wave patterns to compare with the Queney solution. The simplest case that we tested with the linear model is a single-layer atmosphere over two isolated hills with nonhydrostatic mountain waves. Two Witch of Agnesi ridges with half-widths of 10 and 2.5 km, positioned 150 km upwind and downwind of the center of the domain respectively, were computed in a 3200-km-wide domain at a horizontal grid spacing of 400 m. This inner domain was imbedded in an outer domain of 32 678 grid points or 13 107.2 km for the purpose of Fourier calculations. Other parameters used were set to the same as those used in the Queney solution: *U* = 10 m s^{−1}, *N* = 10^{−2} s^{−1}, and *h _{m}* = 500 m. The displacement field for linear case 1 is shown in Fig. 1b.

The wavelet cospectrum and quadrature spectrum of the domain containing the two Witch of Agnesi mountains are shown in Fig. 4. The two wave packets are clearly separable in their wavelet signatures. Both mountains show a single area of positive *p*′*w*′ cospectrum and weak quadrature associated with the purely up-propagating waves. The wavelet signature for the broad mountain is consistent with that expected from the Queney solutions. The narrower mountain shows a skewed cross-spectrum. With nonhydrostatic dispersion, the wavelet signature for the narrower mountain is skewed slightly into a “chirp” pattern with rising wavenumbers on the downwind side. The onset of dispersion is controlled by the Scorer parameter.

### b. Linear case 2: Targeted reflection of dispersive waves

A more complex mountain-wave pattern was generated to test the ability of wavelets to distinguish collocated up- and down-propagating waves with different wavelengths (see observations in section 6). Two separate topographic features were used: a Gaussian-modulated cosine terrain profile, which itself is in the form of a Morlet wavelet, positioned 50 km upwind of a single smooth Gaussian ridge. Both terrain features were set with maximum amplitudes of 1 km and half-widths of 20 km. The wavelength of upstream cosine terrain was adjusted to target the ray path. The Gaussian ridge and cosine-Gaussian are calculated using

respectively, where *x _{c}* is the center coordinate of the cosine-Gaussian terrain,

*a*is the half-width of the Gaussian,

*h*is the amplitude of the Gaussian mountain, and

_{m}*k*is the wavenumber. To aim the reflected wave, we use the slope of the ray path (Smith 2002):

where *U*_{0} and *N*_{0} are the wind speed and buoyancy frequency in the lower layer, *z _{d}* is the total vertical distance traveled by the wave, and

*x*is the distance downstream that the wave will travel. For example, this linear model case was calculated with

_{d}*N*

_{0}= 0.01 s

^{−1}and

*U*

_{0}= 10 m s

^{−1}. We wanted to spatially collocate the up-propagating and down-propagating waves so that the first reflected wave hit at

*x*= 50 km downstream after traveling up and then down 13 km (for a total of

_{d}*z*= 26 km), so we used a terrain wavelength

_{d}*λ*= 2

*π*/

*k*= 7.1 km.

A partial reflector at *z* = 13 km was used to reflect the shorter up-propagating wave packet. In our case, we decreased the Scorer parameter by one-half by increasing the wind speed above 13 km. While partial reflection is also expected and observed for the longer wave, it is much weaker since the upper Scorer parameter is set to a value that allows propagation of the longer wave. To isolate the down-propagating wave packets, we established complete absorption at the lower boundary condition so that no wave energy would reflect off the ground. For details on the application of the lower boundary, see Smith et al. (2002). Figure 1c shows the displacement pattern for the partial reflection case with total surface absorption.

By taking a transect at an altitude of 2 km, we can observe the wave field just above the ridge top but close enough to the ground so that the up-propagating wave’s dispersion is minimal and the down-propagating waves have dispersed enough to allow spatial separation. Figure 5a shows the low-level vertical velocity wavelet power spectrum, which indicates two distinct shortwave packets near the center of the domain. One of the packets is centered over the cosine-Gaussian terrain and the other is centered over the Gaussian ridge. The longer wave launched from the Gaussian ridge is below the shading scale in this case because of the high variance from the dispersive waves, which match very closely the form of the Morlet wavelet.

Figure 5b shows the *p*′*w*′ wavelet cospectrum, with the sign indicating the direction of vertical energy transport. Two up-propagating patches and one down-propagating patch are seen. The up-propagating long wave centered over the Gaussian ridge is now apparent. The two short wave packets evident in the vertical velocity power spectrum (Fig. 5a) have opposite signed energy fluxes, with the up-propagating wave located over the cosine-Gaussian terrain and the down-propagating wave carrying less energy downward over the Gaussian ridge. There is no subsequent up-propagating wave because of the absorbing lower boundary condition.

In this setup, both the short wave and long wave are free to propagate in the lower layer and the reflections above and below are controlling. When the Scorer parameter is reduced by a factor of 0.5 in the upper layer, the short wave evanesces because its horizontal wavelength of ∼6.5 km is shorter than the critical value of ∼8.9 km set by the Scorer parameter in the upper layer. On the other hand, the longer wave with a horizontal wavelength of ∼60 km is long enough to freely travel across the interface. Additionally, while resonant modes are possible in this configuration, their amplitude would be small because of the large thickness of the lower layer and the absorptive lower boundary.

### c. Linear case 3: A resonant trapped wave

In this final example, we ran a two-layer simulation of a trapped wave with the parameters used for the trapped wave in Smith (2002). The solution was computed with a Gaussian ridge height of 500 m and half-width of 2 km. The layer interface was at an altitude of 2 km and the lower boundary has a reflection coefficient of 0.9. The lower and upper layers were set to wind speeds of 10 and 22 m s^{−2} and buoyancy frequencies of 12 × 10^{−3} and 8 × 10^{−3} s^{−1}, respectively. This resulted in a trapped wave that decayed gradually downstream (Fig. 1d). The wave node for this trapped wave was at about *z* = 1.4 km. At its node the wave has a *w*′ maximum and *p*′ and *u*′ are 0. We sampled the wave field in the lower layer at *z* = 1.0 km and in the upper layer at *z* = 1.75 km.

The *p*′*w*′ cospectrum has two distinct regions of *p*′*w*′ power (Figs. 6a,c). The patch of positive *p*′*w*′ power over the mountain ridge is associated with the up-propagating mountain wave while the tail of negative cospectrum downstream shows the decaying trapped wave. The slight negative energy flux in the lee wave is due to the slight absorption at the lower boundary (i.e., *q* = 0.9). It is worth noting that the trapped wave is spectrally narrower than the up-propagating wave due to the nature of the resonant wave.

The *p*′*w*′ quadrature spectrum associated with the trapped wave exhibits a long strong tail of quadrature. The sign of the quadrature signal indicates that the trapped wave is being sampled above (positive tail) or below (negative tail) its node, which is at *z* = ∼1.4 km. Figure 6b shows a positive tail at *z* = 1.75 km, indicating that we are sampling above the node. Below the node at *z* = 1.0 km (Fig. 6d) the quadrature tail is negative (see appendix C).

## 5. Wavelet-derived group velocity

As a further experiment, we attempted to use wavelet energy fluxes to derive the field of fluid-relative group velocity in linear case 3. As there is no wind shear within each layer, we can assume that the group velocity is the ratio of the pressure-velocity work to the wave energy density. The issue of wave action need not be considered (Bretherton 1966). Since wavelets draw information from more than one period, we can calculate the group velocity in wavelet space by

where *u*′* _{i}* is a single directional component of the wind speed. The group velocity diagram for the middle of the waveguide layer at

*z*= 1 km is shown in Fig. 7, where areas of insignificant energy flux are masked. The vectors represent the fluid-relative horizontal group velocity and the color shading represents the vertical group velocity. The vertically propagating wave is moving upward at approximately 1 m s

^{−1}and has a leftward group velocity of nearly 10 m s

^{−1}. The trapped lee wave has a very small downward group velocity and an upwind group velocity of about 5 m s

^{−1}.

The horizontal components can be understood as follows. With a mean flow speed of 10 m s^{−1}and a leftward fluid-relative group velocity of 10 m s^{−1}, the horizontal earth-relative group velocity is 0. As expected, the group velocity vector for the long hydrostatic waves is vertical. The shorter trapped wave has a much smaller but still upwind fluid-relative group velocity. This indicates that the wave energy is propagating upwind with respect to the flow but downwind with respect to the ground. In this way the length of the lee wave increases with time.

While the group velocity can be calculated using wavelets with the linear model output, doing the same for real data is more challenging. Equation (20) will not hold for an atmosphere with rapidly changing *U*, as changes in the intrinsic frequency will change the energy density in order to conserve wave action. To calculate the group velocity, the kinetic energy flux must be normalized by the energy density, which is the sum of the kinetic and potential wave energies. In the linear model, the potential energy can be calculated using the wave displacements and the layer buoyancy frequency. However, because the real atmosphere does not have a constant stability profile, our inability to accurately calculate the potential energy prevents us from calculating the group velocity magnitude even with a constant wind speed, although the propagation direction can still be found using a ratio of the kinetic energy flux components.

## 6. Mountain waves from T-REX

### a. Collocated waves in RF10/IOP 13

S08 documented a case of net negative energy fluxes in the stratosphere of T-REX research flight 10/intensive observation period 13 (IOP 13). In this case, a down-propagating wave of 20-km wavelength was observed in Fourier spectra while a longer 30–40-km up-propagating wave disappeared as it crossed the tropopause. In the stratosphere the net fluxes were reversed. In addition, a trapped wave with a horizontal wavelength of 12 km was observed riding along the tropopause. In this section we repeat that analysis with wavelets.

The Fourier cross-spectrum components of *p*′ and *w*′ for the first aircraft transect at each vertical level on the south leg of IOP 13 and all parcel displacement curves for all aircraft transects are shown in Fig. 8. Each transect shows the down-propagating wave (negative *p*′*w*′) present with a 20-km wavelength. Likewise, the lower two transects both show the 30–40-km upward propagating wave, which unexpectedly is missing at the highest level. Instead, at the highest level we see a ∼12-km oscillation with no cospectrum and negative quadrature, indicative of a trapped wave being sampled below its node. Unfortunately, these Fourier techniques provide no spatial information. Wavelets, on the other hand, present a useful tool for examining the spatial and spectral structure.

Using the wavelet techniques from section 4, we can reveal added details about the structure of the wave field (Figs. 9 and 10). Figure 9 gives the pressure and vertical velocity cospectrum for the same aircraft flight legs shown in Fig. 8. The relatively short length of the T-REX flight legs becomes an immediate concern because of the encroaching boundary effects. The up-propagating waves present at the lowest two levels are barely able to emerge from the “cone of influence” (COI) and so we are limited in our ability to make any statements regarding their presence any farther east than 20 km from the center of the valley, which roughly corresponds to the position of the Inyo Mountains.

There is clearly a down-propagating wave present over the Owens Valley with a wavelength of about 20 km. This down-propagating wave is evident and statistically significant in the same position at all three flight levels. At *z* = 13 km, the down-propagating energy exceeds the magnitude of the up-propagating energy, leading to the observed flux reversal (S08). Referring to both Figs. 8 and 9 it is clear that the wave with a horizontal wavelength of 12 km in the lower stratosphere does not appear strongly in the *p*′*w*′ cospectrum and therefore is carrying little energy flux.

Both the 12- and 20-km wavelength signals appear strongly in the vertical velocity spectra of the lower stratosphere in Fig. 10. There is just a slight hint of the 30–40-km wavelength up-propagating wave, which is insignificant at all levels. The presence of the 12-km wavelength signal in the vertical velocity and quadrature spectra but not in the energy flux spectrum leads us to classify this as a trapped wave being sampled below the wave node, possibly riding the tropopause inversion layer (TIL). In total, three different wave systems were present on this day over the Sierras.

The Scorer parameter profile derived from balloon data is shown in Fig. 11 and supports our interpretation of three distinct wave packets. If a 20-km wavelength down-propagating wave was generated in the stratosphere, it would be free to propagate down to the three flight levels of the GV, but not any farther, before it evanesced. Likewise, the 12-km trapped wave could propagate at *z* = 13 km and above but would evanesce below that level. While we can see the lower reflector of a potential wave duct, the limited depth of our balloon sounding precludes us from making any statement about a potential upper reflection in IOP 13b that would be needed for wave ducting (Lindzen and Tung 1976; Fritts and Yuan 1989). However, the reason for the disappearance of the 30-km wavelength up-propagating wave is not known. The up-propagating wave should be free to propagate into the stratosphere according to linear theory.

### b. More down-propagating waves in IOP 6

For a second example of collocated up- and down-propagating waves, we examine RF5 in T-REX IOP 6 on 25 March 2006 over the Sierra Nevada. Coincidently, IOP 6 was one of the most interesting cases of the entire campaign for those studying rotor circulations. The beginning of the event was characterized by a rotor circulation that transitioned into a severe downslope windstorm in the evening (Grubišić et al. 2008).

Figure 12 shows four consecutive aircraft flight legs along two parallel flight paths, which all show a down-propagating wave with a horizontal wavelength of about 20 km collocated with an up-propagating wave over the Owens Valley. Even the smallest magnitude of the down-propagating waves observed in these four flight legs is statistically significant at the 95% level and lies well outside the COI. The long wave up-propagating patch extends into the COI but has a strong part clear of boundary influences. The down-propagating EF is less than half the up-propagating EF. The vertical velocity wavelet power spectrum (not shown) shows a shorter wavelength wave collocated with the broad vertically propagating wave over the valley.

### c. Trapped waves in IOP 10

There are numerous examples of trapped waves present in the T-REX dataset, both over the Sierras and en route between California and Colorado. The *p*′*w*′ cross-spectra patterns for these trapped waves cases agree qualitatively with the linear simulations of trapping discussed earlier in this paper. Figure 13 shows an example of a trapped wave riding the TIL during RF 8 during T-REX IOP 10. The vertical velocity power (Fig. 13c) shows a spatially long but spectrally isolated streak of power, which confirms the presence of a trapped wave with a horizontal wavelength of ∼13 km. The vertical velocity wavelet spectrum also narrows in frequency with increasing downstream distance as the atmosphere selects out a natural resonant scale. The similarity in shape of the Morlet wavelet and a trapped wave contributes to the very clean representation of trapped waves.

The *p*′*w*′ cospectrum (Fig. 13a) pattern is very different than the *w*′ power. The dominant patch is an up-propagating patch at *x* = −10 km at the same wavelength as the *w* power. This patch is very similar to the patch seen in the ideal lee wave case (linear case 3). It too is followed by a long weak patch of down-propagating wave energy. The *p*′*w*′ quadrature spectrum and *w*′ power also confirm the expectations from the linear model. The quadrature spectrum makes a spatially elongated and spectrally confined streak that, when combined with the cospectrum, shows the trapped wave to have a spatially varying phase relationship between pressure and vertical velocity.

The *p*′*w*′ quadrature spectra from RF8/IOP 10 (Fig. 13b) is qualitatively similar to the *w* power spectra. At an altitude of *z* = 12.2 km and below, it consistently shows significant negative quadrature, whereas above *z* = 13.1 km the quadrature is consistently positive. Figure 14 shows an example of the positive and negative quadrature signature. Statistics are given in Table 2 (see also section 4c and appendix C). The observed quadrature indicate that the TIL trapped wave node is somewhere between *z* = 12.2 km and *z* = 13.1 km. Referring to the IOP 10 balloon sounding in Fig. 11 and the Scorer parameter reference line corresponding to 12 km, there is a strong indication of a wave duct capable of supporting a horizontal wavelength of 13 km centered at approximately *z* = 12.5 km.

### d. Reflected or distorted long wave in IOP 13 (ferry)

A case of a partially reflected mountain wave was seen in IOP 13 (RF10 ferry) over the White Mountains (Fig. 15). Unfortunately, the leg length was quite short (∼85 km), causing COI problems. Ignoring the pulse inside the COI with wavelength of 40 km, we see only one pulse of mountain waves with a wavelength of 20 km. Comparing the wavelet diagram with the spatial plot shows that the wavelet analysis has broadened the pulse considerably. While the broadening has pushed the edges into the COI, we have confidence in the center part of the pulse. It has the positive cospectrum of an up-propagating wave. The quadrature spectrum is also significant (positive) at the same wavelength, suggesting partial reflection. Alternatively, the quadrature spectrum could be caused by a distortion of the wave structure due to nonlinearity in the wave. The large amplitude of the wave (i.e., +6 and −5 m s^{−1}) supports this idea.

### e. Up and down waves plus turbulence in IOP 6

In IOP 6/RF5 a long ferry leg found a mountain wave over the Wasatch Range in Utah in Fig. 16. Because of the long leg length (900 km) we avoid the COI problem, but we also increase the total variance in the series. This can have the unintended consequence of causing real wave signals to not show as statistically significant. We focus on the strong wave near *x* = −130 km. The *w*′ power wavelet diagram shows three discrete pulses at wavelengths of 15, 23, and 40 km. The longer two (23 and 40 km) have narrow spectra while the short one (15 km) is broad. The cospectrum diagram shows that the longest wave (40 km) is up-propagating while the 23-km wave is down-propagating. The magnitude of this down-propagating energy is nearly 3 times the magnitude of the up-propagating wave. The narrow 15-km pulse has little or no cospectrum or quad spectrum. We interpret the lack of phase relationship between pressure and velocity as evidence of turbulence. Its coincidence with wave propagation suggests turbulent gravity wave breaking.

## 7. Conclusions regarding wavelet analysis of mountain waves

Using linear theory and aircraft data, we have evaluated the utility of energy flux wavelet analysis in mountain waves. As in other wavelet applications, the primary advantage is the simultaneous presentation of spectral and spatial information. The Queney hydrostatic mountain-wave solution illustrates another important characteristic of wavelet analysis. While the wave phase lines tilt upstream with height, they do not disperse, and the wavelet power and cross-spectra diagrams are the same at all altitudes.

If the wave disperses aloft because of nonhydrostatic dynamics, the spectral content is unchanged but a chirp develops on the wavelet diagram, with shorter wavelengths trailing behind. Thus wavelets reveal the important wave dispersion characteristics while ignoring the less important detail of the wave pattern changing with height.

In this work, we emphasize the use of complex wavelets and cross-spectra to represent the phase relationship between pressure and velocity. With wavelet cross-spectra the wave location, wavelength, and propagation direction can be simultaneously represented. Complicated multiscale wave fields require such tools to reveal their true nature. As given in Table 3, each wave phenomenon has its own wavelet energy cross-spectrum characteristics.

Wavelet analysis of mountain waves requires some care and judgment. Tests for statistical significance allow differentiation of gravity waves from background noise. In legs with multiple wave regions, care must be taken that large-amplitude waves do not make smaller-amplitude gravity waves fail the significance test. Boundary influences were found to be an important limiting factor that prevents detailed analysis of mountain waves present near the edge of a dataset. Long leg lengths are recommended. Signals limited to the cone of influence must be ignored. Two wavelet artifacts must be recognized: the broadening of the wave field in physical space due to our *ω*_{0} = 6 wiggly Morlet wavelet and the weak false quadrature in regions with strong spatial gradient.

In principle, wavelet *u*′*w*′ cross-spectra analysis could be carried out for momentum flux estimation as well. In steady waves, the linearized Bernoulli equation relates *u*′ and *p*′ so that the two methods should be equivalent. In practice, however, the perturbation longitudinal velocity is noisier than pressure, probably due to the vertical displacement of pre-existing shear layers in the environment (S08). In addition, in unsteady waves, the energy flux is more useful than the momentum flux for indentifying wave source regions.

## 8. Conclusions regarding mountain waves near the tropopause

In this paper, we have shown that the secondary wave event on 17 April 2006 reported by S08 was not a unique occurrence. Although our database is limited (i.e., five cases with strong mountain waves) we have herein documented several more occurrences (Table 3). In our selection of cases, we have tried not to bias any particular pattern; however, research flights were operated only when large-amplitude wave activity was forecast. Surprisingly, no cases of pure vertical propagation were seen. The strongest wave case (RF10 ferry) shows upward wave energy flux with a moderate quad spectrum. The quad spectrum could be caused by partial reflection or a nonlinear distortion of the wave structure.

The most common type of wave pattern in the T-REX dataset is the collocated long up-propagating and short down-propagating wave pair (Tables 3 and 4). This was seen in RF4 (not shown), RF5, RF5-ferry, and RF10. Typically, the up-propagating long wave has a wavelength of 30–40 km while the down-propagating wave has a wavelength of 15–20 km. The RF10 case was discussed earlier by S08. It is unique in that the energy flux in the down-propagating wave exceeds the up-propagating wave at some levels.

An important property of the down-propagating waves is their horizontal phase speed (*C*). Using the aircraft data, there are three potential methods for estimating *C*. First, if the wave crests and troughs are in the same location during each pass, we can conclude that the waves are stationary and *C* = 0. In practice, as the short down-propagating waves are usually mixed in with other waves, such a determination is difficult (but see Fig. 8). Second, the relationship between vertical wind speed and potential temperature [i.e., *ik*(*U* − *C*)(*θ̂*) + *ŵ*(*dθ*/*dz*) = 0] for any Fourier component could be used to find *C*, but the local lapse rate (*dθ*/*dz*) is poorly known. The third method uses the linearized Bernoulli equation, written in coordinates moving with the wave [i.e., *ρ*_{0} (*U* − *C*)*u*′ = −*p*′]. With this generalization, the EP theorem can be rewritten and solved for the speed at which air flows through the wave [i.e., (*U* − *C*) = −EF/**MF**]. Values of (*U* − *C*) are given in Table 4 for a number of down-propagating waves. In general, the up-propagating waves have (*U* − *C*) similar to the wind speed *U*, suggesting stationary waves (i.e., *C* = 0). The down-propagating wave speed values have a bit more scatter, but still no systematic evidence for a nonzero phase speed. We conclude that the down-propagating waves are probably stationary. This conclusion may not be compatible with a Kelvin–Helmholtz mechanism for secondary wave generation.

Two examples of trapped waves on the tropopause inversion layer are reported here in cases RF8 and RF10 (Table 3). The wavelength in both cases was about 15 km. In RF8 aircraft sampling straddled the wave node, while in RF10 we sampled only below the wave node. From the Scorer parameter profile (Fig. 11), it is very unlikely that this short wave energy could have reached the tropopause from a ground source. We suggest instead that the wave was generated in situ by a nonlinear steepening process. In both cases the Scorer parameter profile indicates a lower reflector, consistent with Lindzen and Tung (1976) and Fritts and Yuan (1989). However, an upper reflector is only indicated in RF8; the sounding is insufficiently deep to indicate such a reflector in RF10.

One case of turbulent wave breaking was found (RF5 ferry). This was also a case of collocated up- and down-propagating waves. The turbulent region was spatially compact (∼20 km) but with a broad spectrum (wavelengths from 8 to 18 km).

In future work, mountain waves in different environments should be documented with wavelet energy flux methods. Gravity waves from transient sources such as convective clouds are also suitable for wavelet energy flux analysis. At the present time it is not clear whether these newly discovered “secondary” short waves in the stratosphere are unique to the Sierra region or whether they were missed or ignored in earlier projects elsewhere. Secondary wave generation, including down-propagating and tropopause trapped waves, was nearly ubiquitous in the high-amplitude wave cases from T-REX. Unfortunately, the physics of secondary wave generation is not well documented during the T-REX cases as the aircraft did not seem to penetrate the “active” nonlinear regions. Instead, the aircraft documented just the propagating secondary waves, with the possible exceptions of the RF5 ferry and *z* = 13 km legs of RF10.

Secondary gravity wave generation has been numerically modeled by numerous authors. Some studies used two-dimensional models (e.g., Franke and Robinson 1999; Satomura and Sato 1999; Lane and Sharman 2006) while others argue that three-dimensional simulations are needed to accurately capture the full physics (Andreassen et al. 1994; Fritts and Alexander 2003). The secondary wave generation we observed appears on both transects in the T-REX dataset and is consistent with both hypotheses. Some models predict long wave generation (Vadas et al. 2003); others predict short wave generation (McHugh 2009). The short wave harmonics generated by McHugh resemble neither the short down-propagating waves we find above the tropopause nor the short waves trapped on the inversion.

The down-propagating waves we observed carrying significant energy fluxes are consistent with simulations showing secondary waves contributing significantly the total wave momentum fluxes (Bacmeister and Schoeberl 1989; Vadas et al. 2003). While the waves we observed were forced by topography, secondary generation is also possible from convection (Holton and Alexander 1999; Lane and Sharman 2006). The generation mechanism for the secondary waves that we observed is unknown. Modeling studies have found wave breaking due to convective instability (Walterscheid and Schubert 1990; Andreassen et al. 1994), shear instability (Mahalov and Moustaoui 2009), and with contributions from both convective and shear instabilities (Liu et al. 1999; Satomura and Sato 1999). Even in cases for which we believe the aircraft sampled in situ wave breaking, the horizontal transects do not provide any vertical information of either the buoyancy or shear forcing terms to support either mechanism. More work will be required to reveal the generation mechanism for our T-REX short waves.

## Acknowledgments

Wavelet software was provided by C. Torrence and G. P. Compo and is available online (http://paos.colorado.edu/research/wavelets/).

Analysis was performed in NCAR Command Language (NCL) and scripts were modified from functions included in NCL. The authors thank the entire NCL team, especially Mary Haley and Dennis Shea. T-REX data were collected with the help of the NCAR flight facility staff and the whole T-REX staff. This research was supported by a grant to Yale University from the National Science Foundation (ATM-112354).

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Wavelet Cross-Spectra Confidence Limits

The confidence interval for wavelet cross-spectra was computed by

where *χ _{υ}*

^{2}(

*p*/2) is the chi-square distribution for

*υ*degrees of freedom (

*υ*= 1 for real wavelets,

*υ*= 2 for complex wavelets),

*k*= 0, … ,

*N*/2 is the frequency index,

*P*is the theoretical Markov red-noise spectrum for series

_{k}^{X}*X*at

*k*, and

*p*is the significance (TC98).

### APPENDIX B

#### Wavelet “False” Quadrature

When using pressure–velocity wavelet cross-spectra to identify up-propagating, down-propagating, and trapped waves, it is important to verify that the wavelet transform retains the phase relationships in the original signal. A counterexample was found in the Queney hydrostatic wave in section 3. In this case, every Fourier component comprising Eqs. (14)–(16) is upward propagating and has a perfect correlation (i.e., 0 phase lag) between pressure and velocity. The wavelet transform of Eqs. (14)–(16), however, showed a strong cospectrum but also had a small quad spectrum (Fig. 3b). We refer to this as a false quad spectrum.

Our definition of a false quadrature signal is as follows. Consider any two functions *F*(*x*) and *G*(*x*) whose Fourier transforms are parallel (in the complex plane) for all wavenumbers, for example, Arg [*F̂*(*k*)] = Arg [*Ĝ*(*k*)] or [*F̂*(*k*)*Ĝ**(*k*)] = 0. The wavelet quad spectra for these two functions are defined as false (i.e., they are artifacts of the wavelet method). A simple example of false quadrature is developed by choosing *F*(*x*) and *G*(*x*) to be similar symmetric functions of different widths; e.g., two cocentered Gaussians. Their Fourier transforms are parallel. Their Fourier quadrature is zero. Their wavelet quadrature is nonzero. In a similar fashion, a false cospectrum could be defined from a pair of functions whose Fourier transforms are orthogonal.

Trials using the Queney solution reveal that the magnitude of the false quadrature spectrum is related to the spatial gradient of the cospectrum and the scale according to

where Qu is the quadrature spectrum, Co is the cospectrum, *A* is a nondimensional coefficient, and *λ* is the horizontal wave scale (Fig. 3b). In addition to the Queney solution, we examined linear model solutions at different altitudes with Gaussian and Witch of Agnesi mountain shapes and tested half-widths ranging from 1 to 15 km. A few numerical mountain-wave simulations using the idealized 2D hill case in the Advanced Research WRF model were also consistent with Eq. (B1). A constant *A* ≅ 0.011 ± 0.001 was found for all these trials for the scales of interest for mountain waves, although there was some deviation at shorter scales where little energy was present. Other arbitrary function pairs will not satisfy Eq. (B1), but the ratio of quad to cospectrum will be small if the spatial gradients are small.

Two conclusions follow from Eq. (B1). First, if the wave is slowly varying in *x*, the false quadrature will be small, in part because of the small derivative and in part because the constant *A* is small. Second, if one integrates Eq. (B1) with respect to *x* over an isolated wave field, the false quadrature vanishes, as would a Fourier quadrature of perfectly correlated pressure-velocity signals. The false quadrature is probably too small to significantly degrade real cases, and it was not detected in our analysis of observational data.

### APPENDIX C

#### Phase Reversal across a Trapped Wave Node

A 2D trapped mountain gravity wave consists of a steady pattern of vertical and horizontal velocity perturbation *w*′(*u*, *z*), *u*′(*u*, *z*) along with pressure perturbations *p*′(*x*, *z*). Fourier transforming in *x*, we write each perturbation in the form *w*(*x*, *z*) = *ŵ*(*z*, *k*)exp(*ikx*). The continuity equation *ikû* + *dŵ*/*dz* = 0 and the steady Bernoulli equation with no background shear *ρUû* = −*p̂* can be combined to give a relationship between perturbation pressure and vertical velocity:

With no vertical propagation, the phase lines are vertical so the ratio *r*(*z*) = (1/*ŵ*)(*dŵ*/*dz*) is real. From Eq. (C1), pressure and vertical velocity are in quadrature. All trapped waves have one or more nodes where the vertical velocity is a maximum and *p̂* = *û* = *dŵ*/*dz* = 0. These nodes typically occur at or near levels of large Scorer parameter (e.g., stable layers). Above the wave node *r*(*z*) < 0 and (assuming *U* and *k* are positive), *p̂* ∼ *iŵ*; accordingly, moving downwind, pressure leads vertical velocity, leading to positive quadrature. Below the wave node *r*(*z*) > 0 and *p̂* ∼ −*iŵ*; thus, moving downwind, pressure lags vertical velocity leading to negative quadrature.

The physics of this phase shift across the wave node is illustrated in Fig. C1. Above the wave node, the streamline amplitude decreases aloft. The flow speed is larger, and the pressure lower, where streamlines pinch together. Positive displacement coincides with faster wind and low pressure. Below the wave node, the streamline amplitude increases with height. Positive displacement coincides with slower wind and higher pressure. The phase relation between vertical velocity and pressure reverses across the node.

## Footnotes

*Corresponding author address:* Bryan K. Woods, Yale University, KGL 106D, P.O. 208109, New Haven, CT 06520–8109. Email: bryan.woods@yale.edu

This article included in the Terrain-Induced Rotor Experiment (T-Rex) special collection.