Determining Gravity Wave Sources and Propagation in the Southern Hemisphere by Ray‐Tracing AIRS Measurements

Gravity waves (GWs) are key drivers of atmospheric circulation. Understanding their sources and propagation is essential to improving weather and climate models. We apply a 3D Stockwell Transform to 1 month of stratospheric temperature data from NASA's Atmospheric InfraRed Sounder to obtain 3D GW measurements and parameters. We use ray‐tracing methods to determine the sources and propagation characteristics of these GWs over the entire Southern Ocean. We trace 1.28 million GW measurements per day for the month of June 2010. Our analysis suggests that ground‐based sources around the Andes, Antarctic Peninsula, and Kerguelen play major roles, and that the GWs generated by these and other sources travel large zonal distances. We show evidence that GWs propagate into the 60°S belt, a possible source of “missing momentum flux” in GCMs at this latitude. These results emphasize the need for models to incorporate the possibility that GWs can exhibit large horizontal propagation.

. In general circulation models, these waves are assumed propagate near-vertically upwards in a model column; the GW activity routinely seen above these features is directly consistent with these assumptions. However, a significant "tail" of GWs, starting at the Andes and extending in an eastward direction almost the entire distance around the Earth (e.g., Hindley et al., 2015), has proven much more challenging to explain.
One approach to investigate the origin of this tail is to use ray-tracing. In ray-tracing, a GW with an initial set of properties is propagated either forwards or backwards in time through background wind and stability fields. Assuming that the initial parameters are chosen well and that the wind fields are accurate, backwards-traced rays should pass through the source of the wave, while forwards-traced waves show the path they take from the initial time point onwards. For early studies (e.g., Eckermann & Marks, 1996;Zhong et al., 1995), wind fields were usually coarse in resolution and the GW properties estimated from theory, but with modern reanalyses more accurate wind fields can be used, and more accurate initial parameters can now be estimated from observations. Consequently, several previous studies have used a ray-tracing approach to understand the sources and propagation characteristics of GWs over all or part of the Southern Ocean. These studies have used GW parameters and background fields sourced in various ways. Common approaches have included idealized and/or isotropic GW source spectra (Preusse et al., 2006(Preusse et al., , 2009Wells et al., 2011) or the estimation of the necessary GW parameters from model or reanalysis output (Jiang et al., 2013;Preusse et al., 2014;Sato et al., 2012;Vosper, 2015), and propagated those GWs through model or reanalysis wind fields. A less common approach has been the use of observationally derived parameters; (e.g., Pulido et al., 2013), who backwards-traced radiosonde-derived GW parameters from an Argentina-based launch campaign through winds estimated from a high-resolution WRF model simulation over the Andes. Combined, these idealized, model-based and observationally based studies suggest that (i) the downstream propagation of GWs from the Andes is proportional to their horizontal wavelength (Wells et al., 2011), (ii) that waves in this region typically propagate southwards or southeastwards into the polar jet at around 60°S (Ehard et al., 2017;Jiang et al., 2013;Sato et al., 2012;Vosper, 2015) and (iii) that many of these waves have orographic sources (Pulido et al., 2013).
This body of previous work exhibits some limitations, which the present study addresses. One key limitation is spatial coverage of the GW observations. Due to the significant computational cost of such studies and/or limited GW parametric information, almost all of the above works focus on specific limited areas, such as around small islands (Vosper, 2015) and the Andes and Antarctic Peninsula (Pulido et al., 2013;Sato et al., 2012). Studies which do cover large areas (e.g., Preusse et al., 2014) have used model-derived or idealized GW parameters. Key questions therefore remain about whether the conclusions of this earlier work (i) reflect the observed properties of real GWs and (ii) accurately describe GW propagation across the entirety of this large and important region. Accordingly, we here for the first time ray-trace observed GWs over the entire Southern Ocean region. We measure the GWs and the parameters needed to ray-trace them using satellite data from NASA's Atmospheric InfraRed Sounder (AIRS), then apply the Gravity Wave Regional Or Global Ray Tracer (GROGRAT) ray-tracer to these parameters to trace the waves back to their origins (Eckermann & Marks, 1996;Marks & Eckermann, 1995). We derive our wind fields from the ERA5 reanalysis, which is observationally well-constrained at large scales in this region (Wright & Hindley, 2018) to provide suitable accuracy for this task. We trace over 30 million measurements taken at a height of about 40 km altitude back to their sources, which we believe represents the largest observationally based GW ray-tracing study ever performed.
Section 2 describes the datasets used, and Section 3 the GW detection methods and basic concepts underlying GROGRAT. We then describe the GW measurements made and the results of ray-tracing these GWs in Section 4, before drawing conclusions in Section 5.

AIRS
The Atmospheric InfraRed Sounder (AIRS) is a nadir-sounding hyperspectral imager on NASA's Aqua satellite, part of the A-Train sun-synchronous constellation. AIRS scans a ± 49° swath either side of the orbital track (Aumann et al., 2003). This swath is processed into 90 across-track elements, with a horizontal resolution of 13.5 km at swath center falling to ∼41 km at swath edge. The data are stored every 6 min as a "granule" of 135 consecutive along-track elements, producing 240 granules of 90 × 135 elements per day.
We use data generated using the AIRS L2 retrieval of Hoffmann and Alexander (2009); this produces temperature data at the full horizontal resolution of the instrument on a 3 km vertical grid throughout the stratosphere. Retrieval error and noise are lowest in the altitude range 20-60 km (Hindley et al., 2019;Hoffmann & Alexander, 2009), so we only consider these data. In this range, retrieval error in the temperature measurements is 1.6-3.0 K, and vertical resolution varies between 7 and 15 km (Hindley et al., 2019). AIRS data are sensitive only to a portion of the GW spectrum with horizontal wavelengths less than 1,400 km and vertical wavelengths larger than 15 km (Ern et al., 2017). GWs in this spectral regime contribute heavily to total GWMF (Ern et al., 2004;Meyer et al., 2018), and thus although we are seeing only a subset of the GW spectrum, it is a subset with a significant impact on atmospheric dynamics.

ERA-5
ERA-5 is the most recent reanalysis produced by the European Center for Medium-range Weather Forecasting (ECMWF) (Ehard et al., 2018;Hersbach et al., 2020). We use ERA-5 temperatures and zonal and meridional wind reanalysis fields, at a 1.5 × 1.5° (latitude × longitude) resolution on 137 levels from the surface to 0.01 hPa at a 3 h time resolution. We vertically interpolate these global fields to a regular 1 km spacing from the surface to 75 km altitude. We further smooth the fields using a (7 × 7 × 7)-voxel boxcar in space; this is intended to remove the impact of localized temperature gradients in the background field caused by gravity waves which are directly resolved on the reanalysis grid but which may not match those in our satellite observations.

Obtaining Wave Properties-The 3D Stockwell Transform and AIRS
To provide the information needed to ray-trace measured GWs, we must both estimate their physical parameters and spatially localize where they occur. We therefore require a spectral analysis technique that localizes the wave event in both spectral and position space. We here use the Stockwell Transform (S-Transform) (Stockwell et al., 1996), which has been widely used in previous studies (Alexander et al., 2010;Hindley et al., 2016;Wright et al., 2017;Wright & Gille, 2013).
We first interpolate T′ to a regular distance grid. Following Wright et al. (2017), we then remove the exponential increase of wave amplitude with height, based at a reference point at 39 km altitude, as exp((z − z 0 )/2H). This compensates for the dominating effect of larger amplitudes at the higher altitudes of the granule on the S-Transform, and is reversed in a later step in order to ensure no net effect on measured amplitudes. We then detrend the data with a fourth-order cross-track polynomial filter independently for each height level and along-track row, consistent with previous studies using AIRS. This produces estimates of temperatures perturbations T′, which include systematic perturbations due to GWs, and of the smoothly varying background temperature T . We also apply a horizontal 3×3 voxel boxcar smoother in the horizontal plane, to reduce the effects of small-scale noise in the AIRS data. The horizontal wavelengths measured from this method range from 40 km to around 1,000 km (Hindley et al., 2019).
We next apply the 3D S-Transform to the temperature perturbation data. This produces a 6D output S(τ, where f x,y,z are spatial frequencies (inverse of wavelength) in the x, y, z directions, respectively. Following Hindley et al. (2016), we reduce this 6D field to several 3D fields. First, for each location in x, we determine the single wave with the largest spectral amplitude, giving us a complex object ( ) x  . The locations of ( ) x  within S((X), f) correspond to the wavenumbers of the dominant wave at each x. This gives us ( ), ( ), ( ) , corresponding to the f x , f y , f z frequencies that are dominant at each location in x. These four 3D fields are the same size as the input data h(x).
Using the cross-track and along-track azimuths at each point on the granule, we then in local geophysical coordinates project the components of the spatial frequencies into the three-dimensional spatial wave vector (inverse of frequency) (k, l, m), and use the magnitude of  to obtain the amplitude, from which we remove the scaling applied above. As the AIRS granule is just a single snapshot in time, there is a 180° ambiguity about the propagation direction of measured gravity waves. To break this ambiguity, we assume that all gravity waves measured are upwardly propagating (Ern et al., 2017;Wright et al., 2016Wright et al., , 2017, by requiring m to be negative. After directly fitting these parameters from the S-Transform, we compute estimates of wave phase speed and temporal frequency from the spectral data, following the method of Wright et al. (2017).

Ray-Tracing of Gravity Waves-GROGRAT
To ray-trace, we use GROGRAT. Initial GW parameters are obtained as described above, and background meteorological fields are obtained from ERA-5. Introduced by Marks and Eckermann (1995), GROGRAT is a GW ray-tracer based on the dispersion relation where U and V are zonal and meridional wind,  is the GW intrinsic frequency (i.e., the frequency an observer moving with the background wind would observe), ω is the ground-based frequency, N is the Brunt-Väisälä (buoyancy) frequency, f is the Coriolis parameter, (k, l, m) is the GW wavenumber vector, and α is 1 4H , where H is the atmospheric scale height (assumed to be 7 km).
The ray-tracing equations used in GROGRAT are derived by Lighthill (1967). The position of the ray and refraction of the wavevector along the ray in time, respectively, are computed as: Via the dispersion relation, solving these equations is based on U, V, and N background values, which we estimate using ERA-5. The version of GROGRAT used in this study also implements the spherical geometry corrections derived by Hasha et al. (2008), allowing rays to travel along great circles. This is especially important for rays traveling large horizontal distances or close to the poles. The time step dt is modulated throughout the trace to select the longest time step which produces an error below a threshold of 1%, striking a balance between computational speed and accuracy.
We backwards-trace our rays from observation to origin, but it is important to note that the source of the GW could be anywhere along the path rather than where the ray is calculated to terminate by GROGRAT (Preusse et al., 2014). To increase our certainty about the origins of our rays, we examine derived GW phase speeds and frequencies alongside measured wavelengths in combination with the inferred raypaths to better contextualize possible sources.
For the month of June 2010, we backward-trace all AIRS data in the Southern Hemisphere, from the layer of voxels at 39 km downwards toward their origins. This is at the center of our vertical measurement range, allowing us to measure long vertical wavelengths that would not be accessible to the S-Transform technique at lower or higher altitudes due to edge-truncation effects. For each voxel with a measured temperature perturbation  1.5K T , we use the 3DST to measure the wave's initial position, 2D horizontal wavenumber (k, l), ground-based frequency ω via (1), and horizontal wind amplitude  u , which we compute (Hu et al., 2002) as For each granule, we can trace a maximum of 90 × 135 = 12,150 individual measurements. On average after allowing for the above amplitude filtering, we actually trace paths from 1.28 million voxels per day.
PERRETT ET AL.
10.1029/2020GL088621 Figure 1 demonstrates the raypaths generated from tracing a single AIRS granule observed near the Drake Passage, measured on June 4, 2010 at around 04:10 UTC. GW T′ at 39 km is shown in Figure 1a. We see a chevron-like pattern, indicative of potential GW sources in both the Antarctic Peninsula and Southern Andes. Figure 1b shows rays produced by ray-tracing each voxel in this 39 km layer backwards in time. We see rays from the northern half of the granule terminating near the surface around the Andes, and rays from the southern half of the granule terminating around the Antarctic Peninsula. This hence demonstrates simultaneous northward propagation of GWs from the Antarctic Peninsula and southward propagation of GWs from the Andes into the Drake Passage region. This has been previously been theoretically demonstrated by Sato et al. (2012), using idealized background wind and initial GW parameters. We obtain a similar height of convergence as Sato et al. (2012), who saw this occurring at ∼40 km altitude. The raypaths become increasingly vertical with height, as has previously been seen using 3D group velocity measurements over the Andes by Wright et al. (2017).

Case Study over Southern Andes
Using simulated GWs in a high-resolution model, Sato et al. (2012) also suggested that GWs over the Andes with horizontal wavenumber angles around 220° will propagate roughly 10° in longitude from the Andes by the time they reach 40 km. For the northern portion of this granule over the Andes, we measure a horizontal wavenumber angle of 215°, anticlockwise from east. We see similar propagation characteristics as Sato et al. (2012) in this case study.

Ray-Tracing Over the Southern Ocean-June 2010
Building upon the case study of section 4.1, we proceed to ray-trace all AIRS granules over the Southern Ocean in June 2010. This month was chosen as the wave activity was strongest for 2010 austral winter (Hindley et al., 2019). This allows us to examine GW sources and propagation characteristics over the entire region. We again note that we are inferring GW sources from a combination of AIRS derived GW phase speeds and frequencies, AIRS measured wavelengths and the characteristics of the raypaths. Figure S2 shows the distribution of mean ray launch and termination points. The much denser regions over the Antarctic continent are due to a much greater number of overpasses being averaged, due to the polar orbit of AIRS. As this region is in near constant darkness, there are also many more observations included in our average per satellite overpass.
We first present the 3DST-derived GW parameters used as input for the ray-tracing analysis. These are shown in Figure 2. Data are shown as the monthly mean of all voxels, on a 1.5 × 1.5° grid. Voxels with an PERRETT ET AL.
10.1029/2020GL088621 5 of 10 amplitude less than 1.5 K have been removed from all variables before both plotting and ray-tracing, consistent with the instrument noise analysis presented in Figure 2 of Hindley et al. (2019). We also remove all data poleward of 85°S, as Wright and Hindley (2018) showed very low correlation (r < 0.4) between ERA5 and satellite temperature observations at the initial 39 km height this close to the pole. Since ray-tracing relies on the accuracy of the background fields, using data from this region may provide results much less reliable than for other regions. Figure 2a shows average temperature amplitudes. Note that these temperature amplitudes are converted to wind amplitudes before the ray-tracing analysis; however, we show the temperature amplitudes here since they are the measured quantity. We see a primary peak over the southern Andes and Antarctic Peninsula, extending ∼80° downstream to the east. We also see a secondary maximum over Kerguelen (49°S, 69°E), extending ∼60° eastwards. There is also some evidence of increased activity relative to the surrounding region over South Georgia (54°S, 36°W) and New Zealand; however, both maxima represent only a small increase above the background, and the maximum from South Georgia in particular is hard to distinguish due to the long tail of Andean and Antarctic Peninsula activity. We also see longer wavelengths over the open ocean, with smaller wavelengths occurring over topography. We observe some indication of an inward spiral of shorter wavelengths over the Southern Ocean starting at the Andes and passing over Kerguelen into the Admiralty Mountains, as seen by previous studies (e.g., Hendricks et al., 2014;Hindley et al., 2019Hindley et al., , 2015McLandress et al., 2012;Sato et al., 2012).  from measured properties using linear GW theory (Fritts & Alexander, 2003;Wright et al., 2017). The values here are consistent with Figures 2a and 2b, with a region of particularly low frequency consistent with orographic generation over the southern Andes. Although for a different year, these values agree quantitatively with Wright et al. (2017). We more clearly see the spiral feature described by previous work, particularly in Figure 2d and, with a secondary minimum over Kerguelen. In Figure 2d, we see regions of lower (<30 m s −1 ) phase speed around the southern Andes, Antarctic Peninsula, South Georgia and the Kerguelen Islands, indicative of orographic GW activity. For true orographic generation, measured phase speeds over these regions would be zero; this difference is likely due to the combined effects of nonorographic GWs from storms and jets being included in our averages (Hoffmann et al., 2016;Preusse et al., 2014), and measurement uncertainties propagating through the calculation (Wright et al., 2017). Figure 3a shows the monthly mean Lowest Traceable Altitude (LTA) reached by each ray, averaged on a 1.5 × 1.5° grid. This is the lowest altitude the raypath reaches when back-tracing, before being terminated by GROGRAT due to either reaching the ground or entering an atmospheric regime where the wave cannot be physically supported. Data shown here have been geolocated at the termination point of the ray. We see distinct regions of low LTA around the southern Andes, the Antarctic Peninsula, and Kerguelen, where the mean LTA is <4 km, and less than <2 km for the Andes. There is also a small region of sub-4 km altitude LTA near South Georgia; this is the only other location to have an LTA this low. These LTA-minima are indicative of GW sources in these regions originating near to or at the surface, averaged together with other nonsurface sources. A similar analysis by Preusse et al. (2014) demonstrated regions of low LTA over the Andes and Antarctic Peninsula, however these were much less localized, and did not display the clear sources over small mountainous islands we see in Figure 3a.

GW Sources and Propagation Characteristics through Ray-Tracing
There are also regions of low LTA over the Antarctic continent. However, when considered in the context of the high phase speeds seen over these regions (Figure 2d, it is less probable that these waves are orographic in origin. When we consider the higher phase speeds measured in this region and the small distances traveled (Figures 3b and 3c), this suggests that this could be due to GWs originating from near-surface nonorographic sources. However, low LTAs over Antarctica have also been seen by Mehta et al. (2017), with higher phase speed measurements originating near the surface around the south pole. Preusse et al. (2014) also observed regions of low LTA over the south pole, possibly due to Katabatic winds.  geographically local GW sources are to where the waves are measured in the midstratosphere. Comparing to Figure 2b, we see a positive correlation between horizontal wavelength and zonal propagation, agreeing with previous studies (Wells et al., 2011). We see two patches where GWs have on average traveled large zonal distances.
Over Kerguelen, the zonal distance traveled is close to zero, indicative of orographic waves generated in this region that are ascending vertically upwards through the atmosphere.
In contrast, higher zonal distances are seen for GWs originating over the Andes. The underlying cause of this effect is demonstrated in Figure 1: examination of the direct ray-traces (omitted for brevity) shows that many rays in this region travel long horizontal distances in the troposphere and lower stratosphere before traveling vertically upwards toward their measurement altitude, rather than directly upwards over the Andes.
All GW traces move eastward in the Southern Hemisphere ( Figure S3). In austral winter, the Southern Ocean is known for strong, westerly winds, regularly exceeding 60-80 m s −1 in the stratosphere ( Figure S1); therefore, under our initial assumption that measured GWs are upwardly propagating, this is consistent with the waves being refracted by these strong, westerly winds as they travel from their lowest point to the measurement altitude. Figure 3c shows that GWs measured in the patches of larger GW activity propagate southwards, converging around 60°S, south of which they are replaced by waves which have on average traveled only very short distances both zonally and meridionally with high LTAs. Previous studies have implied this result through measurements of converging wave direction at this latitude, both in observations and models (Hindley et al., 2019;Holt et al., 2017;Sato et al., 2012). Propagation of GWs into this region is hypothesized to be one of the key reasons global climate models underestimate momentum flux around the Southern Ocean (Hindley et al., 2015(Hindley et al., , 2019McLandress et al., 2012), and thus this further demonstrates the need for climate models to model meridional propagation of gravity waves.

Conclusion
We have performed what we believe to be the largest-ever backwards ray-tracing study of measured atmospheric gravity waves, focusing on the region above the Southern Ocean and Antarctica. Building upon previous work, we backwards-trace GW rays from mid-stratospheric measurements to the lower atmosphere, allowing us to determine their likely sources. We show clear evidence for significant surface GW sources around the Andes, Antarctic Peninsula and Kerguelen. In comparison with previous ray-tracing studies, we see much finer localization of sources, especially over orography.
In ray-tracing the waves backwards, we find that the GWs observed by AIRS often travel many hundreds of kilometers zonally from their source location. Since many of these waves will be unresolved in climate models, this suggests that the standard assumption in GW parameterization that subgrid scale GWs travel with no significant horizontal propagation is a poor approximation in this region. Our results clearly demonstrate the need for improvements in these schemes which will allow subgridscale GWs to propagate in a more realistic manner.
We also show that GWs around the Southern Ocean region propagate into and converge around a belt at 60°S. This builds on previous studies, providing further clear evidence that this convergence happens all around the region. Missing propagation of GWs into this belt is hypothesized to be one of the key reasons global climate models underestimate momentum flux around the Southern Ocean.

Data Availability Statement
AIRS data used in this study were originally obtained from NASA's DISC service, and have been preprocessed as described in Hoffmann and Alexander (2009), and are available from https://datapub.fz-juelich. de/slcs/airs/gravity_waves/data/jon/. ERA5 data were downloaded from the Copernicus Climate Change Service Climate Data Store: https://cds.climate.copernicus.eu/. This research made use of the Balena High Performance Computing (HPC) Service at the University of Bath.