Exploring gravity wave characteristics in 3-D using a novel S-transform technique: AIRS/Aqua measurements over the Southern Andes and Drake Passage

Gravity waves (GWs) transport momentum and energy in the atmosphere, exerting a profound influence on the global circulation. Accurately measuring them is thus vital both for understanding the atmosphere and for developing the next generation of weather forecasting and climate prediction models. However, it has proven very difficult to measure the full set of GW parameters from satellite measurements, which are the only suitable observations with global coverage. This is particularly critical at latitudes close to 60 S, where climate models significantly under-represent wave momentum fluxes. Here, we present a novel fully 3-D method for detecting and characterising GWs in the stratosphere. This method is based around a 3-D Stockwell transform, and can be applied retrospectively to existing observed data. This is the first scientific use of this spectral analysis technique. We apply our method to high-resolution 3D atmospheric temperature data from AIRS/Aqua over the altitude range 20–60 km. Our method allows us to determine a wide range of parameters for each wave detected. These include amplitude, propagation direction, horizontal/vertical wavelength, height/direction-resolved momentum fluxes (MFs), and phase and group velocity vectors. The latter three have not previously been measured from an individual satellite instrument. We demonstrate this method over the region around the Southern Andes and Antarctic Peninsula, the largest known sources of GW MFs near the 60 S belt. Our analyses reveal the presence of strongly intermittent highly directionally focused GWs with very high momentum fluxes (∼ 80–100 mPa or more at 30 km altitude). These waves are closely associated with the mountains rather than the open ocean of the Drake Passage. Measured fluxes are directed orthogonal to both mountain ranges, consistent with an orographic source mechanism, and are largest in winter. Further, our measurements of wave group velocity vectors show clear observational evidence that these waves are strongly focused into the polar night wind jet, and thus may contribute significantly to the “missing momentum” at these latitudes. These results demonstrate the capabilities of our new method, which provides a powerful tool for delivering the observations required for the next generation of weather and climate models.


Introduction
Atmospheric gravity waves (GWs) are propagating disturbances which transport energy and momentum, coupling and connecting the atmospheric layers.They are generated by processes including wind flow over mountains, jet stream instabilities, and meteorological sources including weather systems and convection (e.g.Fritts and Alexander, 2003;Alexander et al., 2010).
GWs are ubiquitous in and vital to the Earth's atmosphere.They contribute significantly to driving the strato-Published by Copernicus Publications on behalf of the European Geosciences Union.spheric winds (Alexander and Rosenlof, 1996), the mesospheric Equator-to-pole and pole-to-pole circulations (Andrews et al., 1987), and dynamical processes as diverse as the tropical quasi-biennial oscillation (Baldwin et al., 2001) and by preconditioning the polar vortex for sudden stratospheric warmings (Wright et al., 2010).They affect ozone depletion (Carslaw et al., 1998) and stratospheric water vapour (Kim and Alexander, 2013), and cause clear air turbulence affecting aircraft (Williams and Joshi, 2013).They also couple into the electrically charged ionosphere (Hooke, 1968), where they can generate disturbances which disrupt communication and navigation systems (MacDougall et al., 2009).
Characterising their distribution and behaviour is thus vital to understanding the atmospheric system and to advancing weather and climate modelling, for both weather forecasting and longer-term climate change predictions.In particular, models are believed to significantly under-represent waves fluxes near the 60 • S latitude band (McLandress et al., 2012), which is in turn believed to contribute to the "cold pole problem", arguably the most significant bias in existing models (Butchart et al., 2011).
GWs can be observed using a wide range of techniques (e.g.Fritts and Alexander, 2003;Alexander et al., 2010).These include ground-based methods such as radars, lidars and imagers and in situ methods such as tracer balloons and radiosondes.However, such techniques are intrinsically restricted to specific locations, leading to biases in our measurements, and to limited coverage over remote regions such as the open ocean and deserts.Only satellites have the geographic coverage and consistency needed to constrain GW effects in global models.
Since the mid-1990s, advances in satellite technology and methods have let us begin to characterise the global GW distribution (e.g.Fetzer and Gille, 1994;Wu and Waters, 1996;Preusse et al., 1999;Eckermann and Preusse, 1999;Hoffmann and Alexander, 2009).In particular, more recent methods have been able to provide first-order estimates of the absolute (i.e.non-directional) momentum flux (MF) transported by GWs (e.g.Ern et al., 2004;Alexander et al., 2008;Wright and Gille, 2013).This MF, which transfers to the background wind when the wave breaks, is arguably the main way GWs act upon the dynamics of the atmosphere, and is thus of fundamental geophysical importance.
However, until very recently, satellite GW measurements have been restricted to 1-D and 2-D by the scanning patterns that the instruments use.This is a critical issue, because the high travel speed of satellites in low Earth orbit prevents us from monitoring the time-evolution of observed GWs.As a result, the only way to measure the directions that GWs travel is by measuring their 3-D structure and hence their wavenumber vector.Directional information is key to understanding how the MF associated with GWs acts upon the atmosphere -without it, for example, we usually cannot even determine if a particular wave will accelerate or decelerate the local winds when it breaks.Three-dimensional measure-ments can also give us access to many other GW properties, including their phase speeds, group velocities and frequencies.These properties are fundamental to how models simulate waves, and are extremely poorly observationally constrained at the global scale.
Since 2015, two new techniques (Alexander, 2015;Schmidt et al., 2016;Wright et al., 2016a) have been proposed partially addressing this problem.These techniques overlap 1-D and 2-D observations from multiple satellite instruments to produce pseudo-3-D observations, allowing the measurement of 3-D wave structure in restricted planes.This has allowed preliminary climatologies of GW direction to be produced.However, these methods are heavily restricted in spatiotemporal coverage, require major assumptions to be made about the intercomparability of the datasets used, and are either vulnerable to significant aliasing issues (Alexander, 2015, due to the issues described by, for example, Faber et al. (2013)) or very significantly underestimate wave amplitudes relative to other observations and to theory (Wright et al., 2016a).
Since these studies, computational developments have finally made available high-resolution 3-D satellite measurements of atmospheric temperature.Hoffmann and Alexander (2009) demonstrated that radiance measurements from the nadir-sounding Atmospheric Infrared Sounder (AIRS) on NASA's Aqua satellite could be used to produce 3-D atmospheric temperature measurements capable of resolving GWs.However, the high computational cost of their analysis (∼ 2000-4000 core hours per retrieved day on 2008 hardware) made retrieval of the extended measurement periods necessary for climatological GW studies prohibitive at the time.With improved hardware performance per unit price, this level of computational investment has become feasible for large volumes of data, making available 3-D AIRS temperature estimates with global stratospheric coverage.Ern et al. (2017) were the first to exploit these data for 3-D GW measurements.They used the existing S3-D wavecharacterisation method (Lehmann et al., 2012) to produce a global directional GW MF climatology.The S3-D analysis method is relatively computationally cheap and has been previously tested extensively on model data (e.g.Preusse et al., 2014).However, it does not fully exploit the capabilities of the AIRS data.In particular, (i) the spatial scales of measured waves are dependent on the spatial scales of the analysis cubes which must be defined before the analysis, and (ii) the method assumes spatial homogeneity of the wave field over these cubes, which becomes increasingly unlikely as cube size increases.
Here, we develop and apply a new 3-D spectral analysis method, exploiting the AIRS data to a much fuller extent.Our technique is based around a generalised multidimensional extension of the existing 1-D and 2-D Stransforms (STs; Stockwell et al., 1996;Hindley et al., 2016).The use of an ST allows us to systematically analyse the data across the full range of length scales within the AIRS-resolvable wavelength spectrum, while simultaneously estimating wave amplitudes as they vary across individual wavepackets.Our method is capable of automatically detecting and characterising GWs and their associated properties, including their amplitudes, frequencies, directional MFs, and phase and group velocities.Using the single assumption, in common with previous studies, that the waves resolved by our data propagate upwards (Wright et al., 2016a;Ern et al., 2017), we are then able to unambiguously identify the propagation direction of all vector quantities to within measurement uncertainties.
Our measured properties can be geolocated to the voxel level, allowing us to study spatial resolutions finer than ∼ 1 • of latitude or longitude and allowing us to study individual GW packets in 3-D.We believe our results to be both the first height-resolved directional GW MFs measured by a single satellite, and also the first satellite measurements of 3-D GW phase and group velocity.
Our analysis method has several other benefits.It allows for extension to higher dimensions, for example in the analysis of time-varying model or ground-based datasets.It also includes a correction ameliorating the amplitude loss associated with standard ST methods.Finally, although we do not do so due to computational constraints, our implementation can be easily extended to measure an arbitrary number of multiple overlapping wave modes in an analogous way to the 2-D study of Wright and Gille (2013).
Section 2 describes the geography of the Andes/Drake Passage region, where we test our analysis.Section 3 then describes the 3-D AIRS temperature data we use.We next introduce our 3-D ST analysis (Sect.4), and discuss some implementation choices that may affect our results.We then describe a case study of a single pass of AIRS data known from previous work to contain GW signatures in Sect.5, then extend our analysis to produce a regional climatology of two winters (2014 and 2015) in Sect.6, including estimates of MF, wave direction and short-timescale wave variability.Section 7 then uses our new data to demonstrate observationally that GWs generated over the Andes are focused into the zonal wind jet in this region, an example of a result which could not be achieved with previous satellite observations.Finally, we draw brief conclusions in Sect.8.

Gravity waves and the Drake Passage region
To demonstrate our method, we study the geographic region around the Drake Passage, 80-50 • W, 70-35 • S.This region is a major wave source region, and GWs have been observed here using a vast range of methods including ground-based radar (e.g.Fritts et al., 2010Fritts et al., , 2012;;Wright et al., 2016b), in situ super-pressure balloons (e.g.Hertzog et al., 2008;Plougonven et al., 2013), and satellite limb (e.g.Eckermann and Preusse, 1999;Jiang, 2002;Alexander and Teitelbaum, 2011) and nadir sounders (e.g.Alexander and Barnet, 2007;Hindley et al., 2016;Hoffmann et al., 2016).These results strongly suggest that the region is the most significant GW source worldwide, with peak GW amplitudes (Yan et al., 2010) and MFs (Geller et al., 2013) at least a factor of 2 and potentially an order of magnitude greater than any other known source.
We consider here the consecutive autumns and winters (May to September) of 2014 and 2015.This portion of the year is the main period for GW activity in this region, with previous studies using a range of instruments measuring negligible GW activity during the rest of the year above 30 km altitude (e.g.Wright et al., 2016b).
Figure 1 shows the geography of the region.The area has a mix of flat open ocean and a variety of land regions, including both the high orography of major mountain ranges and broad lowlands such as the Pampas.Two mountain ranges dominate the region topographically: (1) the Antarctic Peninsula mountains, which lie along a NE-SW axis in the southern part of the region and (2) the Southern Andes, which form a N-S barrier along the western edge of mainland South America shifting to a NW-SE barrier along the southern coast of Tierra del Fuego.This southern extension of the Andes, while much lower in altitude than the main part, can still form a significant barrier to surface winds.
GWs in the region are believed to be primarily orographic in source, arising from fast zonal surface winds flowing across mountain ranges lying near-perpendicular to the prevailing wind.Once generated, the waves can propagate downstream across the Southern Ocean, where they contribute strongly to a leeward "tail" of GW activity (e.g.Sato et al., 2012;Hindley et al., 2015).This tail is believed to contribute very large quantities of MF to the climate system, and is thus of great geophysical significance.

Data
We use data from the Atmospheric InfraRed Sounder (AIRS) on NASA's Aqua satellite (Aumann et al., 2003;Chahine et al., 2006;Olsen et al., 2007).Launched on 4 May 2002, Aqua flies in a Sun-synchronous near-polar orbit as part of the A-Train satellite constellation, with an ascending-node Equator-crossing local solar time of 13:30.Aqua completes 14.55 orbits per day, with a 16-day repeat cycle.
AIRS scans across-track from +49.5 to −49.5 • offnadir axis, measuring a continuous swath of radiances over 2378 spectral channels.The data are processed as 90 parallel scan tracks, with horizontal resolution varying from ∼ 13.5 km × 13.5 km at nadir to 41 km × 21.4 km at track edge.
The scan track is split into arbitrary 135 along-track element sections, referred to as granules, which provide the basic unit of data storage.These correspond to 6 min of data collection.Approximately 650 AIRS measurement granules overlap our analysis region each month, providing sufficient data for statistical analysis of our results.
We specifically analyse 3-D stratospheric temperatures derived from AIRS radiance measurements in the 4.3 and 15 µm CO 2 wavebands, using the retrieval scheme developed by Hoffmann and Alexander (2009).Retrievals are carried out for each satellite footprint independently, at the full horizontal sampling capacity of the instrument.The horizontal resolution of the retrieval is enhanced by a factor of 3 × 3 in the along-and across-track directions compared with the AIRS operational data, which is an important asset for GW analyses.The total retrieval error of the individual temperature measurements varies between 1.6 and 3.0 K in the altitude range from 20 to 60 km.Retrieval noise is the leading error and varies between 1.4 and 2.1 K in the same vertical range.The vertical resolution of the temperature retrievals is 7 to 15 km.Further validation of the retrievals is discussed by Hoffmann and Alexander (2009) and Meyer and Hoffmann (2014).We note that the retrieval method uses different channels for daytime and nighttime; we combine both periods here.
Temperature fluctuations due to GWs are derived from the 3-D AIRS dataset by subtracting a 4th-order polynomial fit for each across-track scan (Wu, 2004;Alexander and Barnet, 2007;Hoffmann et al., 2014).This detrending method effectively removes slowly varying background signals due to large-scale temperature gradients or planetary wave activity.The AIRS sensitivity function to GWs is defined by both the averaging kernels of the temperature retrieval and the detrending method.The dataset discussed in our study is sensitive to GWs with vertical wavelengths larger than about 15 km.The sensitivity cutoff at short horizontal wavelengths depends on the footprint size and varies between 30 km at nadir and 80 km for the outermost tracks.For long horizontal wavelengths, sensitivity drops below 90 % at horizontal wavelengths of 730 km and below 10 % at 1400 km; thus, longer horizontal wavelengths will be more strongly attenuated than short ones, modifying the observed spectrum.For a more detailed discussion see Hoffmann et al. (2014) and Ern et al. (2017).

The 3-D S-transform
The ST (Stockwell et al., 1996) is a widely used tool for localised time-frequency (or distance-wavenumber) data analysis.In particular, the 1-D ST has been extensively used in recent years for GW detection and characterisation in atmospheric datasets (e.g.Alexander et al., 2008;Fritts et al., 2010;Hertzog et al., 2012;Wright and Gille, 2013).More recently, the 2-D extension of the ST has been similarly applied to atmospheric datasets including single-altitude stratospheric AIRS retrievals (Hindley et al., 2016) and mesospheric airglow measurements (Stockwell et al., 2011).Here, we introduce a new three-dimensional implementation of the ST for full 3-D GW analysis of our AIRS measurements, which we generalise to N dimensions for future use in other contexts.
For a smoothly varying, continuous and one-dimensional function of time h(t), the 1-D ST S(τ, f ) is defined as where τ represents translation in the time (or spatial) domain and f frequency.Here, c is a scaling parameter which is usually set to 1, but which can be adjusted to improve temporal localisation at the expense of frequency localisation or vice versa (Mansinha et al., 1997;Fritts et al., 1998;Pinnegar and Mansinha, 2003;Hindley et al., 2016).Equation (1) uses as an apodising function a Gaussian window whose standard deviation is scaled as the inverse of frequency, although we note for completeness that any suitable such function may be used as long as it has a spatial integral equal to unity (Hindley et al., 2016).
The ST in Eq. ( 1) can easily be extended to higher dimensions.However, such an extension becomes quite unwieldy beyond the 2-D case.We therefore use a more compact vector notation to define an N-dimensional ST.Specifically, for any function h(x), where x = (x 1 , x 2 , . .., x N ) is a column vector describing an N-dimensional coordinate system, we can write the ST S(τ , f ) as Here, τ = (τ 1 , τ 2 , . .., τ N ) and f = (f 1 , f 2 , . .., f N ) are column vectors denoting translations and spatial frequencies (inverse of wavelength) in the x 1 , x 2 , . .., x N directions, and f denotes the transpose of f .The scaling parameter c n is a scalar quantity in each dimension n that can be adjusted in order to spatially/spectrally tune the localisation capabilities of the N-dimensional ST for each dimension independently (Fritts et al., 1998;Hindley et al., 2016), and hence represents a tunable parameter set which could be used in future work to emphasise different wave properties.
In the 3-D AIRS analyses presented here, we have empirically selected c n = (c x , c y , c z ) = (0.25, 0.25, 0.1); this produces physically plausible results in tests, and is justifiable on the basis that GWs at the length scales visible to AIRS typically occupy a much larger fraction of the granule in the vertical than in the horizontal, allowing us to boost spectral localisation in this direction at a cost of less-critical spatial localisation.

Methodological outline
Our 3-D ST-based GW analysis involves seven key steps.These steps are outlined in Sects.4.2-4.6.We then discuss some possible methodological restrictions and how we expect them to affect our results in Sect.4.7.
1. Firstly, we detrend the raw temperature using a fourthorder cross-track polynomial, as described in Sect. 3 above.This leaves us a three-dimensional temperature perturbation field T (x, y, z).
2. We next interpolate T (x, y, z) onto a regular spatial grid.While in principle the resolution of this grid is arbitrary, our computational implementation of the 3-D ST makes heavy use of the fast Fourier transform, which requires regularly spaced data and which is most efficient if the number of elements in input dimensions are powers of 2. We thus interpolate the original 90 × 135 × 21 element data grid for each granule to a grid of 64 × 128 × 16 elements in the (cross-track, along-track, vertical) directions.
3. We then temporarily remove the exponential increase of wave amplitude with height, arising from the decline of atmospheric density with height.This was done to mitigate the effect of much larger wave amplitudes at the top of the granule unfairly dominating the S-transform spectrum via spectral leakage.We define a reference height level z 0 = 41 km, representing approximately the centre altitude of the original data, and apply a scaling factor exp (−(z − z 0 )/2H ) to all other height levels.Here, H is the atmospheric scale height, which we define as 7 km.
4. Next, we compute the 3-D ST.We then project the measured wavenumbers from the alongtrack/across-track/vertical frame of reference into the zonal/meridional/vertical wavenumbers k, l and m to provide a consistent geometric basis for our results.Section 4.3 below describes these steps.
5. The wave amplitudes estimated by the 3-D ST analysis exhibit significant amplitude reduction.This is due to fundamental limitations arising from the finite number of wave cycles present in any real wavepacket.While this problem is minor enough to be neglected for the 1-D ST (Wright, 2010), in the 3-D case it can reduce measured amplitudes by as much as 70 % of their true value.Accordingly, we next apply an empirically derived correction to restore the measured wave amplitudes to values close to their initial values in the AIRS retrieval.This is discussed in Sect.4.4 below.
6.We next reverse the amplitude scaling with height we applied above.This restores the true height-scaling of the measured wave amplitudes, typically exponentially increasing with height.
7. Finally, we use the wave amplitudes, wavenumbers and geolocations provided by the above steps to compute additional wave properties, including their phase speeds, group velocities and momentum fluxes.Sections 4.5-4.6 describe this step.

Computing the 3-D ST
For each granule of 3-D AIRS data analysed, we apply the vectorised ST S(τ , f ), Eq. ( 2), producing a six-dimensional complex-valued object, i.e. S (τ , f ) ≡ S(τ x , τ y , τ z , f x , f y , f z ).Following Hindley et al. (2016), we then collapse this 6-D object into a more computationally manageable 3-D object by considering only the spectral peak for each location in the localised ST spectrum.This peak corresponds physically to the single wave with the largest spectral amplitude at that location in (τ x , τ y , τ z ). Figure 1 of Hindley et al. (2016) visualises this for the 2-D case.
Once we have identified the spectral peak of the localised (f x , f y , f z ) spectrum for each location in (τ x , τ y , τ z ), we record the values of the complex coefficients at these spectral peaks in a 3-D object A(τ x , τ y , τ z ).The magnitude of the complex coefficients of A correspond to the underlying wave amplitude of the largest-amplitude (i.e.most dominant) wave at each location in (τ x , τ y , τ z ).As described by Hindley et al. (2016), taking the real part R(τ x , τ y , τ z ) = Re A(τ x , τ y , τ z ) also provides a useful "reconstruction" of the input granule using only the dominant wave at each location.This reconstruction can be used to visually assess the effectiveness of our 3-D ST analysis, and is also used for our amplitude correction in Sect.4.4 below.
The exact location of the spectral peak in each localised (f x , f y , f z ) spectrum is also recorded.This denotes the f x , f y and f z frequencies that correspond to the dominant wave at that location in (τ x , τ y , τ z ).We record each of these in the 3-D objects F x (τ x , τ y , τ z ), F x (τ x , τ y , τ z ) and F z (τ x , τ y , τ z ) respectively.
Conveniently, in our application the (τ x , τ y , τ z ) domain corresponds exactly to locations in the regular grid onto which we interpolated each granule of our input AIRS data, i.e. (τ x , τ y , τ z ) = (x, y, z).This effectively gives us A(x, y, z), F x (x, y, z), F y (x, y, z), R(x, y, z) and F z (x, y, z) relative to each granule.Thus, we measure and record the underlying amplitude and frequency of the dominant wave at every location in our interpolated 3-D granule of AIRS measurements.
We then follow the method of Hindley et al. (2016) to project the spatial wavenumbers f x ≡ F x (x, y, z), f y ≡ F y (x, y, z), f z ≡ F z (x, y, z) into zonal, meridonal and vertical wavenumbers k, l and m, by considering the geographic azimuths of the along-track and cross-track directions at each location on the granule.

3-D ST amplitude restoration
The N-dimensional ST analysis has many advantages.However, one major drawback arises which we must correct for.In the real atmosphere, GWs typically exist as wavepackets, i.e. a finite number of wavecycles contained within an amplitude envelope.Analytically, if we consider the wave inside the packet to be a perfect sinusoid with wavenumber k, the presence of a packet-like windowing function in the spatial domain effectively applies numerical prefactors to wavenumber voices in the spectral domain.This means that, when the ST is evaluated at the "correct" wavenumber k, the measured spectral amplitude at that wavenumber is reduced.Since the spectral power is not really lost, merely spread across other wavenumber voices and not attributed to the true wavenumber k, this effect is similar to spectral leakage.
For waves with a large number of wavecycles within the wavepacket, this effect is negligible, as shown by, for example, Wright (2010) for the 1-D ST.However, for higher dimensions and for waves with only a few wavecycles, this amplitude reduction can be significant.Stockwell et al. (2011) and Hindley et al. (2016) observed wave amplitudes reduced by 40-60 % of their initial values using the 2-D ST, and our investigations suggest 3-D ST-measured wave amplitudes can be reduced by as much as ∼ 70 % of their initial values.
Here we apply a simple correction to restore our 3-D STmeasured amplitude values to close to their initial values in the AIRS measurements.To do this, we make use of the initial interpolated 3-D AIRS measurements T (x, y, z) (after removal of their exponential amplitude trend with height) and the 3-D ST reconstruction R(x, y, z) of these measurements, which is found by taking the real part of the complex-valued object A(x, y, z) as described above.
These reconstructed wave amplitudes R(x, y, z) are reduced by a factor dependent on how many wave cycles were present in each wavepacket.Since complete knowledge of the wave field, specifically the number of wave cycles present in each wavepacket, is fundamentally very challenging, it is thus difficult to correct for each and every wave exactly and independently.
Our approach is to first take the 3-D Fourier transform of both the AIRS measurements and the reconstruction.We take the median of the linearised quotients of all the absolute values of the complex spectral coefficients, i.e.
This gives us a scalar value ξ for the general amplitude reduction factor of the dominant waves in the granule.We then multiply the complex-valued object A(x, y, z) by ξ in order to get a new, restored underlying wave amplitude and reconstruction for each location in the granule, by taking the magnitude and real parts of ξ A(x, y, z) respectively.
Although this is a relatively simplistic approach, it typically restores wave amplitudes of the dominant waves in the granule to within ∼ 10 % of their original values.It should be noted, however, that the use of a single scalar factor for the whole granule could introduce errors, where some less dominant wave features are rescaled by an incorrect amount.As a precaution, we only apply factors between 1 and 5. Future work will refine this correction by estimating the spatial extent of observed wavepackets in the granule, allowing us to measure the number of wave cycles present in the packet.This will allow us to better constrain the attenuation and separately apply the correct scaling factor for every wave in the granule independently.

Frequencies, phase speeds and group velocities
By simultaneously characterising k, l and m, our data allow us to estimate GW intrinsic (i.e.wind-relative) frequencies ω, phase speeds ĉp and group velocities ĉg .Further, by using ancillary background wind data, we can convert these properties to their ground-relative values ω, c g and c p .
Following Fritts and Alexander (2003), we calculate the GW intrinsic frequency ω from the GW dispersion relation where u and v are the background zonal and meridional wind speeds respectively, N B is the Brunt-Väisäla (buoyancy) frequency, H is the atmospheric scale height, (∼ 7 km for the stratosphere), and f = 2 sin(φ) is the Coriolis parameter for the Earth's rotation rate at a latitude of φ.This in turn allows us to compute the GW group velocity , where c gx , c gy , c gz are the components of group velocity in the zonal, meridional and vertical directions.
Finally, GW intrinsic phase speed ĉp is given by and is defined in the direction of phase propagation (k, l, m).

Momentum flux
We can also compute the MF associated with an observed wave.Under the midfrequency approximation, which is valid for our data, MF can be computed as (Ern et al., 2004) where M z and M m are the zonal and meridional projections respectively of the MF, ρ, and T are the local atmospheric density and background temperature, and g is the acceleration due to gravity.Equation ( 7) neglects several terms in the full treatment of Ern et al. (2004) involving quantities AIRS does not measure.Ern et al. (2017) show that, for these data, neglecting these terms typically leads to errors < 30 % provided the ratio m/k > 1.5.The range of observable wavelengths imposed by the resampling step in our analysis, described below, means that our measured m/k never falls below ∼ 2.5 in the worst case, and is typically much greater than this.
Note that our values of k and l are estimates of true values rather than lower bounds arising from 2-D projections of the wave field.Our estimate of MF is thus an estimate of the true value rather than a lower bound as in most previous satellite studies.

Assumption of upward propagation
Since the satellite measurements are instantaneous, there is a 180 • ambiguity in the wavevector (k, l, m) which cannot be directly broken from the observations.As in Wright et al. (2016a), we break this ambiguity by assuming that the GWs propagate vertically upwards, thus requiring that the vertical wavenumber m be negative.In order for the data to remain internally consistent, this implicitly constrains the other two components of the wavevector, providing a full-3-D wavevector at each point in the data.
Upwards propagation is very likely to be the case for the great majority of GWs in this region.In particular, the assumption is consistent with the model study of Sato et al. (2012), where GW energy flux across the entire region was upwards, with the exception of an area immediately downwind of Tierra del Fuego and below 30 km altitude.This assumption may lead to the misclassification of some GWs in our sample, especially any generated by high-altitude sources such as the edge of the polar vortex.

Handling of along-track granule-edge waves
From Sect. 6 onwards, we consider bulk averages of many AIRS granules.The arbitrary boundaries between these granules can cause GWs near the boundaries to be poorly resolved.For example, the GW studied by Fig. 6  To accommodate this issue, we process the data in two distinct passes, (i) as each granule individually and (ii) as the merged second half of each granule and the first half of the next, i.e. the original granule boundaries stepped along-track  by 1/2 granule.The results of both passes are then stored independently and treated as separate observations.This ensures that all GWs will be measured as accurately as the data permit in at least one of the two passes.Since every point in the input data is considered twice, such partially truncated waves may have a very small value in one pass and an accurate value in the other pass.Thus, the wave may contribute less to bulk averages than it would if fully resolved in both individual processing passes.This may affect our results slightly.

Computational limitations
Unlike the temperature retrieval, we implement our 3-D ST analysis on desktop hardware.Runtime in this context is reasonable (∼ 2 min per granule, equating to a total runtime ∼ 8 days for the period and region considered), but memory restrictions impose some limitations on our analysis.Specifically, we (1) use unsigned 8 bit integers, scaling all T values into this range; (2) interpolate the data to a spatial resolution of (128 × 64 × 16) points in the (along-track × across-track × vertical) direction, chosen as powers of 2 to optimise discrete FFT performance; and (3) do not consider spatial frequencies above P /15 in the horizontal, where P is the number of points in that direction after interpolation.
We do not expect these choices to affect our results significantly, for the following reasons: 1.An 8 bit integer allows 256 discrete values to be stored.
Once the background is removed, our range of retrieved temperatures is around ±40 K, though many granules do not nearly have as large a range as this, particularly as we remove the exponential increase in wave amplitude for this step.Converting to 8 bit integers means our precision is at worst reduced to around 0.3125 K for this range, whereas retrieval noise can vary from around 1.4 to 2-3 K (Sect.3).Thus, any introduced noise from converting to 8 bit integers does not significantly contribute to further uncertainty.
2. The spatial resolution values used are similar to those of the original data -128 vs. 135 points along-track, 64 vs. 90 points across-track, and 16 vs.13 points in the vertical.
3. The frequencies we discard correspond to signals spanning only low-single-digit numbers of points, and are thus potentially vulnerable to noise -in previous work using AIRS data (Wright et al., 2016a, b), we have presmoothed the input fields to remove these frequencies.
Nevertheless, it is important to note that these choices may affect our results, particularly across-track where we discard ∼ 1/3 of the input resolution.We therefore expect future work to be capable of better wave characterisation, even without changes in the analysis technique or the input data.In particular, our current analysis will discard very small GWs near the Nyquist limit of the original data, as observed by, for example, Alexander et al. (2009).
Furthermore, as described above and in common with Hindley et al. (2016), we discard all signatures in the data associated with frequencies other than the signature corresponding to the largest-amplitude wave at each spatial gridpoint.While this is necessary to implement the 3-D ST on our current hardware, we expect that future implementations will be able to use this discarded data to measure distinct overlapping wave signatures, as done in the 2-D case using HIRDLS data by Wright et al. (2015).Figure 2a shows a map of the region, with the AIRS scan track indicated by the red outline.AIRS is travelling from NE to SW, as indicated by the red arrows at track centre.Figure 2b-e then show T measured by AIRS at four height levels.We see distinct phase fronts of hot (red) and cold (blue) air tightly localised over Patagonia, changing with height in a manner consistent with vertically stacked phase fronts.The spatial distribution of this T is very similar to Fig. 1b of Wright et al. (2016a).This is unsurprising since the data share a common source, but it provides a useful cross-check on the accuracy of the 3-D temperature retrieval.T amplitudes are, however, much larger here, due to our use of retrieved kinetic temperature rather than brightness temperature.Temperature amplitudes are reduced at the highest level plotted (55 km), perhaps due to the reduced vertical resolution of the retrieval at these altitudes.
Figure 3 shows the same data, as a 3-D view from the northwest.This viewing angle has been chosen empirically as that from which the distinct GW phase fronts are visually clearest.T is plotted as 3-D semi-transparent nested surfaces of constant T .We see that the phase-front volumes are mostly continuous in 3-D.Assuming that the GW propagates from an orographic source in the Southern Andes, we also see the phase fronts tilt with height, shifting from a nearhorizontal alignment for the earliest phase fronts (at low altitude) of the wave to a near-vertical alignment for later ones (at high altitude).

Wave amplitudes and wavelengths
Figure 4a-e show the direct output of our 3-D ST analysis for this GW.We show the 30 km altitude level only, but similar diagnostics are computed by our analysis at all other heights within the measurement volume.
Figure 4a shows T and Fig. 4b T .Several wide phase fronts are visible over the Southern Cone, overlaid by a very fine set of narrow phase fronts centred at (50 • S, 74 • W).In what follows, we assume these two features represent two distinct GW signatures.We then define the large-amplitude region of the larger-area GW centred at (50 • S, 70 • W) as "Wave A" (outlined approximately by the thick dotted line) and the tightly defined narrower feature at (50 • S, 74 • W) as "Wave B" (outlined approximately by the thin solid line).We note, however, that other wave features exist in the same data, and that Wave A extends over most of the granule to some degree.
Wave B is very close to the data's horizontal resolution limit and could conceivably be spurious.However, the signal is visible across a range of heights, and analysis gives physically plausible results for most metrics (with the possible exception of phase speed and temporal frequency, for reasons discussed below).We thus believe it to be a real GW signal.GWs at similarly small scales in 2-D AIRS data have been studied by, for example, Alexander et al. (2009) and, while small in area, can dominate the MF distribution of a granule due to the strong inverse dependence of MF on λ h (see Eq. 7, below).This feature was invisible to Wright et al. (2016a) due both to the coarseness and relative spatial location of the MLS scan track, and is very close to the resolution limit of our 3-D AIRS data.
Figure 4c and d show λ h and wave propagation directions θ .λ h is consistent with the phase fronts seen in Fig. 4a.Propagation is mostly southwestwards.This is consistent with orographic waves generated by surface wind flow over the local NW-SE topographic orientation, which must be directed upwind in order to maintain a constant ground-based location.An exception is a region centred at (47 • S, 78 • W) where the waves propagate in a near-westward direction.This is due to the mountain ridge here being nearer to a N-S alignment.
Figure 4e shows measured λ z .Due to the small number of vertical levels used, the number of distinct λ z values measur- able is restricted (e.g.Wright et al., 2015), and thus we see a "speckled" pattern of only two values over the great majority of the region, λ z ∼ 14 km and λ z ∼ 18 km.For an orographic GW under the hydrostatic relation (e.g.Eckermann and Preusse, 1999), λ z ∼ 2π u h /N B , where u h is the background horizontal wind speed.Taking N B = 0.02 s −1 and u h = 50 ms −1 , the latter derived from ECMWF operational analyses for this day, this predicts λ z ∼ 16 km.Given our restricted height resolution prevents the measurement of λ z values between 14 and 18 km, this prediction is fully consistent with our observations across the region, and the true value almost certainly lies somewhere between these two.The values are also consistent with visual inspection of vertical sections through the region (Figs.5d, g, discussed below).

Frequencies, phase speeds and group velocities
To compute ω, c p and c g , we first pre-smooth the λ z field using a 3 × 3 voxel boxcar smoother in the horizontal plane only before calculating these properties.This is due to the "speckled" nature of the λ z distribution (Fig. 4e), which transfers directly through to all subsequent figures.Given we suspect the true value to lie between the values observed, this smoothing represents a better estimate of the true λ z value, and thus of derived quantities.No similar smoothing is applied to λ h .
Figure 4f-h show maps of intrinsic ω, ĉgh and ĉph , calculated after this smoothing.These values do not require us to use ancillary background wind data.They are shown for completeness only, and we do not discuss them further.Corresponding vertical speed maps are not shown since they are by definition identical to the ground-based maps on the bottom row of the figure.
Figure 4i-m then show ground-based properties.These are calculated using u and v values derived from ECMWF operational analyses.These estimates are typically insufficiently accurate for use at the level of localised wind perturbations u , v , but are well-suited to characterising the bulk behaviour of the local atmosphere u, v. Our ground-based values will, however, exhibit additional uncertainty due to the use of reanalysis data to derive them.In particular, in the region centred on 47 • S, 76 • W a large angle exists between the measured wave fronts and reanalysis winds, leading to an equivalent region of very large values in Figs.4j, l induced by the projection of the wave into the wind frame which may be unphysical.
Figure 4i shows measured ω. ω = 0.5×10 −3 s −1 for Wave A and ω = 2 × 10 −3 s −1 for Wave B. Both of these features are consistent with midfrequency GWs f = 2 ×10 −4 ω N B = 0.02 s −1 .Testing of a range of similar 3-D AIRS granules (not shown) indicates that this is usually the case for GWs detected in these data using this method.
Figure 4j and k show horizontal group speed c gh = c 2 gx + c 2 gy and vertical group speed c gz .Directionality is omitted in the horizontal plane for brevity, but is generally southeastwards across the region, consistent with group velocity perpendicular to the GW phase fronts.c gh ∼ 0-10 ms −1 for the majority of the region, again consistent with midfrequency waves (e.g.Wright et al., 2011).c gz is 1 ms −1 for Wave A and 3 ms −1 for Wave B. Values of c gh increase over the Drake Passage; this is consistent with focusing of the waves into the centre of the southern zonal wind jet in this region (Sato et al., 2009(Sato et al., , 2012;;Hindley et al., 2015), and is discussed further in Sect.7.
Finally, Fig. 4l and m show horizontal phase speed c ph = c 2 px + c 2 py and vertical phase speed c pz .Horizontal directions are primarily southwestwards (not shown).This is theoretically consistent with the observed group velocities and with the phase fronts depicted in Fig. 3f: the phase speed is directed across the phase fronts and downwards, while the group velocity is directed along the phase fronts and upwards.c ph ∼ 0 ms −1 for Wave B, consistent with orographic wave generation.Values vary across Wave A, but are typically low, consistent again with uncertainties in u. c pz is nearzero for Wave A, but ∼ 3 ms −1 for Wave B.
The relatively large values of ω and c pz for Wave B are somewhat surprising.Given the close spatial correspondence between this wave and the significant topography around Mt San Valentin, we would assume this wave to be orographic in source and thus to have near-zero values for both these quantities, as with Wave A. Accordingly, we suspect that this wave is indeed orographic in source, and thus that the nonzero magnitudes of these quantities arise from measurement uncertainties associated with our coarse-grained input ver-tical structure.While the same problem nominally exists for the larger Wave A, the much larger spatial extent of this wave provides more data to constrain our analysis, allowing more accurate determination of all quantities.

Momentum fluxes
Figures 5 and 6 shows the MF associated with our example wave, split into zonal M z and meridional M m components.Figure 5a-c show maps at 30 km altitude, Fig. 5d, e, f show height-longitude sections at 50 • S, and Fig. 5g-i show height-latitude sections at 70 • W. Black shading indicates topography for that section, scaled by 5 for ease of comparison.T , shown previously, is reproduced in Fig. 5a, d, g for ease of spatial comparison.Zonal MFs M z are shown in Fig. 5b, e, h and meridional MFs M m are shown in Fig. 5c, f, i. Figure 6a and b show the 3-D structure of M z and M m , from the same perspective as Fig. 2f.All plots are signed, i.e. show both positive and negative-valued data.Regions with T < 3 K in panel (a) and T < 3 K in panels (b) and (c) have been set to zero.
As with Fig. 4, the dominant features in Figs. 5 and 6 are Waves A and B. Wave A has a peak MF of 50 mPa, and Wave B 230 mPa (saturating the M z colour scale used in Fig. 5).Note that, while we believe our measurements of Wave B to be more uncertain than those of Wave A due to the limited vertical resolution of our measurements, even a halving of measured λ z would still cause the measured MF to saturate this colour scale.Wave B is thus clearly the largest MF signal in the region regardless of this uncertainty.Both waves extend in height over the majority of the stratosphere.
MF for Wave B is very tightly localised directly over the region around Mt San Valentin at (73 • S, 50 • W), with the flux forming a narrow column directly above the region.Wave A is much more spatially diffuse, leading to heightened MF values across the entire Southern Cone and out over the Drake Passage and Atlantic Ocean.Note, however, that the much larger volume of this wave allows the transport of more total MF than Wave B.
MF is seen to decline at both low and (particularly) high altitudes.While at high altitudes this is consistent with wave dissipation, we cannot decouple this effect easily from data limitations.At both the highest and lowest altitudes, (i) temperature retrieval accuracy is significantly reduced, (ii) the vertical resolution of the retrieval is also reduced, and (iii) long vertical wavelengths unavoidably experience some edge-truncation due to the limited vertical height window in our analysis (e.g.Wright et al., 2015), with effect (iii) being particularly important.In particular, small regions of northeastwards flux are seen at high and low altitudes at (50 • S, 75 • W); these may be due to (i) limited data quality (ii) our assumption of upward ascent being invalid or (iii) they may genuinely be GWs propagating in a northeastward direction.

Climatological MF
Having demonstrated our method on the above example, we now move on to produce a climatology of directional MF in the region.Only MF is considered for brevity; however, similar climatologies can be produced for all the variables discussed.For reference, amplitude, frequency, wavelength and group velocity distributions are typically smooth, while phase speed distributions exhibit significant noise.We assume this noise is due to the direct k −1 dependence of phase speed, which makes measurements of this extremely sensitive to any measurement accuracies.We see two dominant regions of high flux: the Southern Cone and the Antarctic Peninsula.Consistent with theory and with previous work, zonal MF in both of these regions is primarily negative (i.e.directed in an eastward direction), while meridional MF is negative (i.e.southward) over the Southern Cone and positive (i.e.northwards) over the Antarctic Peninsula.Measured MFs in both regions are largest in July, with the Southern Cone region showing activity across all five months considered and the peninsula primarily active in July and August only.In July we also see enhanced MF east (downstream) of the Southern Cone, consistent with GWs being advected downstream by the background winds.
In no month and at no height do we observe significant quantities of MF over the Drake Passage itself, with the exception of small regions close to coastlines.This observation is consistent with Wright et al. (2016a), who also observed high-horizontal resolution GWs in this region using a 3-D approach, and lends further credence to the hypothesis that the vast majority of MF in this region during autumn and winter is orographically generated rather than arising from non-orographic storm activity.However, we note that our observational filter preferentially selects for waves with long vertical wavelengths (>∼ 7 km), which are less likely to be observed over the Drake Passage if generated in the Andes due to their propagation characteristics.
Finally, Fig. 9 shows an example 3-D plot of |MF|, in this case for August 2014.Consistent with the maps above, we see two major |MF| peaks, one over the Southern Cone and one over the peninsula.Certain features are, however, clearer in this format.We see that the |MF| peak over the peninsula is extremely spatially well defined, forming a very narrow spike directly over the mountains.We also see that the en- hanced region over the Southern Cone has significant internal structure.

Time variation and direction
We now define three subregions for detailed analysis, indicated by coloured boxes on Fig. 1.Two of these subregions, the Southern Andes (purple) and the Antarctic Peninsula (or-ange), are chosen as the regions of highest GW activity over the period studied.The third, upstream in the open Southern Ocean (red), is a background region chosen for comparison.
Figure 10 shows time series of observed MF averaged over these three geographic regions.Data are shown as 7-day rolling means.Note that there is very significant variability at timescales shorter than this, which we quantify in Sect.6.3 below.
Figure 11, meanwhile, shows the direction of observed MF for our subregions as rose plots.Each row shows three height levels for each subregion, with the shaded boxes showing the proportion of measured MF in each direction.The data are binned into 30 • bins, with different altitude levels offset slightly to enhance legibility.Each height level distribution sums to 100 %, and each column shows a specific month, with the final column showing the combined results for all 5 months.Both 2014 and 2015 have been averaged together.
At the seasonal scale, the three regions exhibit very different MF magnitudes.Zonal (meridional) MF values over the Andes peak at ∼ 30 mPa (10 mPa), over the peninsula at ∼ 8 mPa (5 mPa), and over the Southern Ocean at ∼ 2 mPa (2 mPa).There are also significant differences in the relative directionality of the three regions; when considered at the seasonal scale (i.e.Fig. 11f, l, r) both the Andes and Peninsula regions at all heights exhibit stronger directional preferences, with ∼ 30 % of MF in this regions directed into a specific 30 • arc compared to a peak of ∼ 15 % over the Southern Ocean Over the Andes, MF peaks from late June to mid-August (between days 180 and 230); this continues into September for 2015.Values are consistently higher in 2015 than 2014, demonstrating significant interannual variability.Andean MFs are primarily westward at the 30 km level, shifting towards the southwest at the 40 km level before becoming more generally westward again at the 50 km level.The degree of directionality here is strongest in July, when 45 % of MF at 30 km is directed southwestwards, and weakest in May when no direction exceeds 30 % of observed MF.
Peninsula MFs peak over a much narrower window.In 2014, this is from day 180-210, and in 2015 from day 200 to 250.MFs are strongly northwestward at low altitude, becoming more variable in direction at higher altitudes.
Finally, over the Southern Ocean, MFs are much less tightly directional.With the exception of the very lowest height levels, which as described above are less reliable than others, a strong northwestwards spike is seen in the directionality distribution for July (Fig. 11o), but investigation shows this peak to be primarily due to a single large-amplitude wave in July 2014, visible in the time-smoothed Fig. 10i as a small region of enhanced MF around day 200.This single observation dominates the directional distribution for both Julys due to the very low levels seen at all other times.
In no region is a significant quantity of eastward MF seen.

Intermittency
A fundamental property of the GW field is its intermittency, or short-timescale variability.This variability is significant, with diurnal variations in MF of hundreds of percent not atypical (Hertzog et al., 2012;Plougonven et al., 2013;Wright et al., 2013).This has major implications on the background atmosphere.An MF of > 100 mPa at 35 km altitude, as seen in Fig. 5, would produce a net short-term acceleration of order hundreds of m s −1 day −1 if dissipated by a critical level near this height level, dramatically altering local atmospheric dynamics.On many other days, however, the same region of atmosphere is devoid of waves.Further, since the height at which GWs break is determined by their structure and amplitude, the driving induced by an intermittent GW field can be dramatically different to that induced by a uniformly distributed wave field.
Previous studies (Plougonven et al., 2013;Wright et al., 2013) have proposed the use of Gini's coefficient of concentration G (Gini, 1912) as a metric for this intermittency.G provides a single-number estimate of the evenness of a distribution, with 0 representing a perfectly even distribution (i.e.all wavepackets carry equal MF) and 1 representing a perfectly uneven distribution (i.e. a single wavepacket carries all MF observed).Values of G nearer to 1 thus represent regions where the total measured MF is dominated by a small number of large-magnitude wavepackets, and vice versa for values nearer 0. While G can be deficient when characterising unusual distributions, the GW MF distribution of our observations is near-lognormal in form, and thus use is robust in this context.Accordingly, we characterise the intermittency of our data in this manner, with G calculated using the method of Glasser (1962).
Figure 12a-o shows the results of this analysis as maps for each calendar month at three height levels from May to September, with the Gini coefficient calculated independently for each gridbox on a 1 • × 1 • spatial grid.We again show combined data from both 2014 and 2015.
We first observe the very high degree of spatial correspondence between Fig. 12 and Figs.7 and 8.This indicates that the large time-mean values of MF we see across the region are associated not with continuously high values of MF, but instead with large-magnitude individual wavepackets overlaid on a relatively flat background.This has previously been observed, but the very fine spatial resolution of our results at the upper height levels in our analysis reinforces this finding strongly.
Our tight spatial localisation also allows us to observe clearly that high-intermittency regions spread with height.We note, however, that the variable height discrimination of our results may introduce small biases when different height levels are compared, and such comparisons should be treated with caution at the very top and bottom of the analysed range.In particular, while the relative spatial distribution of G at   any individual height is internally consistent, precise numerical variation with height may be suspect at the very top and bottom.
At the 30 km level, regions where G > 0.5 are very tightly localised over the Southern Cone, particularly the Southern Andes, and (between June and September) the Antarctic Peninsula.At the 40 km level, these regions become much more diffuse, spreading out both zonally and meridionally.At the 50 km level, the geographic spread of regions with G > 0.5 is similar to the 40 km level, but with reduced G over the peak regions.
Figure 12p shows the same data as height series for each month, averaged over our three defined subregions.Two main clusters are seen, one at high G made up of the Andes region in all 5 months and the peninsula in July and August, and the other at low G representing the Southern Ocean in all months except July and the peninsula in May and June.A decline in G is seen at the top and bottom of the analysed range, but while this may be geophysical it may also be due to our varying height resolution and/or edge truncation effects, as discussed above, and thus it is difficult to draw conclusions based on this.
The series for August over the Southern Ocean and the peninsula in September lie between these clusters.The Peninsula series lies in this region due to the decline in GW activity in this region in September.The distribution for the Southern Ocean July is more surprising, but on detailed investigation is again due to the single large GW packet in July 2015 mentioned above, which skews the combined distribution for both Julys.This thus illustrates again the extreme intermittency of GW activity, and reinforces that while G is useful in the general case, no composite metric can fully describe the unevenness of the distribution in all cases.

Wave focusing into the stratospheric jet
As GWs propagate, they are refracted by wind shear.This refraction can very significantly change where the waves propagate and eventually break.
Previous studies using both ray-tracing (e.g Dunkerton, 1984;Eckermann, 1992) and fully GW-resolving model studies (Sato et al., 2009(Sato et al., , 2012) ) suggest that the atmospheric wind jets are a major cause of such refraction.An especially strong example is believed to be the polar night jet in our region of study, which interacts with the large orographic sources below to redirect very significant wave fluxes.In particular, Sato et al. (2012) showed that, in the GW-resolving Kanto climate model, northward GW fluxes over the Antarctic Peninsula and southward GW fluxes over the Southern Andes were consistent with latitudinal wave focusing into the polar night jet, using both ray tracing and calculated energy flux vectors.Hindley et al. (2015), meanwhile, measured GW potential energy using the COSMIC GPS-RO constellation.By careful analysis, they were able to identify a 3-D column of potential energy corresponding to the path such propagation would take.However, due to the lack of GW directional information they were unable to confirm directly that this column corresponded to the refraction seen in models.
Our measurements of full 3-D GW group velocity allow us to directly observe this refraction in observational data.This is illustrated by Fig. 13, which shows zonal wind speeds (shading, derived from ERA-Interim) and GW group velocity vectors (arrows, derived from AIRS data as described above) for two selected months, August 2014 and July 2015.These months are chosen as different wind jets are dominant in each case; the effect is, however, clearly visible throughout our dataset.
All data are averaged onto a 2 • × 2 • grid in latitude and longitude.Vertically, winds are shown on the standard ERA-Interim 60 Sigma-level grid and GWs on the 16 levels of our analysis.Data are shown as monthly means, but the effect is also seen at shorter timescales.
Considering first wind, we see clear zonal jets in both months.In August 2014, the main polar jet is centred at approximately 60 • S and 40 km altitude, with the higheraltitude midlatitude jet maximum visible at 37 • S and 60 km altitude.In July 2015, the midlatitude jet dominates, centred at 45 • S, 55 km altitude, with a very weak polar jet visible at 40 km altitude and 60 • S. Wind speeds are positive (i.e.eastward) throughout the region at all heights, falling away from their maximum towards the Equator and the pole.In both months, the jets are slightly faster on the eastern side of the region, increasing in speed by ∼ 10 ms −1 across the meridional width of the region at the peak latitude.
The first column, Figs.13a, d, shows sections along the 70 • W meridian.This longitude is chosen for consistency with Fig. 5, but is geographically close enough to the 55 • W meridian studied by Sato et al. (2012) and the 65 • W meridian studied by Hindley et al. (2015) in their equivalent figures to allow comparison.It is also close to the latitude of the 60 • S "gap" in GW MF noted by, for example, McLandress et al. (2012), and thus results here are of geophysical significance.The second column (Figs.13b, e) shows maps at the 40 km altitude level and the third column (Figs.13c, f) a 60 • S section along the Drake Passage.The 40 km height level is again chosen for consistency with Fig. 5, while the zonal section is selected as being close to the polar jet centre latitude.
In both months we see low vertical GW group speeds at the very lowest altitudes in both the zonal and meridional sections.We presume this to be due to the issues with resolution and edge-truncation described above.Above ∼ 30 km, vertical group velocities increase sharply in both months shown, maximising at the highest altitudes.
In August 2014, our results are strikingly similar to those of both Sato et al. (2012) and Hindley et al. (2015).We see a clear path in c g upwards, southwards and eastwards from the region over the Southern Andes directly into the polar jet core.This pathway is consistent with the waves being vertically refracted by the positive wind shear into the jet core while being advected downwind in the eastward direction (Dunkerton, 1984).In July 2015, the waves are instead primarily directed into the now-dominant midlatitude jet core, but otherwise behave in a similar fashion.In the south of the region in both months, c g vectors are directed weakly northwards at altitudes below the jet core, but are travelling almost entirely upwards and eastwards by the time they pass the 40 km level.
The maps show a similar result.This is most clear in August 2014, where we see c g vectors directed from the Andean regions directly into the 60 • core, but is also visible in the north of the region in July 2015.This is visually weaker, however, since the dominant jet core here is at a much higher altitude.
Finally, in the meridional sections, we again see refraction into regions of higher wind speed.This is again best visible in August 2015, but is a much weaker effect since the meridional variation of the jet speed is much smaller than the zonal or vertical.

Summary and conclusions
In this study, we develop and describe a new analysis method for detecting and characterising atmospheric gravity waves using a 3-D S-transform.We apply this method to 3-D atmospheric temperature measurements from the AIRS satellite instrument.We use these data to study in detail wave activity over the Andes/Drake Passage GW "hot spot" region, arguably the largest GW source in the world.
Our method exploits the available data to a much more complete extent than any previous technique.This lets us characterise GWs across the full range of length scales present in the input data.The method includes a correction significantly ameliorating the amplitude suppression inherent to ST methods, allowing us to estimate true wave amplitudes rather than lower bounds.The method is also extensible to higher dimensions, of use for example in the analysis of time-varying model or ground-based data.Additionally, while in this study we consider only the largest-amplitude wave at each point, the method can be easily extended to the study of multiple overlapping GWs.
Using our method, we are able to derive a full range of GW properties, including their wavevectors, amplitudes, phase and group velocities, temporal frequencies and momentum fluxes.We can associate these properties spatially with wave variations at the voxel level, i.e. ∼ tens of kilometres in the horizontal, at the time resolution of a single measurement swath rather than statistically over a period of time.Aside from a single tie-breaking assumption on the vertical propa-gation direction of waves, we are able to measure the waveintrinsic-frame values of these properties from AIRS data alone, i.e. without ancillary information.Some limitations remain in the edge-truncated regions at the very top and bottom of the height range with regards to determining precise amplitudes and vertical wavelengths, but these issues are general to the analysis of any height-restricted dataset.
Our observations show clearly that directional MF in this region is primarily directed directly upwind and is greatest above mountainous regions, consistent with previous studies (Plougonven et al., 2013;Wright et al., 2016a;Ern et al., 2017).The spatial correspondence with topography at low altitudes is extremely close, particularly over the Antarctic Peninsula.The directions of these fluxes becomes less focused with increasing height.The MF is carried primarily by a very small fraction of observed waves.Very little MF is observed over the Drake Passage.This may suggest either that non-orographic generation of waves is low in this season and region, that what waves are present lie outside the observational filter of our dataset, or some combination of these factors.We observe declining levels of |MF| above 50 km altitude, consistent with wave dissipation, but due to falling vertical resolution at these heights are unable to separate this effect from a reduction in retrieval quality.
Finally, wave group velocity vectors observed over the region are highly consistent with wave focusing into the atmospheric jets.This strongly reinforces the results of previous model studies which suggest that GWs are refracted by wind shear, and reinforces the suggestion that a large fraction of the 60 • S "wake" of strong GW activity observed downstream of the Andes in previous studies is orographic in origin (e.g.Hindley et al., 2015, and references therein).
Collectively, our results demonstrate the great potential offered by 3-D analyses of GWs.This technique thus provides a powerful tool for the study of GW physics at the global scale, providing key information needed to constrain the next generation of weather prediction and climate models.
of Hindley  et al. (2016)  and Fig.7ofHoffmann et al. (2014) transports significant MF, but is almost exactly split in half by a granule boundary.It is thus not properly resolved in either individual granule spanning the GW.

Figure 2 .
Figure 2. (a) Map of the region covered by AIRS granules 56-58 on 6 May 2008.Red solid lines indicate the AIRS scan track outline, and red arrows indicate the direction of satellite travel.(b-e) T measured by AIRS at four height levels above this region.The exponential increase in wave amplitude with altitude has been removed in (b-e) for consistency with Fig. 3.

Figure 3 .
Figure 3. 3-D plot of the AIRS-derived T measurements shown in Fig. 2, plotted as semi-transparent nested surfaces of constant T .Viewing perspective has been rotated relative to Fig. 2 to better highlight 3-D T phase fronts present in the signal, and the exponential increase in wave amplitude with altitude has been removed in order to highlight the vertical structure.

Figure 4 .
Figure 4.Estimated (a) temperature perturbations T (see text for details of annotations), (b) wave amplitude T , (c) horizontal wavelength λ h , (d) horizontal propagation angle θ , (e) vertical wavelength λ z , (f) intrinsic frequency ω, (g) intrinsic horizontal group speed ĉgh , (h) intrinsic horizontal phase speed ĉph , (i) ground-based frequency ω, (j) ground-based horizontal group speed c gh , (k) vertical group speed |c gz |, (l) ground-based horizontal phase speed c ph , and (m) vertical phase speed |c pz | at 30 km altitude for the AIRS data shown in Fig. 2. Regions with T < 5 K have been masked from (c-m).Black dotted lines indicate the edge of the AIRS measurement granule; some data extend slightly beyond this region due to the data being interpolated onto a regular grid.Intrinsic vertical group and phase speeds are not shown separately, since they are equal to the ground-based values (see Eqs. 5 and 6).Panels (g) and (j) share a colour bar, as do (h) and (l).The exponential increase in wave amplitude with altitude has been removed in (a) for consistency with Fig. 3.

Figure 5 .
Figure 5. (a, d, g) T , (b, e, h) zonal MF and (c, f, i) meridional MF associated with the waves shown in Fig. 2. Panels (a-c) show maps at 30 km altitude.(d-f) show data at 50 • S, and (g-i) 70 • W, both against height.Black shading in panels (d-i) shows surface topography, scaled by 5 for clarity.As in previous figures, the exponential increase in wave amplitude with altitude has been removed in (a, d, g) in order to highlight the vertical structure.

Figure 6 .
Figure 6.(a) Zonal and (b) meridional MF associated with the waves shown in Fig. 2.

Figure 7 .
Figure 7. Zonal MF monthly means at three height levels.

Figure 8 .
Figure 8. Meridional MF monthly means at three height levels.

Figure 10 .
Figure 10.Time variability and directionality of measured MF for three geographic subregions and 2 years.Data are daily means boxcarsmoothed by 7 days.Note different scales for each row.

Figure 11 .
Figure 11.MF directional distribution by height and month for three geographic subregions.

Figure 12 .
Figure 12.Gini coefficient G over the region of interest.(a-o) Maps of measured monthly G at (a-e) 50 km (f-j) 40 km (k-o) 30 km altitude.(p) Height series of G for each month over the three regions defined in Fig. 2.

Figure 13 .
Figure 13.Filled contours: monthly mean ECMWF-derived zonal wind (black arrows) observed GW ground-based group velocity for (a, b, c) August 2014; (d, e, f) July 2015 at (a, d) 70 • W against height and latitude; (b, e) 40 km altitude, mapped; (c, f) 50 • S against height and longitude.Length of arrows is proportional to group velocity magnitude in that direction.