A two-dimensional Stockwell transform for gravity wave analysis of AIRS measurements

Abstract

Gravity waves (GWs) play a crucial role in the dynamics of the earth's atmosphere. These waves couple lower, middle and upper atmospheric layers by transporting and depositing energy and momentum from their sources to great heights. The accurate parameterisation of GW momentum flux is of key importance to general circulation models but requires accurate measurement of GW properties, which has proved challenging. For more than a decade, the nadir-viewing Atmospheric Infrared Sounder (AIRS) aboard NASA's Aqua satellite has made global, two-dimensional (2-D) measurements of stratospheric radiances in which GWs can be detected. However, one problem with current one-dimensional methods for GW analysis of these data is that they can introduce significant unwanted biases. Here, we present a new analysis method that resolves this problem. Our method uses a 2-D Stockwell transform (2DST) to measure GW amplitudes, horizontal wavelengths and directions of propagation using both the along-track and cross-track dimensions simultaneously. We first test our new method and demonstrate that it can accurately measure GW properties in a specified wave field. We then show that by using a new elliptical spectral window in the 2DST, in place of the traditional Gaussian, we can dramatically improve the recovery of wave amplitude over the standard approach. We then use our improved method to measure GW properties and momentum fluxes in AIRS measurements over two regions known to be intense hotspots of GW activity: (i) the Drake Passage/Antarctic Peninsula and (ii) the isolated mountainous island of South Georgia. The significance of our new 2DST method is that it provides more accurate, unbiased and better localised measurements of key GW properties compared to most current methods. The added flexibility offered by the scaling parameter and our new spectral window presented here extend the usefulness of our 2DST method to other areas of geophysical data analysis and beyond.

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/ doi:10.5194/amt-9-2545-2016 Author(s) 2016. CC Attribution 3.0 License.

Neil P. Hindley1, Nathan D. Smith1, Corwin J. Wright1, D. Andrew S. Rees2, and Nicholas J. Mitchell1

1Centre for Space, Atmosphere and Ocean Science, University of Bath, Bath, UK

2Department of Mechanical Engineering, University of Bath, Bath, UK

Correspondence to: Neil P. Hindley ([email protected])

Received: 7 December 2015 Published in Atmos. Meas. Tech. Discuss.: 18 January 2016 Revised: 5 May 2016 Accepted: 6 May 2016 Published: 9 June 2016

Abstract. Gravity waves (GWs) play a crucial role in the dynamics of the earths atmosphere. These waves couple lower, middle and upper atmospheric layers by transporting and depositing energy and momentum from their sources to great heights. The accurate parameterisation of GW momentum ux is of key importance to general circulation models but requires accurate measurement of GW properties, which has proved challenging. For more than a decade, the nadir-viewing Atmospheric Infrared Sounder (AIRS) aboard NASAs Aqua satellite has made global, two-dimensional (2-D) measurements of stratospheric radiances in which GWs can be detected. However, one problem with current one-dimensional methods for GW analysis of these data is that they can introduce signicant unwanted biases. Here, we present a new analysis method that resolves this problem. Our method uses a 2-D Stockwell transform (2DST) to measure GW amplitudes, horizontal wavelengths and directions of propagation using both the along-track and cross-track dimensions simultaneously. We rst test our new method and demonstrate that it can accurately measure GW properties in a specied wave eld. We then show that by using a new elliptical spectral window in the 2DST, in place of the traditional Gaussian, we can dramatically improve the recovery of wave amplitude over the standard approach. We then use our improved method to measure GW properties and momentum uxes in AIRS measurements over two regions known to be intense hotspots of GW activity: (i) the Drake Passage/Antarctic Peninsula and (ii) the isolated mountainous island of South Georgia. The signicance of our new 2DST method is that it provides more accurate, unbiased and better localised measurements of key GW properties compared to most current methods. The added exibility offered by the

A two-dimensional Stockwell transform for gravity wave analysis of AIRS measurements

scaling parameter and our new spectral window presented here extend the usefulness of our 2DST method to other areas of geophysical data analysis and beyond.

1 Introduction

Gravity waves are a vital component of the atmospheric system. These propagating mesoscale disturbances can transport energy and momentum from their source regions to great heights. They thus are a key driving mechanism in the dynamics of the middle atmosphere through drag and diffusion processes (e.g. Fritts and Alexander, 2003, and references therein).

The accurate parameterisation of unresolved gravity waves in global climate models (GCMs) has proven to be a longstanding problem in the modelling community. One example of this is the cold pole bias (Butchart et al., 2011) suffered by nearly all GCMs, which has been suggested to be due to a deciency of resolved and parameterised gravity wave drag near 60 S. One reason for this is that these gravity wave parameterisations continue to be poorly constrained by observations (Alexander et al., 2010). The accurate measurement of gravity wave properties is thus critical for the development of the next and current generation of climate models. In the last decade, satellite-based remote-sensing has greatly increased our capability to make gravity wave observations on a global scale, but large discrepancies between observed and modelled uxes still remain (Geller et al., 2013). This highlights the need for more accurate global and regional gravity wave measurements.

Published by Copernicus Publications on behalf of the European Geosciences Union.

2546 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

The Atmospheric Infrared Sounder (AIRS) (Aumann et al., 2003) is a nadir-sounding spectral imager on board the Aqua satellite, launched in 2002. Part of the A-Train satellite constellation, AIRS scans the atmosphere over the range [notdef]49 from the nadir of the satellite in a 90 pixel

( 1800 km) wide swath, using 2378 infrared channels along

a sun-synchronous polar orbit. This continuous swath is archived in granules, usually 135 pixels ( 2400 km) along-

track. Stratospheric gravity waves can be detected in these granules as radiance perturbations in the 15 and 4.3 m CO2 emission bands (e.g. Alexander and Barnet, 2007; Hoffmann and Alexander, 2009; Hoffmann et al., 2013). AIRS measurements enable the study of stratospheric gravity waves at unprecedented horizontal resolution. In order to fully exploit these observations, accurate and easily reproducible analysis methods for the measurement of gravity wave properties must be developed.

The Stockwell transform (S-transform) (Stockwell et al., 1996; Stockwell, 1999) is a widely used spectral analysis technique for providing time-frequency (or distance-wavenumber) localisation of a time series (or spatial prole).This capability makes the S-transform well suited to gravity wave analysis of a variety of geophysical data (e.g. Fritts et al., 1998; Alexander et al., 2008; McDonald, 2012; Wright and Gille, 2013). The S-transform has also been used in a variety of other elds, such as the planetary (Wright, 2012), engineering (Kuyuk, 2015) and medical sciences (Goodyear et al., 2004; Yan et al., 2015).

Alexander and Barnet (2007) developed a method for measuring gravity wave amplitudes, horizontal wavelengths and directions of propagation from AIRS granules using the one-dimensional (1-D) S-transform. In their method, the S-transform is computed for each cross-track row, and co-spectra between adjacent cross-track rows are used to obtain spectral information in the along-track dimension. To nd the dominant waves in each granule, these cross-track co-spectra are averaged together and up to ve peaks are located in each averaged spectrum. The method of Alexander and Barnet (2007) provides good rst-order measurement of the properties of the (up to ve) dominant wave features in a granule, but it can introduce unwanted biases as discussed further in Sect. 5.

Here, we present a new analysis method. AIRS radiance measurements are two-dimensional (2-D) images; thus a gravity wave analysis method using a two-dimensional Stockwell transform (2DST) is a more logical approach. In this study we present a 2DST-based method for the measurement of gravity wave amplitudes, horizontal wavelengths and directions of propagation from AIRS measurements. Our method takes advantage of the spatialspectral localisation capabilities of the S-transform in both dimensions simultaneously, equally and without bias.

South Georgia and the Antarctic Peninsula, together with the southern tip of South America, lie in a well-known hotspot of stratospheric gravity wave activity during austral

winter, which has been extensively studied both observationally (Eckermann and Preusse, 1999; Jiang et al., 2002; Alexander and Teitelbaum, 2007; Baumgaertner and Mc-Donald, 2007; Hertzog et al., 2008; Alexander et al., 2009; Alexander and Teitelbaum, 2011; Alexander and Grimsdell, 2013; Hindley et al., 2015) and with numerical modelling techniques (Hertzog et al., 2008; Plougonven et al., 2010; Shutts and Vosper, 2011; Hertzog et al., 2012; Sato et al., 2012; Plougonven et al., 2013) in the last decade. These mountainous regions are subjected to a strong wintertime circumpolar ow in the troposphere and stratosphere and, as a result, are major orographic gravity wave sources (e.g. Hoff-mann et al., 2013). Despite this, discrepancies between ob-served and modelled gravity wave uxes in this region are the largest anywhere on the planet (Geller et al., 2013). This unique geography of this important region provides a natural laboratory in which to make clear gravity wave measurements from space and is thus an ideal region in which to test our 2DST methodology on AIRS measurements.

In Sect. 2 we introduce AIRS data and the 1-D and 2-D S-transforms. In Sect. 3 we apply the 2DST to a specied wave eld, describing our methodology for spatial localisation of the dominant spectral components. In Sect. 4 we present a new alternative spectral window for use in the 2DST. In Sect. 5 we apply the 2DST to two selected AIRS granules measured over the Antarctic Peninsula/Drake Passage and South Georgia and discuss our results. Finally, in Sect. 6 we summarise the key results of this study and discuss the advantages of our 2DST method in the context of previous work.

2 Data and spectral methods

2.1 AIRS data

Gravity waves can be detected in AIRS radiance measurements as perturbations from a background state. Here, we use AIRS Level 1B radiance measurements from the 667.77 cm1 channel. These Level 1 radiances have considerably higher horizontal resolution than operational Level 2 temperature retrievals due to retrieval choices imposed on the latter (Hoffmann and Alexander, 2009). We compute brightness temperature T directly from radiance R as

T =

2hc2 3R + 1[parenrightbigg][parenrightbigg]1, (1)

where h is Plancks constant, c is the speed of light, kB is Boltzmanns constant and = 667.77 cm1 is the wavenum

ber of the specied channel. Temperature perturbations T [prime] from the local background state T are then extracted via a fourth-order polynomial t (Wu, 2004; Alexander and Barnet, 2007). This t removes limb-brightening and other large-scale uctuations. These brightness temperature perturbations are a more useful physical quantity with which to de-

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

hc kB

ln

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2547

ne gravity wave amplitudes, since gravity wave energies and momentum uxes are easier dened in terms of temperature perturbations (e.g. Ern et al., 2004). Common limitations in temperature perturbation-based momentum ux calculations, such as the assumption of a monochromatic wave and observational ltering effects, do remain, however, and are discussed in more detail in Sect. 5.3.

The weighting function of the 667.77 cm1 channel peaks near 3 hPa ( 40 km), with a full width at half maximum

of 12 km (Alexander and Barnet, 2007, also illustrated in

Fig. 1 of Wright et al., 2016a). Gravity waves with vertical wavelengths shorter than 12 km are thus unlikely to be resolved and vertical wavelengths close to this limit will be strongly attenuated.

If the vertical wavelength is known, it is possible to correct for this attenuation by dividing the amplitude by an appropriate rescaling factor (Alexander and Barnet, 2007, their Fig. 4). Although methods for measuring long vertical wavelengths using multiple AIRS channels have been developed (e.g. Hoffmann and Alexander, 2009), we do not have direct measurements of vertical wavelength from our single AIRS channel, and so we do not apply such a correction to brightness temperature perturbations at this stage. The true amplitude of some waves in our initial analysis may therefore be between 2 and 5 times greater than the values shown. Later, however, for the estimation of momentum ux in Sect. 5.3, we do apply the attenuation correction described in Alexander and Barnet (2007).

2.2 The Stockwell transform

In its analytical form, the 1-D Stockwell transform (Stock-well et al., 1996) closely resembles a continuous wavelet transform (CWT) with a complex sinusoidal mother wavelet windowed with a scalable Gaussian window (Gibson et al., 2006). For time series data, this scalable Gaussian localises wave perturbations in the time domain through spectral localisation in the frequency domain.

For a smoothly varying, continuous and one-dimensional function of time h(t), the generalised analytical form of the S-transform S(,f ) (e.g. Pinnegar and Mansinha, 2003) is given as

S(,f ) =

where is the standard deviation. A key aspect of the Gaussian window in Eq. (3) is that the standard deviation is scaled for each frequency as =

c

|f [notdef] , where c is a scaling parame

ter usually set to 1 (Mansinha et al., 1997a). This window is often referred to as the voice Gaussian, which provides localisation of a specic frequency voice (Stockwell, 1999). Another key aspect of the Gaussian in Eq. (3) is the normalisation factor 1/p2, which ensures that the integral of the window over all t is equal to unity, a requirement for any windowing function used in the S-transform. Substituting Eq. (3) into Eq. (2) allows us to write the S-transform more explicitly as

S(,f ) = [notdef]

f [notdef] cp2

1 h(t)e

(t)

2f 22c2 ei2ftdt. (4)

Typically, the scaling parameter c is set to 1 (e.g. Stock-well et al., 1996; Alexander et al., 2008; Wright and Gille, 2013), but it may also be set to other values to achieve more specic time-frequency localisation requirements (e.g. Mansinha et al., 1997b; Fritts et al., 1998; Pinnegar and Mansinha, 2003). Setting c > 1 provides enhanced frequency localisation at the expense of time localisation, and contrarily setting c < 1 achieves enhanced time localisation at the expense of frequency localisation. This effect is discussed in more detail in Sect. 4.

To compute the S-transform using the form in Eq. (4), it seems we must compute a convolution involving the voice Gaussian and the time series for each frequency voice f , which can become quite computationally intensive. Stock-well (1999, their Sect. 3.5.1) showed that under the convolution theorem (Brigham, 1974), the time-domain convolution in Eq. (4) could be written as a frequency-domain multiplication as

S(,f ) =

1 H( + f )e

22c2 2f 2 ei2 d , (5)

where H( + f ) is a shifted version of H( ), which is in

turn the frequency analogue of H(t). The frequency-domain form of the voice Gaussian, denoted by !g( ,f ), is given as

!g( ,f ) = e

22c2 2f 2 . (6)

The standard deviation of this frequency-domain Gaussian window in Eq. (6) scales with frequency as = [notdef]f [notdef]/c.

Note that this voice Gaussian is unnormalised; its peak value is equal to 1 in the frequency domain.

In this frequency-domain form, the S-transform is computed for each frequency voice f as the inverse Fourier transform of the product of H( + f ) and the correspond

ing frequency-domain voice Gaussian !g( ,f ) in Eq. (6).

Crucially, writing the S-transform as the frequency-domain multiplication in Eq. (5) enables computationally efcient

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

1 h(t)!g(t ,f )ei2ftdt, (2)

where is translation in the time domain, f is frequency and !g(t ,f ) is a windowing function, scaled with frequency,

that provides spatial and spectral localisation. Traditionally, !g(t ,f ) takes the form of the normalised Gaussian win

dow

!g(t ,f ) =

1 p2

e

(t)

222 , (3)

2548 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

(fast) discrete Fourier transform (DFT) algorithms and simple multiplication operations to be used. The S-transform is most commonly implemented in this manner within the atmospheric sciences.

The S-transform has a number of desirable characteristics for geophysical data analysis. Unlike a CWT, the absolute magnitudes of the complex-valued S-transform coefcients in S(,f ) are directly related to the true underlying amplitude of the corresponding frequency voice f at each location . Information regarding wave amplitude is not strictly recoverable from a CWT, since the corresponding CWT coefcients are pseudo-correlation coefcients between the signal and the analysing wavelet.

One disadvantage to using fast DFT algorithms in an S-transform implementation is the familiar coarse wavelength resolution at low frequencies, a limitation not encountered by the CWT. Since both the S-transform and DFT algorithms are easily extended to higher dimensions, however, the reduced computational expense of a DFT-based S-transform makes this a practical tool for large 2-D data sets. Retention of the wave amplitude information in the S-transform is another key advantage.

2.3 The two-dimensional Stockwell transform

The S-transform is easily extended to higher dimensions.

For a two-dimensional image h(x,y), the 2DST is given by (Mansinha et al., 1997a; Stockwell, 1999)

S(x,y,fx,fy) =

1

[integraldisplay]

1

1

[integraldisplay]

1

This is the 2-D form of the 1-D Gaussian window in Eq. (3). Here, the standard deviations of the 2-D Gaussian window in Eq. (9) are scaled with wavenumber in the x and y directions as 2c/[notdef]kx[notdef] and 2c/[notdef]ky[notdef], where c is a scaling

parameter.

As discussed in Sect. 2.2, greater computational efciency is achieved by computing the 2DST as an operation in the wavenumber domain as

S(x,y,kx,ky) =

1 1 H( x + kx, y + ky)

1 1 [parenrightbigg] ei( xx+ yy)d xd y, (10)

where H( x +kx, y +ky) is a shifted version of H( x, y),

which is in turn the wavenumber analogue of the input image h(x,y). Here wavenumbers kx and ky are used to scale, in x and y directions respectively, the standard deviations of the wavenumber-domain form of the 2-D voice Gaussian Wg( x, y,kx,ky), which is given as

Wg( x, y,kx,ky) = e

22c2 2x

k2x

+ 22c2 2 y k2y

[notdef]

1 42 e

22c2 2x

k2x

+ 22c2 2 y k2y

. (11)

The 2DST is introduced and well described by Mansinha et al. (1997a) and Mansinha et al. (1997b), who demonstrated its promise for pattern analysis. It has since been discussed and applied in a variety of elds (Liu and Wong, 2007; Kocahan et al., 2008; Liu, 2009; Barry et al., 2012) but to our knowledge it has yet to be used for geophysical data analysis in the atmospheric sciences, despite the wide use of the 1-D form. In the following section, we describe our 2DST implementation methodology for the purpose of gravity wave analysis from 2-D data.

3 2DST analysis of a specied wave eld

To assess the capabilities of the 2DST, it is logical to rst apply it to a two-dimensional specied wave eld containing synthetic waves with known characteristics.

We create a specied wave eld h(x,y) with dimensions 100 [notdef] 100 km containing synthetic waves with unit ampli

tudes and known wavelengths. Wave amplitudes are dened as temperature perturbations T [prime] in units of Kelvin. The synthetic waves are localised around their central locations with Gaussian functions (although note that they do overlap). We also add random (salt and pepper) noise up to 10 % of the wave amplitude.

We rst compute the 2-D DFT H( x, y) of our specied

wave eld h(x,y). To recover an estimate of the underlying wave amplitude, we use the familiar symmetry around the zeroth frequency in the Fourier domain to recover a 2-D analogy of the analytic signal, following the approach of Stockwell (1999). A 2-D DFT contains four quadrants that

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

h(x,y)

[notdef] [notdef]

fx[notdef][notdef]fy[notdef]

2c2 e

(xx )2f 2x +(yy )2f 2 y 2c2

ei2(fxx+fyy)dx dy, (7)

where x, y are translation in the x and y directions respectively. Here, fx and fy are simple spatial frequencies (inverse of wavelength) in the x and y directions respectively, following the notation of Stockwell (1999). For the remainder of the present paper, however, we switch to using angular wavenumbers kx = 2fx and ky = 2fy, since this notation

is more commonly used in the atmospheric sciences. Rewriting Eq. (7) in terms of angular wavenumbers kx and ky gives

S(x,y,kx,ky) =

1

[integraldisplay]

1

1

[integraldisplay]

1

h(x,y)

ei(kxx+kyy)dx dy. (8)

The Gaussian windowing term in Eq. 8 describes the 2-D voice Gaussian wg(x x,y y,kx,ky), where

wg(x,y,kx,ky) = [notdef]

kx[notdef][notdef]ky[notdef]

83c2 e

[notdef] [notdef]

kx[notdef][notdef]ky[notdef]

83c2 e

(xx )2k2x +(yy ) 2k2 y 82c2

y y282c2 . (9)

k2

x x2+k

2

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2549

Figure 1. The specied wave eld h(x,y) (background) for which the two-dimensional Stockwell transform (2DST) has been computed. The absolute magnitudes of the localised 2DST wavenumber spectra [notdef] (kx,ky)[notdef] (foreground) are plotted for three separate locations in (a),

(b) and (c).

contain coefcients which are in complex-conjugate pairs with the coefcients in the opposite quadrant. The sum of these pairs always yields a real signal. By setting the coefcients of two of these quadrants to 0, and doubling their opposite quadrants, we obtain a complex-valued image when we take the inverse DFT. The magnitude of this image is analogous to the underlying wave amplitude, while the complex part describes instantaneous phase. All coefcients not in a complex conjugate pair are unchanged. The full 2DST spectrum S(x,y,kx,ky) can then be computed by taking the inverse 2-D DFT of the product of the shifted spectrum H( x kx, y ky) and the corresponding voice Gaussian

Wg( x, y,kx,ky) for each wavenumber voice kx and ky.For our purposes, we do not need to evaluate

S(x,y,kx,ky) for all positive and negative wavenumbers, since evaluating all positive and negative values of ky and only the positive values of kx gives us all the degrees of freedom. There is a residual 180 ambiguity in wave propagation direction which cannot be broken without additional information, which is supplied in Sect. 5.3.

A useful aspect of our implementation is that, like the

CWT, we can compute the 2DST for any individual or range of permitted wavenumber voices by applying the appropriate wavenumber-scaled Gaussian windows. Although the permitted wavenumber voices in the spectral domain are evenly spaced, their corresponding wavelengths are limited to integer fractions (i.e. 1/n where n = 0,1,2,...,N 1) of the

number of elements N in each dimension. This is an unavoidable consequence of using computationally efcient DFT algorithms, which results in the familiar coarse spectral resolution seen at long wavelengths.

The ability to analyse an image at specic wavenumbers is a desirable aspect in geophysical data analysis, where some a priori information regarding the spectral range of wavenumbers detectable in a given data set can be used to reduce the impact of unphysical, spurious or noisy results in 2DST analysis.

3.1 Measuring gravity wave properties

The 2DST S(x,y,kx,ky) of our specied wave eld is a four-dimensional (4-D) object. For each location in h(x,y), a two-dimensional complex-valued image of the localised spectral coefcients (kx,ky) is evaluated.

Figure 1 shows a specied wave eld h(x,y) for which the 2DST has been computed. The absolute magnitude of the localised two-dimensional wavenumber spectrum [notdef] (kx,ky)[notdef]

is plotted for three different example locations. The coefcients of [notdef] (kx,ky)[notdef] can be directly interpreted as the under

lying amplitudes of waves with wavenumbers kx and ky at a given location in the specied wave eld. As discussed in previous studies (e.g. Wright and Gille, 2013; Wright et al., 2015), there are likely to be multiple peaks in [notdef] (kx,ky)[notdef]

corresponding to overlapping waves at the same location in h(x,y). Indeed, in Fig. 1b we examine a location in the spec-ied wave eld where a small, high-wavenumber wave is located at the intersection of four lower-wavenumber waves. The localised spectrum computed by the 2DST shown in the foreground represents this feature well. The maximum spectral response is located in a peak at high kx and ky wavenumbers, with four smaller spectral peaks at lower wavenumbers with lower spectral responses.

A 4-D complex-valued function can be difcult to visualise. A more useful product might be a series of two-dimensional images, the same size as the input image, that contain the characteristics of the dominant wave at each location. In the implementation presented here, we neglect overlapping waves and identify a single dominant wave for each location in h(x,y).

For each such location, we record the complex coefcient of (kx,ky) located at the spectral peak of [notdef] (kx,ky)[notdef]. This

yields one complex-valued image (x,y), with the same dimensions as the specied wave eld h(x,y), which contains the amplitude and phase of the dominant wave at each location.

The location of the spectral peak in [notdef] (kx,ky)[notdef] also gives

us the wavenumbers kx and ky to which this peak coefcient corresponds. Hence, we can produce two further im-

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

2550 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

Figure 2. The specied wave eld (a), containing synthetic waves numbered 18, for which our two-dimensional Stockwell transform analysis has been performed. The reconstructed wave eld T [prime]2DST, underlying wave amplitudes [notdef]T [prime][notdef]2DST, horizontal wavelengths H and

directions (measured anticlockwise from the positive x direction) are shown in (b), (c), (d) and (e) respectively. Distances and wavelengths have units of kilometres and amplitudes have units of Kelvin. Panel (f) compares input and measured wavelengths (blue numbered dots) and input and measured propagation angles (red numbered dots) for the eight synthetic waves. The dashed grey line in (f) shows 1 : 1 correspondence.

ages Kx(x,y) and Ky(x,y) which contain the dominant wavenumbers at each location in the specied wave eld to which the coefcients of (x,y) correspond.

Thus, in the three images (x,y), Kx(x,y) and

Ky(x,y), we can measure the amplitudes, phases, wavelengths and propagation directions of the dominant wave features at each location in our specied wave eld.

Figure 2a shows our specied wave eld h(x,y). The central locations of the eight synthetic waves with unit amplitudes and known wavelengths are numbered 18.

By taking the real part of the complex-valued image (x,y) containing the dominant coefcients, we can recover a reconstruction of the specied wave eld, which is shown in Fig. 2b. This is made possible by the approach described above in Sect. 3.1.

The 2DST identies the different spectral regimes of the specied wave eld very well, but the reconstructed wave amplitudes are reduced by comparison to their original values.

We suspect the main reason for the reduced amplitudes relates to the spreading of spectral power in the transform. Here, as is often the case for gravity waves in the real world, our simulated waves form small wave packets, where wave amplitude decreases around a central location. Such wave packets are usually represented in the spectral domain as

some combination of wavenumber voices, in addition to the dominant wavenumber of each of the packets, in order to accurately describe their spatial properties. This means that the spectral power of a single, non-innite wave packet can be spread across multiple wavenumber voices. Spectral leakage can further contribute to this effect.

The Gaussian window in the 2DST is equal to 1 at its central location but immediately falls away with increasing radius. This means that any spectral power contained in adjacent wavenumber voices, which is required to fully reconstruct the wave, is reduced. When the inverse DFT is computed, the recovered wave amplitude at this location is thus often diminished.

A further reason for the diminished amplitude recovery in

Fig. 2b is due to wave undersampling. This undersampling effect is worse for longer wavelengths, since fewer wave cycles are present in the same-sized region of the image. The wave undersampling limitations of the S-transform are well understood in one dimension (Wright, 2010; Wright et al., 2015).

Figure 2c shows the absolute magnitude of the complex-valued image [notdef](x,y)[notdef], which corresponds to the full un

derlying amplitude of the dominant wave at each location.

This output is useful for dening regions of the specied

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2551

wave eld that do or do not contain clear and obvious wave features (McDonald, 2012).

The horizontal wavelength H (x,y) = (Kx(x,y)2 +

Ky(x,y)2)1/2 of the dominant wave at each location is shown in Fig. 2d. Again, the different regimes of each wave in the specied wave eld are clearly distinguished.

The direction of wave propagation (x,y), measured anticlockwise from the x axis, is found as tan1 KxKy and plotted as Fig. 2e. Note that (x,y) is subject to a [notdef] radian am

biguity, which is reconciled with a priori information in our AIRS analysis in Sect. 5.3.

To assess the effectiveness of our spectral analysis of the specied wave eld, we compare the known wavelengths and propagation angles of the synthetic waves in the test image with the 2DST-measured wavelengths and propagation angles in Fig. 2f. For each wave numbered 18, blue dots show the input wavelength IN against measured wavelength OUT, indicating that the 2DST measures the horizontal wavelengths and propagation angles in the test image very well. Generally, shorter wavelengths are well resolved but longer wavelengths are slightly underestimated. This may be due to the coarse spectral resolution of DFT-based methods for waves with wavelengths that are a large fraction of the image size, since such waves can be more susceptible to spectral leakage problems.

4 An alternative spectral window

The use of the Gaussian window in the S-transform has some convenient mathematical advantages; it is analytically simple and has a denite integral over an innite range. However, when it is used for 2DST analysis an unfortunate side effect of the Gaussian window is the poor recovery of wave amplitude, discussed in the previous section. Although a Gaussian is traditionally used, any suitable apodizing function may be used, so long as its spatial integral is equal to unity (Stock-well, 2007). For example, Pinnegar and Mansinha (2003) used an asymmetric hyperbolic time-domain window for enhanced measurement of the onset times of one-dimensional time series components.

In this section, we introduce a new spectral windowing function for the 2DST. This new function takes the shape of an ellipse in the wavenumber domain, and a rst-order Bessel function of the rst kind J1(z) function with a scaled 1/z envelope in the spatial domain (for denition of z see Eq. (14) below). For this reason we refer to this window as the EllipticBessel window. We nd that when AIRS measurements are analysed with the 2DST using this new EllipticBessel window in place of the traditional Gaussian, the measurement of gravity wave amplitudes is greatly improved. Spectral resolution is also improved slightly, without adversely compromising spatial resolution.

4.1 The EllipticBessel window

As discussed in Sect. 3.1, the spectral peaks in a DFT spectrum have a characteristic width, where the spectral power is spread in a broad peak around the central wavenumber. This spectral power is slightly reduced when a Gaussian window is applied, due to the immediate decrease in the Gaussian function around the central location. The effect can be mitigated, but not fully reconciled, by decreasing the scaling parameter c, which broadens the Gaussian window in the wavenumber domain. However, this decreases the width of the spatial window, which increases the effect of wave undersampling for low wavenumbers.

One solution to this problem is to use a window that is an ellipse in the wavenumber domain. Here we introduce an EllipticBessel window Web, dened in the wavenumber domain as the ellipse

Web( x, y,kx,ky) =

8 >

<

>

:

0 for

x a

2+

y b

2[greaterorequalslant]1

,

1 for

x a

2+

y b

2 < 1

(12)

where a = [notdef]kx[notdef]/2c and b = [notdef]ky[notdef]/2c are the widths in

the x and y directions. We see that the semi-major and semi-minor axes of this voice ellipse scale with angular wavenumbers and are equal to the standard deviations of the equivalent voice Gaussian window in Eq. (11).

A key feature of this new window is that, in the wavenumber domain, it does not immediately decrease with displacement from the central location but rather has a scalable elliptical region within which the function is equal to unity. Thus, the window captures a much greater extent of the targeted spectral peak at kx and ky, which can greatly improve wave amplitude recovery compared to the traditional Gaussian window. The width of Web can also be more carefully adjusted using the scaling parameter c. This window can then be used in place of the Gaussian windowing term in Eq. (10).

One requirement for any apodizing window used in the Stockwell transform is that the spatial integral of the function must be equal to unity. This is so that the spatial integral of the Stockwell transform is equal to the Fourier transform H(kx,ky) (Mansinha et al., 1997a), namely

[integraldisplay]

1

1 1 S(x,y,kx,ky)dx dy = H(kx,ky), (13)

which has the useful result of making the 2DST fully invertible (Stockwell, 1999, 2007).

The normalisation term [notdef]kx[notdef][notdef]ky[notdef]/83c2 in the traditionally

used Gaussian window in Eq. (9) ensures that the Gaussian window satises this requirement. To check that the EllipticBessel function is admissible as an apodizing function, we must rst nd its spatial form, then check that its spatial integral is also equal to unity.

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

2552 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

Figure 3. Illustrative surface plots of the wavenumber-domain (top row) and spatial-domain (bottom row) forms of the traditional Gaussian (a, c) and new EllipticBessel (b, d) windowing functions used in the two-dimensional Stockwell transform for arbitrary wavenumbers kx and ky and scaling parameter c. The semi-major and semi-minor axes of the EllipticBessel window in (b) are equal to the corresponding standard deviations of the Gaussian window in (a), where both windows have a central value equal to unity. In the spatial domain, the EllipticBessel window in (d) has a central value equal to half that of the Gaussian window in (c), both of which have spatial integrals equal to unity. For details, see text.

The EllipticBessel window Web( x, y,kx,ky) is easily dened in the wavenumber domain as an ellipse, but its spatial form, which we denote as web(x,y,kx,ky), is given as

web(x,y,kx,ky) = [notdef]

kx[notdef][notdef]ky[notdef]

83c2

The surfaces in Fig. 3a and c show the wavenumber-domain and spatial-domain forms of the Gaussian window, for arbitrary wavenumbers kx and ky. As discussed above, the maximum value of the Gaussian is equal to unity in the wavenumber domain but equal to [notdef]kx[notdef][notdef]ky[notdef]/83c2 in the spa

tial domain such that its spatial integral is equal to unity.

This is a requirement of any windowing function in the S-transform. The standard deviations of the wg and Wg scale with wavenumbers kx and ky as described in Sect. 2.3 (Eqs. 9 and 11), providing the voice Gaussian.

Likewise, Fig. 3b and d show the wavenumber-domain (Web) and spatial-domain (web) forms of the EllipticBessel window, for the same arbitrary wavenumbers kx and ky as used for the Gaussian windows in panels a and c. The semi-major and semi-minor axes of the elliptic region in Fig. 3b are scaled with wavenumbers kx and ky and are equal to the standard deviations of the equivalent Gaussian in Fig. 3a, providing the voice ellipse.

The spatial-domain form of the EllipticBessel window web, described by a Bessel-shaped function within an envelope, is shown in Fig. 3d and described by Eq. (14). The maximum value of web is [notdef]kx[notdef][notdef]ky[notdef]/163c2, which ensures

that its spatial integral is equal to unity. This is equal to half of the maximum value of the equivalent Gaussian in Fig. 3c, since the terms involving z in Eq. (14) tend to 1/2 as x ! 0 and y ! 0. The width of the central region of web

is very slightly wider than the equivalent Gaussian, resulting in slightly coarser spatial resolution. However, the ability to

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

J1(z)

z , (14)

where J1 is the rst-order Bessel function of the rst kind (Abramowitz and Stegun, 1964) and

z =

1

2c

qk2xx2 + k2yy2.

A short derivation of the function in Eq. (14) is provided in Appendix A1. Fortunately, the spatial integral of Eq. (14) is indeed equal to unity, proof of which is presented in Appendix A2. This conrms that the EllipticBessel window is admissible as an apodizing function in the 2DST and validates its use in this study and beyond.

To recap our notation in this study, we have described two windowing functions for the 2DST: the traditional Gaussian and the new EllipticBessel windows, which we denote in the spatial domain as wg(x,y,kx,ky) and web(x,y,kx,ky) respectively and in the wavenumber domain as Wg( x, y,kx,ky) and Web( x, y,kx,ky) respectively, where the Wg and Web are the Fourier transforms of wg and web.

Figure 3 shows three-dimensional surface plots of the spatial and wavenumber domain forms of the traditional 2-D Gaussian and new EllipticBessel windows used in the 2DST here.

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2553

Figure 4. Orthographic projection of AIRS brightness temperature perturbations T [prime]

AIRS (a) from a granule over the Southern Andes at 05:30 UTC on 24 May 2008, with reconstructed temperature perturbations T [prime]

2DST

and horizontal wavelengths H computed using the 2DST with three different windowing approaches: (b, e) a Gaussian window with scaling parameter c = 1, (c, f) a Gaussian window with

c = 0.25 and (d, g) the EllipticBessel window with c = 0.25.

tune the 2DST with the scaling parameter c ensures that this effect can be compensated by a reasonable trade-off.

In the next section we show that the use of the Elliptic

Bessel window in the 2DST, in place of the traditional Gaussian window, signicantly improves wave amplitude recovery. This is very useful for our analysis of AIRS data in Sect. 5.

4.2 Invertibility

A very convenient aspect of the S-transform is its invertibility. Since we have shown here that both the traditional Gaussian and new EllipticBessel windows have spatial integrals equal to unity, the 2DST can be completely inverted to recover the original 2-D image, whichever of these windows or real non-zero positive values of the scaling parameter c are used. Note that a traditional 1-D or 2-D CWT does not necessarily have this capability. The fact that we can achieve such exibility in spatialspectral resolutions by swapping windows or by adjusting c, yet still retain the capability of inversion, further highlights the strength of the 2DST as a tool for spatialspectral analysis of geophysical data.

Unfortunately, to take full advantage of DFT algorithms and the inversion capability of the 2DST for AIRS data, we must compute the 2DST using all permitted wavenumber voices in both dimensions. This requires nearly 12 000 inverse DFT calculations for each AIRS granule using the traditional voice-by-voice implementation described here, the computational load of which could be quite impractical for

large-scale studies. Interpolating AIRS measurements to a coarser resolution with fewer pixels could be one solution to reduce computational cost, but this will obviously undersample short horizontal wavelengths in the data. Faster methods for computing the S-transform have been developed (Brown et al., 2010) which may increase practicality in the future. Other steps, such as avoiding programming loops and ensuring that any 2-D objects to be transformed have dimensions that are powers of two, may also reduce relative computational expense.

4.3 The effect of window choices on AIRS analysis

Figure 4 shows an AIRS granule over the Southern Andes measured on 24 May 2008, analysed using the 2DST with three different windowing approaches.

In Fig. 4b and e, we use a Gaussian windowing function with the scaling parameter c set to 1. This is the window usually used in 1DST implementations. We see that, as discussed above, this choice of window is only able to recover the very general, long-horizontal wavelength features of the granule, with poor spatial localisation and signicantly reduced amplitude. This is due to a large proportion of the spectral response being lost by the windowing Gaussian when applied to two dimensions.

We can reduce the impact of this by decreasing the scaling parameter c, which broadens (narrows) the spectral (spatial) window. This provides improved amplitude recovery and improved spatial localisation at the expense of spectral localisa-

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

2554 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

tion (Fritts et al., 1998). Since we only select a single dominant spectral peak for each location on the granule, this is acceptable for our purposes. The reconstructed perturbations and horizontal wavelengths (Fig. 4c and f) are now much more representative of the wave features in the granule.

One problem remains, however. By decreasing c, we narrow our spatial window. In regions where wave amplitudes are low, such as the bottom-left corner of Fig. 4a, this narrow Gaussian window starts to undersample long wavelengths, such that only very short wavelengths are attributed to the region. The EllipticBessel window used in Fig. 4d and g performs better at recovering the underlying larger-scale structure of the granule, without defaulting to the small-scale noisy variations. Amplitude recovery at all wavelengths is also improved over either of the Gaussian approaches.

In the general case, these low-amplitude, small-scale variations are unlikely to be due to gravity waves with vertical wavelengths visible to AIRS, so their recovery is something we try to avoid. Furthermore, such wavelengths are very close to or at the Nyquist limit for these data. Our condence in their measurement is thus very low, yet the momentum uxes they transport can dominate. We discuss this further in Sect. 5.4.

For the windowing functions considered, it is clear from Fig. 4 that the scaling parameter c has a signicant effect in determining the spatialspectral localisation capabilities of the 2DST. The EllipticBessel windowing function, with a scaling parameter of c = 0.25, was selected for our AIRS

analysis in the next section. This choice provided the best trade-off between spatial and spectral localisation of different wave regimes in AIRS measurements.

As discussed in Sect. 3.1, the 2-D images (x,y), Kx(x,y) and Ky(x,y) contain the dominant measured wave amplitudes and wavelengths at each location on the granule. These images are computed on a pixel-by-pixel basis, selecting a single monochromatic wave with the largest amplitude in the localised spectrum for each pixel.

As a result, the reconstructed images shown in Fig. 4b d, computed by taking the real part of the complex image (x,y), will never be perfect representations of the input data, but they will provide a best guess of the dominant features of the granule.

Since we have shown that the 2DST is fully invertible for both the Gaussian and EllipticBessel windowing approaches (Sect. 4 and Appendix A), a complete reconstruction of the input image is of course producible by taking the inverse of the full 4-D 2DST object, but here we desire 2-D maps of wave properties, so a best guess method is used.

A possible quantitative metric to assess the rst-order effectiveness of our 2DST analysis in Fig. 4 could be to compare the variance of the input image with the variances of each of the reconstructions. However, since the reconstructions are computed as a best guess method on a pixel-by-pixel basis, their total variance is not readily related to the total variance of the input image and thus may not be mean-

ingful as a comparison. Furthermore, such use of the image variance would only be appropriate if the distribution of perturbations was unimodal and ideally Gaussian, which is not the case for an image of a sinusoidal wave. In practice, however, we generally expect the variance of the reconstruction not to exceed the variance of the input image, since wave amplitudes computed on a pixel-by-pixel basis from a localised spectrum will usually be underestimated for the reasons given in Sect. 3.1.

It is not impossible that in some rare cases the total variance of the reconstruction could exceed the total variance of the input image, for example due to the spatial extent of a wave feature being slightly overestimated. If the localised spectrum for one pixel is affected a larger-amplitude wave feature in one of its neighbouring pixels, this can result in subtle articial borders between different wave regimes in the reconstructions. This is not a limitation of the 2DST itself, but arises in the somewhat forced extraction of localised gravity wave parameters contained in the 4-D Stock-well transform object S(x,y,kx,ky) in order to produce the 2-D image. This effect should be carefully considered in future work to ensure wave properties are not over-represented.

5 AIRS gravity wave analysis using the 2DST

In this section, we use our 2DST-based method to perform gravity wave analysis on two-dimensional granules of AIRS radiance measurements, comparing our analysis to that of previous studies. We use the 2DST to measure gravity wave amplitudes, horizontal wavelengths, and directions of propagation. We then use ECMWF-derived wind speeds and the assumption of an orographic wave source to infer vertical wavelengths and make estimates of gravity wave momentum ux (the vertical ux of horizontal pseudomomentum) by closely following the method of Alexander et al. (2009).

5.1 AIRS granule selection and pre-processing

The rst AIRS granule selected for our study is granule 32 of 6 September 2003, over South Georgia. The second granule is a 135-pixel swath over the intersection between granules 39 and 40 on 2 August 2010, located over the Antarctic Peninsula and Drake Passage.

Alexander et al. (2009) and Hoffmann et al. (2014) performed an analysis of these AIRS granules over South Georgia and the Antarctic Peninsula respectively. Their studies measured wave amplitudes, horizontal wavelengths and wave propagation directions using a 1-D S-transform method, as described by Alexander and Barnet (2007). In their method, the 1-D S-transform is computed for each cross-track row. Then, covariance spectra are computed between pairs of adjacent cross-track rows to measure phase shifts in the along-track direction, from which along-track wavelengths can be inferred. To nd the dominant wave fea-

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2555

tures in a granule, these co-spectra are averaged together and up to ve spectral peaks are found in this averaged spectrum.This approach can provide computationally fast, rst-order gravity wave analysis of AIRS granules, and it has been used in numerous other studies (e.g. Alexander and Teitelbaum, 2011; Alexander and Grimsdell, 2013; Wright et al., 2016a).

One limitation of this method is that the phase difference measurements required to recover along-track wavenumbers can introduce a strong cross-track bias in resolved features, since the S-transform is only computed in the cross-track direction. In addition, waves which occupy only small regions of the granule in the along-track direction may also be under-represented in the averaged co-spectrum. Furthermore, selecting no more than ve dominant waves in the averaged co-spectrum implicitly limits the maximum number of available along-track wavenumber voices to no more than ve for each location on the entire granule. The use of a 2DST is a logical solution to each of these problems. With the increased convenience of computational power since the study of Alexander and Barnet (2007), the 2DST now represents a more practical alternative to the one-dimensional method.

Before implementing the 2DST, each granule of brightness temperature perturbations is interpolated onto a regularly spaced grid with approximately 17.7 and 20.3 km separating adjacent pixels in the along-track and cross-track directions respectively. In the centre of the AIRS swath, the resolution of this regularly spaced grid closely matches the spatial resolution of AIRS, so very little if any information is lost. Towards the edge of the swath, this grid is ner than the spatial resolution than AIRS, but the grid points will not exactly match the location of the AIRS footprints. A useful graphic of typical AIRS footprints can be found in Hoffmann et al. (2014, their Fig. 2).

As a result of using DFT algorithms, the maximum numbers of permitted wavenumber voices available in the along-track and cross-track directions are limited to NAT 1 and

NXT 1, where NAT and NXT are the number of pixels in

the along-track and cross-track directions (135 and 90 respectively). These wavenumber voices have corresponding wavelengths that are integer fractions of the total along-track and cross-track dimensions of the granule. Here, we compute the 2DST for wavelengths greater than around 40 km. This is just over twice the Nyquist-sampling distance between AIRS pixels after interpolation onto our regular grid. The zeroth frequencies are omitted. Increased along-track spectral resolution at low wavenumbers can be obtained by applying the 2DST to two or more adjacent granules, thus increasing the number of along-track pixels.

5.2 AIRS gravity wave properties measured by the 2DST

The results of our 2DST analysis of the selected AIRS granules over South Georgia and the Antarctic Peninsula are shown in Figs. 5 and 6 respectively.

In both Figures, panel a shows the brightness temperature perturbation measurements calculated as described in Sect. 2.1. Note that the colour scale is chosen so as to make wave perturbations clearer by eye, but at some locations it is saturated.

Clear wave-like perturbations are observed in both granules directly over and to the east of the mountain ranges. As in previous work, such clear wave-like perturbations are attributable to gravity waves with a high degree of certainty.

Reconstructed 2DST temperature perturbations T [prime]2DST

are shown in panel b. These are found by taking the real part of the complex 2DST object (x,y) as described in Sect. 3.1.

The image T [prime]2DST shows the dominant wave features in

the granule reconstructed using only the pre-dened range of permitted wavenumber voices in the 2DST. Since we only consider the coefcients of the dominant wavenumber at each location, this reconstruction cannot be perfect, but it provides a visual inspection of how well the 2DST outputs represent the dominant wave characteristics of the granule.The T [prime]2DST image can be used to ne-tune the 2DST by

changing the windowing function, by adjusting the scaling parameter c or by redening the range of frequency voices until the desired outcome is achieved. Such ne-tuning exibility cannot be so easily achieved using the 1DST method.

Generally, the agreement between reconstructed wave features in Fig. 5b and AIRS measurements in Fig. 5a is very good, but some uncertainties remain. As discussed in Sect. 4.3, there is some discrepancy regarding the spatial extent of some wave features, such as a small positive wave crest located just south-west of South Georgia which appears to be located slightly east, with an apparently slightly overestimated amplitude, than is observed in the AIRS measurements. Conversely, at 63 S, 28 W, a positive bow-shaped wave crest is observed in the AIRS measurements but is underestimated in the reconstruction. As mentioned in Sect. 4.3, these small misrepresentations are not a limitation of the 2DST itself but rather the forced extraction of gravity wave parameters from the 4-D S-transform object in order to create the 2-D reconstruction, where only one single wave feature with the largest localised spectral amplitude is assigned at each location. The overall agreement is still very good, but future work to improve the extraction of gravity wave parameters from the 4-D S-transform object may help to resolve some of these discrepancies.

Figures 5c and 6c show full underlying wave amplitudes

|T [prime][notdef]2DST for each granule. This is found by taking the abso

lute magnitude of the complex 2DST object (x,y) as de

scribed in Sect. 3.1. This property provides us with a useful metric with which to dene regions of the granule which do or do not contain wave-like perturbations, such that we can limit spurious detections (e.g. McDonald, 2012). In Figs. 5h and 6ch, we exclude regions of each granule where the underlying wave amplitude is more than 1 standard deviation below the mean underlying wave amplitude of the granule.

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

2556 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

Figure 5. Orthographic projections of a granule of AIRS brightness temperature measurements T [prime]

AIRS (a) over South Georgia at around 03:00 UTC on 6 September 2003 and selected outputs (be) of our 2DST analysis using the EllipticBessel window. This granule was also analysed by Alexander et al. (2009) (their Fig. 3) using a one-dimensional S-transform method. The 2DST outputs shown here are reconstructed brightness temperature perturbations T [prime]

2DST (b)

, underlying wave amplitudes [notdef]T [prime][notdef]2DST (c), horizontal wavelengths H (d)

and wave propagation directions (e) in degrees anticlockwise from the east. Also shown are mean wind speed parallel to the horizontal wavenumber vectors U[notdef][notdef] (f) from ECMWF operational analyses at z 40 km, vertical wavelengths Z (g) and the magnitude of the horizontal

component of vertical momentum ux [notdef]MF[notdef]x,y (h). Black arrows in (f) and (h) show the horizontal direction of U[notdef][notdef] and [notdef]MF[notdef]x,y respectively.

For details, see text.

In Fig. 5cg, we do not exclude such regions for discussion purposes, so as to provide an example of the data we would otherwise omit.

Figures 5d and 6d show absolute horizontal wavelengths H = 2(k2AT + k2XT)1/2, where kAT and kXT are the along-

track and cross-track angular wavenumbers respectively. We can see that these horizontal wavelengths clearly dene different regimes of the dominant wave features of the granules, as in the test case in Sect. 3.1, though the AIRS data are more complex. In the South Georgia granule in Fig. 5d, we see that the island lies within a wave eld where long horizontal wavelengths are dominant around and to the east of the island over the ocean, with their wavenumber vectors aligned roughly parallel to the direction of the mean ow.This is characteristic of a wing-shaped mountain wave eld (Alexander and Grimsdell, 2013) and is in good agreement with visual inspection of the granule itself.

In Figs. 5e and 6e, we show the orientation of the horizontal wavenumber vector measured anticlockwise from east. is calculated by rst projecting the along-track and cross-track wavenumber vectors kx and ky into their zonal and

meridional components k and l using the azimuths of the along-track and cross-track directions at each location on the granule, then taking = tan1

l k

. Note that only describes the orientation and not the true horizontal direction of propagation of the wavenumber vectors, which retain a

180 ambiguity that we break below.

In the South Georgia granule (Fig. 5), we see that our 2DST measurements in the southern region of the granule are largely dominated by small-scale, low-amplitude, short horizontal wavelength features with random directions of propagation. Most of these features are likely to be due to noise and not attributable to coherent wave structures. By using a threshold amplitude, such regions are effectively removed, leaving well-dened regions with clear wavelike perturbations. The contribution of small-scale features that remain after this step is discussed further in Sect. 5.4.

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2557

Figure 6. As Fig. 5 but for AIRS measurements over the Antarctic Peninsula around 04:00 UTC on 2 August 2010. These measurements were also analysed by Hoffmann et al. (2014, their Fig. 8) using a one-dimensional S-transform method. Here, regions in (ch) where the underlying wave amplitude [notdef]T [prime][notdef]2DST is less than 1 standard deviation below the mean are coloured white.

5.3 Momentum uxes

Here we make estimates of gravity wave momentum ux for the dominant wave-like features measured by the 2DST in our selected granules, following the method of Alexander et al. (2009).

Ern et al. (2004) showed that the zonal and meridional components of gravity wave momentum ux MFx and MFy can be given by

MFx,MFy

[parenrightbig]

=

2

[parenleftBig]

g N

2 T [prime]a T

2[parenleftbigg]

k

m,

l m

, (15)

where is density at a height of 40 km, g is the acceleration due to gravity, N is the buoyancy frequency, T [prime]a is the attenuation-scaled full underlying wave amplitude, T is the background temperature and k, l and m are wavenumbers in the zonal, meridional and vertical directions respectively. Equation (15) assumes the mid-frequency approximation (Fritts and Alexander, 2003, their Sect. 2.1), which is valid for a large portion of the observable gravity wave spectrum. For these waves, m2 (k2 + l2) and Coriolis ef

fects are negligible, as is the case for nearly all waves measured here due to the size of the AIRS beam footprint (Hoff-mann et al., 2014) and the vertical weighting function of the 667.77 cm1 channel (Alexander and Barnet, 2007). When compared with momentum uxes calculated without mak-

ing the mid-frequency approximation, Ern et al. (2004) found discrepancies typically not exceeding around 10 %.

We also have a [notdef]180 ambiguity in direction of propaga

tion, which we break by assuming the waves in our granules always propagate against the mean ow. One last variable in Eq. (15) is vertical wavenumber m = 2/ Z. Under the as

sumption that the waves are upwardly propagating mountain waves with ground-based phase velocity equal to 0, Z is

given as

Z

2U[notdef][notdef]

N , (16)

where U[notdef][notdef] is the component of the mean wind speed parallel

to the waves horizontal wavenumber vector (Eckermann and Preusse, 1999). U[notdef][notdef] is found by projecting the mean wind

vector U in the direction of the wave propagation angle shown in Figs. 5e and 6e.

Figures 5f and 6f show U[notdef][notdef] coincident with each granule at

an altitude of 40 km from ECMWF operational analyses, projected onto each granules regular grid. Orange contours and black arrows show the magnitude and direction of U[notdef][notdef] respec

tively. Vertical wavelength Z is shown in Figs. 5g and 6g.

Towards the south-eastern corner of both granules, mean wind speeds become quite weak. As they fall below around 40 m s1, vertical wavelengths start to drop below the vertical resolution limit of the AIRS channel. The wave eld

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

2558 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

may continue into this region, but the vertical wavelengths may be too short to be resolved such that wave amplitudes are attenuated to below the ambient noise level. This is particularly clear in the Antarctic Peninsula granule, where detectable wave fronts abruptly terminate just as U[notdef][notdef] begins to

fall below 30 m s1.

No further wave amplitude attenuation corrections were applied to regions where vertical wavelengths are inferred to be below 12 km, since this is beyond the resolution limit of the weighting function of the 667.77 cm1 channel. It should be noted, however, that although theoretically valid, the attenuation correction applied here can in some cases increase wave amplitudes by perhaps an order of magnitude or more, despite being signicantly susceptible to errors in wind speed, wave direction or both. This correction should be applied with caution, and the effects of its inclusion are discussed below and in Sect. 5.4.

In the South Georgia granule, peak momentum ux values of more than 500 mPa are associated with a small region of large amplitude and short horizontal wavelength wave features, located toward the south-eastern tip of the island. In the Antarctic Peninsula granule, momentum uxes of a few hundred millipascals are generally co-located with the clearly visible wave structures in the raw brightness temperature perturbations just downwind of the peninsula.

The key strength of the results presented here is the much-improved spatialspectral localisation and resolution capabilities provided by full two-dimensional treatment of the AIRS data. Condence in the accuracy of subsequent measured quantities in our 2DST-based analysis is thus greatly improved over previous 1-D S-transform-based methods.Understandably, the former is more computationally intensive than the latter, and this should be considered if data sets are large or computational resources are limited.

Although the magnitude, direction and distribution of momentum uxes in both granules are broadly in line with previous AIRS gravity wave studies in the region (e.g Alexander and Teitelbaum, 2007; Alexander et al., 2009; Alexander and Teitelbaum, 2011), momentum ux magnitudes are much higher than those observed in studies using limb-sounder data sets (e.g. Ern et al., 2004, 2011; Alexander et al., 2008;Geller et al., 2013; Hindley et al., 2015; Wright et al., 2016a), which generally do not exceed a few tens of mPa in this region. While we would not expect the magnitude of these uxes to be the same as those observed by limb sounders (rstly, nadir-sounding instruments are generally more sensitive to waves with longer vertical wavelengths and hence higher momentum uxes due to the deep vertical weighting function (Alexander and Barnet, 2007); secondly, momentum ux estimates from limb sounders are typically lower-bound estimates due to the projection of horizontal wavelengths (Ern et al., 2004); and, thirdly, wavelength-dependent wave amplitude attenuation corrections are not generally applied to limb-sounder results), we note that our results and those of other AIRS gravity wave studies in this region,

which use a correction factor for wave amplitude attenuation based upon a vertical wavelength estimation, are substantially higher. This effect and its implications are discussed in the next section.

5.4 Small-scale perturbations and the attenuation correction

In the implementation of any spectral image processing, it is important to strike a balance between accurate measurement of the desired properties and the spurious interpretation of noise. One of the advantages of AIRS measurements is the high horizontal resolution of the data. With the close exception of the Infrared Atmospheric Sounding Interferometer (IASI; e.g. Clerbaux et al., 2009), currently no other space-borne instrument can measure stratospheric gravity waves with comparable horizontal resolution to AIRS. Therefore, accurate measurement of resolved waves with horizontal wavelengths close to the AIRS resolution limit are of great importance. Such short horizontal wavelength waves, if reliably resolved, will generally carry higher momentum uxes via Eq. (15). However, as we approach the resolution limits of AIRS measurements and our spectral methods, our condence in the accuracy of our measurement of such waves decreases.

In Fig. 5, our 2DST analysis resolves a very small region of short horizontal wavelengths over the south-eastern tip of South Georgia. Wave perturbations in this small region are just a few pixels across, but their uncorrected brightness temperature perturbations are large, with peaks of order 56 K.If these perturbations were located almost anywhere else on the granule we would likely attribute them to retrieval noise.

However, gravity waves with very short horizontal wavelengths, tightly packed in a region immediately downwind of a mountainous island, are in good agreement with mountain wave theory. Examples of such waves can be found in, for example, the modelling studies of Shutts and Vosper (2011) and Alexander and Teitelbaum (2011). If these waves are indeed real, their accurate measurement is of great importance.

An added complication is introduced as a result of the mountain wave assumption used. Since these waves are only a few pixels across, their directions of propagation are difcult to dene, introducing a random element. The component of the mean wind parallel to the horizontal wavenumber vector can thus be very low, which decreases the vertical wavelength estimate, which in turn increases the attenuation correction applied to the observed temperature perturbations.This attenuation correction can increase temperature perturbations by 400 % or more which, since momentum ux is proportional to the square of wave amplitude, can increase our estimate of momentum ux to extremely large values.

Thus, by applying a correction factor for wave amplitude attenuation, very small-scale low-condence perturbations can yield extremely high momentum uxes which can dominate the momentum budget of the entire granule if the moun-

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2559

As is expected, these values are much lower than those calculated using the attenuation-corrected wave amplitudes.

If these extremely high uxes correspond to real waves, then their measurement is of crucial importance. However, if such perturbations are simply instrument noise and their uxes are spurious, then the biases and errors introduced by their inclusion in broader studies could be very large.

Without a priori knowledge of the wave environment, which is most readily gained by visual inspection of the AIRS measurements, it would be unwise to include the uxes from these small-scale perturbations in any automated analysis. Indeed, Wright et al. (2015) suggest that when AIRS granules are pre-smoothed with a boxcar of width 3 pixels, resolved momentum uxes calculated using the method of Alexander et al. (2009) can be reduced by a order of magnitude. This suggests that these small-scale, high-momentum ux features are reasonably common and can impact larger-scale momentum ux estimates, either realistically or spuriously.

Further work investigating this problem is encouraged.

Pre-smoothing granules so as to exclude these perturbations (Wright et al., 2015, 2016b) is one solution. Excluding even more wavenumber voices corresponding to short horizontal wavelengths from our 2DST analysis is another solution.However, in cases such as the South Georgia granule presented in Fig. 5, where such perturbations may well be physical, this exclusion can reduce peak momentum uxes by orders of magnitude, introducing a systematic low bias and thus further uncertainty.

6 Summary and conclusions

In this study, we have applied the 2DST to granules of AIRS measurements, extracting gravity wave amplitudes, wavelengths and directions of propagation. Our 2DST method builds upon the work of Alexander and Barnet (2007), who used the 1-D Stockwell transform for the same purpose.Their method can introduce a strong cross-track bias problem, which we solve by using a full 2DST.

We rst dene our 2DST implementation and test it on a specied wave eld containing synthetic waves with known amplitudes, wavelengths and directions of propagation. We nd that the 2DST provides very good spatial representation of the dominant spectral components of the specied wave eld, accurately measuring wavelengths and orientations of all the synthetic waves.

Due to the spread of spectral power in the spectral domain and wave undersampling in the spatial domain, we nd that localised wave amplitudes as measured by the 2DST are reduced by more than a factor of two when the typical Gaussian windowing function is used in the Stockwell transform.We compensate for this by decreasing the scaling parameter c and by replacing the Gaussian window with a new alternative EllipticBessel window, which we test on a granule of

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

Figure 7. As Figs. 5h and 6h but for estimates of the magnitude of gravity wave momentum ux over South Georgia (a) and the Antarctic Peninsula (b) calculated without correcting for attenuation of wave amplitude as discussed in Sect. 5. Momentum uxes less than 0.5 mPa are coloured white.

tain wave assumption is used, correctly or otherwise. This is evident in the South Georgia granule studied by Alexander et al. (2009) and shown in our Fig. 5, where peak momentum uxes are localised over only a few large-amplitude pixels just over the south-eastern tip of the island. In Fig. 5h, the ux peaks at almost 1000 mPa for one pixel in this small region, though the colour scale is saturated. Whether these uxes are real or not, they can nevertheless be extremely large and should be approached with caution.

Figure 7a and b show gravity wave momentum uxes over the South Georgia and Antarctic Peninsula granules as shown in Figs. 5 and 6, only with uxes calculated without scaling the temperature perturbations for attenuation as described in Sect. 5. The highest momentum ux values in Fig. 7a and b are around 70 and 30 mPa respectively, located within small regions in the wave eld where the colour scale is saturated.

2560 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

AIRS measurements over the southern Andes. We nd that this new spectral window provides a better balance between spatialspectral localisation and the accurate measurement of wave amplitudes. Wave amplitude recovery is thus improved to around 80 to 90 % of input values. We also demonstrate that this new window is indeed valid for use in the 2DST by showing its spatial integral is equal to unity.

Next, we measure gravity wave amplitudes, horizontal wavelengths and directions of propagation in two granules of AIRS measurements over South Georgia and the Drake Passage/Antarctic Peninsula region. Our 2DST method signicantly improves two-dimensional representation of the dominant spectral features of the granules over previous 1DST methods. These spectral features are directly measured in both dimensions simultaneously for each location of the granule, without the introduction of potential biases caused by the use of averaged co-spectra. This is a clear advantage over previous methods.

Another key advantage of our 2DST method is the ability to visually inspect the quality of our spectral analysis. By taking the real parts of the dominant localised spectral coefcients at each location, a reconstruction of the granule can be created. This can be used to ne-tune the adjustable parameters and provide a useful sanity check on the performance of the 2DST. Future work may involve comparing this output to the original data via a variance argument or similar, such that we can obtain a quantitative measure of the quality of the 2DST analysis for quality control purposes in larger-scale studies.

To conclude, our new 2DST-based gravity wave analysis method for AIRS data makes signicant improvements over current methods in several key areas, and we would advocate its use in future work.

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2561

Appendix A: Admissibility of the EllipticBessel window in the Stockwell transform

In Sect. 4 we introduced the EllipticBessel window as a new apodizing function for the 2DST. One requirement for any apodizing function for use in the Stockwell transform is that its spatial sum must be equal to unity. If this condition is satised, the spatial sum of the 2DST is equal to the 2-D Fourier transform, making the 2DST fully invertible.

In this appendix we demonstrate that the EllipticBessel window is admissible as an apodizing function in the S-transform. To do this, we must rst nd the spatial analogue of the wavenumber-domain ellipse we dened in Eq. (12).We must then take the spatial integral of this function to demonstrate that it is equal to unity.

A1 The EllipticBessel window in the spatial domain

The EllipticBessel window is dened in the wavenumber ( x, y) domain as

Web( x, y,kx,ky) =

as

ei( xx+ yy)

=(cos( xx) + i sin( xx)) cos( yy) + i sin( yy) [parenrightbig]

=cos( xx)cos( yy) + i sin( xx)cos( yy)

+i sin( yy)cos( xx) sin( xx)sin( yy). (A4)

We can omit the last three sin() terms in Eq. (A4) since, due to the symmetry of the sine function around (0,0), each term will eventually sum to 0. We can then rewrite Eq. (A3) as

web(x,y,kx,ky) =

4 42

rb2 2x b2 a2

[integraldisplay]0

a

[integraldisplay]

0

8 >

<

>

:

0 for

x a

2+

y b

2[greaterorequalslant]1

,

1 for

x a

2+

y b

2 < 1

cos( xx)cos( yy)d x d y. (A5)

This integral can be further simplied if we switch to polar coordinates using the substitutions x = ar cos() and y =

br sin() after which the expression in Eq. (A5) becomes

web(x,y,kx,ky) =

ab 42

2

[integraldisplay]

0

(A1)

where a = [notdef]kx[notdef]/2c and b = [notdef]ky[notdef]/2c are the half-widths of

the ellipse in the x and y directions (see Fig. 3b). The spatial-domain form of the EllipticBessel window, denoted here by web(x,y,kx,ky), is found by taking the inverse 2-D

Fourier transform of Eq. (A1) as

web(x,y,kx,ky) (A2)

= F1xF1y

[integraldisplay]

0

1

cos(arx cos())cos(bry sin())r dr d. (A6)

Next we substitute A = arx and B = bry and, using mul

tiple angle formulae, rewrite Eq. (A6) as

web(x,y,kx,ky)

=

Web( x, y,kx,ky)

[bracketrightbig]

1

=

1 42

ab 42

[integraldisplay]

0

r

2

[integraldisplay]

0

1 1 Web( x, y,kx,ky)ei( xx+ yy) d x d y.

Since Web( x, y,kx,ky) = 1 within the ellipse and 0 ev

erywhere else, and has double symmetry, we can change the limits of integration to be the boundaries of the ellipse, expressing the total integral as a sum of four equal quadrants:

web(x,y,kx,ky) =

4 42

rb2 2x b2 a2

[integraldisplay]0

cos(Acos())cos(B sin()) d dr

1

=

ab 82

[cos(Acos + B sin)

+cos(Acos B sin)] d dr

[integraldisplay]

0

r

2

[integraldisplay]

0

1

=

ab 82

[integraldisplay]

0

r

2

[integraldisplay]

0

[cos(

pA2 + B2 cos( [Lambda1])

pA2 + B2 cos( + [Lambda1])]d dr, (A7)

where [Lambda1] = tan1(B/A). Here, [Lambda1] is simply an arbitrary

phase due to the periodicity of the cosine function when integrated over 0 to 2, so the integrals of both terms in the square brackets in Eq. (A7) will be equal. Hence we can sim-

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

a

[integraldisplay]

0

+cos(

ei( xx+ yy) d x d y. (A3)

We then recognise that the exponential term in the transform above can be replaced with sine and cosine functions

2562 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

ply add these terms such that we have

web(x,y,kx,ky) =

ab 82

1

d dr. (A8)

Next we recall the integral denition of the zeroth-order Bessel function of the rst kind J0(x) (Abramowitz and Stegun, 1964) given as

J0(x) =

1

2

This spatial-domain form of the EllipticBessel window in Eq. (A14) is plotted in Fig. 3d.

Equation (A14) describes a J1(z) function within a scaled 1/z envelope. Because of this, the terms involving z in Eq. (A14) converge to 1/2 as x ! 0 and y ! 0, such

that the central region of the function has peak value of

|kx[notdef][notdef]ky[notdef]/ 163c2, as shown in Fig. 3d. Interestingly, this

value is equal to half the peak value of the equivalent Gaussian window shown in Fig. 3c. The central peak of the EllipticBessel window is also, for each frequency voice, slightly broader than that of the equivalent voice Gaussian.

A2 Spatial integral of the EllipticBessel window

Now that we have found an analytical expression for the spatial-domain form of the EllipticBessel window (Eq. A14), we can proceed to check that it is admissible as an apodizing function in the 2-D Stockwell transform namely that its spatial sum is equal to unity (e.g. Pinnegar and Mansinha, 2003). The spatial sum of Eq. (A14), denoted here by I, can be written as

I =

[integraldisplay]

0

r

2

[integraldisplay]

0

2cos

[parenleftBig][radicalbig]

A2 + B2 cos[parenrightBig]

2

[integraldisplay]

0

cos(x cos) d (A9)

and substitute into Eq. (A8) and reintroduce our substitutions of A = arx and B = bry to give

web(x,y,kx,ky) =

1

qa2x2 + b2y2[parenrightbigg]

r dr. (A10)

We now use a new substitution that = r

pa2x2 + b2y2

and rewrite Eq. (A10) as

ab 2

[integraldisplay]

0

J0

r

1 1 |kx[notdef][notdef]ky[notdef]

83c2

web(x,y,kx,ky) =

ab 2

pa2x2+b 2y2

[integraldisplay]0

J0() a2x2 + b2y2

d. (A11)

Next we use the standard result (e.g. Abramowitz and Stegun, 1964) that

x2

[integraldisplay]

x1

J1(z)

z dx dy. (A15)

This integral can be simplied if we reintroduce our substitutions a = [notdef]kx[notdef]/ 2c and b = [notdef]ky[notdef]/ 2c and switch to

polar coordinates, using the substitutions x = [Lambda1]cos(')a and y = [Lambda1]sin(')b to give

I =

ab 2

2

[integraldisplay]

0

1

[integraldisplay]

0

J1([Lambda1]) [Lambda1]

[Lambda1]

ab d[Lambda1] d'

x J0(x)dx = xJ1(x) [vextendsingle][vextendsingle][vextendsingle]

x2

x1 (A12)

to rewrite Eq. (A11) as

web(x,y,kx,ky) =

=

1

2

2

[integraldisplay]

0

1

[integraldisplay]

0

J1([Lambda1])d[Lambda1] d'

ab 2

J1() a2x2 + b2y2

[vextendsingle][vextendsingle][vextendsingle]

pa2x2+b 2y2

0

=

2 2

1

[integraldisplay]

0

ab 2

J1([Lambda1])d[Lambda1]. (A16)

Using the standard result (e.g. Abramowitz and Stegun, 1964) that

x2

[integraldisplay]

x1

[parenrightbig]

a2x2 + b2y2[parenrightBig]

pa2x2 + b2y2. (A13)

Finally, recalling that a = [notdef]kx[notdef]/2c and b = [notdef]ky[notdef]/2c are

the half-widths of the original ellipse in Eq. (A1), we now can write the analytical expression for the spatial form of the EllipticBessel window as

web(x,y,kx,ky) = [notdef]

kx[notdef][notdef]ky[notdef]

83c2

J1

[parenleftBig][radicalbig]

=

J1(x)dx = J0(x)

[vextendsingle][vextendsingle][vextendsingle]

x2

x1, (A17)

we see that Eq. (A16) becomes

I = J0([Lambda1])

[vextendsingle][vextendsingle][vextendsingle]

1 0

J1(z)

z , (A14)

= (0) (1)

= 1 (A18) as required. This result conrms that the spatial sum of web(x,y,kx,ky) is indeed equal to unity; thus the Elliptic

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

where

z =

1

2c

qk2xx2 + k2yy2.

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2563

Bessel window is admissible as an apodizing window for the 2-D Stockwell transform.

A3 Admissibility of other windows

In this appendix so far, we have found a useful analytical expression for spatial form of the EllipticBessel window presented in this study. We have then shown that its spatial integral is equal to unity and it is thus admissible as an apodizing function in the 2-D Stockwell transform. In other cases, a quick test may be performed on candidate S-transform windowing functions to check whether their spatial integral is unity.

If we take the spatial integral I of the spatial-domain form of a candidate windowing function w(x,y,kx,ky), namely

I =

1

[integraldisplay]

1

1

[integraldisplay]

1

[vextendsingle][vextendsingle][vextendsingle]k

1 1 w(x,y,kx,ky)ei(kxx+kyy)dx dy

x =ky=0

(A20)

= W( x, y,kx,ky)

[vextendsingle][vextendsingle][vextendsingle]k

(A21)

where W( x, y,kx,ky) is the wavenumber domain form of the candidate window and the notation [notdef]kx=ky=0 denotes that

the function is evaluated at kx = ky = 0. This means that if

the value of W( x, y,kx,ky) evaluated at kx = 0 and ky = 0

is equal to unity, then its spatial integral will also be equal to unity. If it is not, then the candidate window is not admissible for use the Stockwell transform. Figure 3a and b show that both the Gaussian window and the EllipticBessel window are equal to unity at kx = ky = 0 and thus satisfy this

requirement. This short test may be helpful in the design of proposed alternative S-transform windowing functions in the future.

x =ky=0,

w(x,y,kx,ky)dx dy, (A19)

and introduce the factor ei(kxx+kyy), noting that when kx =

ky = 0 this factor is equal to unity, then I can be written as

I =

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016

2564 N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis

Acknowledgements. Neil P. Hindle is funded by a NERC studentship awarded to the University of Bath. Corwin J. Wright and Nicholas J. Mitchell are supported by NERC grant NE/K015117/1.The authors would like to thank the AIRS programme team for many years of hard work producing the data used here and also the handling editor and anonymous reviewers for their very helpful suggestions.

Edited by: L. Hoffmann

References

Abramowitz, M. and Stegun, I.: Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Applied mathematics series, Dover Publications, 1964.

Alexander, M. J. and Barnet, C.: Using satellite observations to constrain parameterizations of gravity wave effects for global models, J. Atmos. Sci., 64, 16521665, doi:http://dx.doi.org/10.1175/JAS3897.1

Web End =10.1175/JAS3897.1 http://dx.doi.org/10.1175/JAS3897.1

Web End = , 2007.

Alexander, M. J. and Grimsdell, A. W.: Seasonal cycle of orographic gravity wave occurrence above small islands in the Southern Hemisphere: Implications for effects on the general circulation, J. Geophys. Res., 118, 1158911599, doi:http://dx.doi.org/10.1002/2013JD020526

Web End =10.1002/2013JD020526 http://dx.doi.org/10.1002/2013JD020526

Web End = , 2013.

Alexander, M. J. and Teitelbaum, H.: Observation and analysis of a large amplitude mountain wave event over the Antarctic peninsula, J. Geophys. Res., 112, D21103, doi:http://dx.doi.org/10.1029/2006JD008368

Web End =10.1029/2006JD008368 http://dx.doi.org/10.1029/2006JD008368

Web End = , 2007.

Alexander, M. J. and Teitelbaum, H.: Three-dimensional properties of Andes mountain waves observed by satellite: A case study, J.Geophys. Res., 116, D23110, doi:http://dx.doi.org/10.1029/2011JD016151

Web End =10.1029/2011JD016151 http://dx.doi.org/10.1029/2011JD016151

Web End = , 2011.

Alexander, M. J., Gille, J., Cavanaugh, C., Coffey, M., Craig, C., Eden, T., Francis, G., Halvorson, C., Hannigan, J., Khosravi,R., Kinnison, D., Lee, H., Massie, S., Nardi, B., Barnett, J., Hepplewhite, C., Lambert, a., and Dean, V.: Global estimates of gravity wave momentum ux from High Resolution Dynamics Limb Sounder observations, J. Geophys. Res., 113, D15S18, doi:http://dx.doi.org/10.1029/2007JD008807

Web End =10.1029/2007JD008807 http://dx.doi.org/10.1029/2007JD008807

Web End = , 2008.

Alexander, M. J., Eckermann, S. D., Broutman, D., and Ma, J.: Momentum ux estimates for South Georgia Island mountain waves in the stratosphere observed via satellite, Geophy. Res. Lett., 36, L12816, doi:http://dx.doi.org/10.1029/2009GL038587

Web End =10.1029/2009GL038587 http://dx.doi.org/10.1029/2009GL038587

Web End = , 2009.

Alexander, M. J., Geller, M., McLandress, C., Polavarapu, S., Preusse, P., Sassi, F., Sato, K., Eckermann, S., Ern, M., Hertzog,A., Kawatani, Y., Pulido, M., Shaw, T. A., Sigmond, M., Vincent, R., and Watanabe, S.: Recent developments in gravity-wave effects in climate models and the global distribution of gravity-wave momentum ux from observations and models, Q. J. Roy.Meteor. Soc., 136, 11031124, doi:http://dx.doi.org/10.1002/qj.637

Web End =10.1002/qj.637 http://dx.doi.org/10.1002/qj.637

Web End = , 2010.Aumann, H., Chahine, M., Gautier, C., Goldberg, M., Kalnay,E., McMillin, L., Revercomb, H., Rosenkranz, P., Smith, W., Staelin, D., Strow, L., and Susskind, J.: AIRS/AMSU/HSB on the Aqua mission: design, science objectives, data products, and processing systems, IEEE T. Geosci. Remote, 41, 253264, doi:http://dx.doi.org/10.1109/TGRS.2002.808356

Web End =10.1109/TGRS.2002.808356 http://dx.doi.org/10.1109/TGRS.2002.808356

Web End = , 2003.

Barry, R. L., Strother, S. C., and Gore, J. C.: Complex and magnitude-only preprocessing of 2D and 3D BOLD fMRI data at 7T, Magnetic Resonance in Medicine, 67, 867871, doi:http://dx.doi.org/10.1002/mrm.23072

Web End =10.1002/mrm.23072 http://dx.doi.org/10.1002/mrm.23072

Web End = , 2012.

Baumgaertner, A. J. G. and McDonald, A. J.: A gravity wave climatology for Antarctica compiled from Challenging Minisatellite Payload/Global Positioning System (CHAMP/GPS) radio occultations, J. Geophys. Res., 112, D05103, doi:http://dx.doi.org/10.1029/2006JD007504

Web End =10.1029/2006JD007504 http://dx.doi.org/10.1029/2006JD007504

Web End = , 2007.

Brigham, E. O.: The Fast Fourier Transform, Prentice-Hall Inc.,1974.

Brown, R. A., Lauzon, M. L., and Frayne, R.: A General Description of Linear Time-Frequency Transforms and Formulation of a Fast, Invertible Transform That Samples the Continuous S-Transform Spectrum Nonredundantly, IEEE Trans. Sig. Proc., 58, 281290, doi:http://dx.doi.org/10.1109/TSP.2009.2028972

Web End =10.1109/TSP.2009.2028972 http://dx.doi.org/10.1109/TSP.2009.2028972

Web End = , 2010.

Butchart, N., Charlton-Perez, A. J., Cionni, I., et al.: Multimodel climate and variability of the stratosphere, J. Geophys. Res., 116, D05102, doi:http://dx.doi.org/10.1029/2010JD014995

Web End =10.1029/2010JD014995 http://dx.doi.org/10.1029/2010JD014995

Web End = , d05102, 2011. Clerbaux, C., Boynard, A., Clarisse, L., George, M., Hadji-Lazaro,J., Herbin, H., Hurtmans, D., Pommier, M., Razavi, A., Turquety,S., Wespes, C., and Coheur, P.-F.: Monitoring of atmospheric composition using the thermal infrared IASI/MetOp sounder, Atmos. Chem. Phys., 9, 6041-6054, doi:http://dx.doi.org/10.5194/acp-9-6041-2009

Web End =10.5194/acp-9-6041-2009 http://dx.doi.org/10.5194/acp-9-6041-2009

Web End = , 2009.

Eckermann, S. D. and Preusse, P.: Global Measurements of Stratospheric Mountain Waves from Space, Science, 286, 15341537, doi:http://dx.doi.org/10.1126/science.286.5444.1534

Web End =10.1126/science.286.5444.1534 http://dx.doi.org/10.1126/science.286.5444.1534

Web End = , 1999.

Ern, M., Preusse, P., Alexander, M. J., and Warner, C. D.: Absolute values of gravity wave momentum ux derived from satellite data, J. Geophys. Res., 109, D20103, doi:http://dx.doi.org/10.1029/2004JD004752

Web End =10.1029/2004JD004752 http://dx.doi.org/10.1029/2004JD004752

Web End = , 2004.

Ern, M., Preusse, P., Gille, J. C., Hepplewhite, C. L., Mlynczak,M. G., Russell, J. M., and Riese, M.: Implications for atmospheric dynamics derived from global observations of gravity wave momentum ux in stratosphere and mesosphere, J. Geophys. Res., 116, D19107, doi:http://dx.doi.org/10.1029/2011JD015821

Web End =10.1029/2011JD015821 http://dx.doi.org/10.1029/2011JD015821

Web End = , 2011. Fritts, D. C. and Alexander, M. J.: Gravity wave dynamics and effects in the middle atmosphere, Rev. Geophys., 41, 1003, doi:http://dx.doi.org/10.1029/2001RG000106

Web End =10.1029/2001RG000106 http://dx.doi.org/10.1029/2001RG000106

Web End = , 2003.

Fritts, D. C., Riggin, D. M., Balsley, B. B., and Stockwell,R. G.: Recent results with an MF radar at McMurdo, Antarctica: Characteristics and variability of motions near 12-hour period in the mesosphere, Geophy. Res. Lett., 25, 297300, doi:http://dx.doi.org/10.1029/97GL03702

Web End =10.1029/97GL03702 http://dx.doi.org/10.1029/97GL03702

Web End = , 1998.

Geller, M., Alexander, M. J., Love, P., Bacmeister, J., Ern, M., Hertzog, A., Manzini, E., Preusse, P., Sato, K., Scaife, A., and Zhou,T.: A Comparison between Gravity Wave Momentum Fluxes in Observations and Climate Models, J. Climate, 26, 63836405, doi:http://dx.doi.org/10.1175/JCLI-D-12-00545.1

Web End =10.1175/JCLI-D-12-00545.1 http://dx.doi.org/10.1175/JCLI-D-12-00545.1

Web End = , 2013.

Gibson, P. C., Lamoureux, M. P., and Margrave, G. F.: Letter to the editor: Stockwell and wavelet transforms, J. Fourier Anal. Appl., 12, 713721, doi:http://dx.doi.org/10.1007/s00041-006-6087-9

Web End =10.1007/s00041-006-6087-9 http://dx.doi.org/10.1007/s00041-006-6087-9

Web End = , 2006. Goodyear, B. G., Zhu, H. M., Brown, R. A., and Mitchell, J. R.:

Removal of phase artifacts from fMRI data using a Stockwell transform lter improves brain activity detection, Magn. Reson. Med., 51, 1621, doi:http://dx.doi.org/10.1002/mrm.10681

Web End =10.1002/mrm.10681 http://dx.doi.org/10.1002/mrm.10681

Web End = , 2004.

Hertzog, A., Boccara, G., Vincent, R. A., Vial, F., and Cocquerez,P.: Estimation of gravity wave momentum ux and phase speeds from quasi-Lagrangian stratospheric balloon ights Part II: Results from the Vorcore campaign in Antarctica., J. Atmos. Sci., 65, 30563070, doi:http://dx.doi.org/10.1175/2008JAS2710.1

Web End =10.1175/2008JAS2710.1 http://dx.doi.org/10.1175/2008JAS2710.1

Web End = , 2008.

Atmos. Meas. Tech., 9, 25452565, 2016 www.atmos-meas-tech.net/9/2545/2016/

N. P. Hindley et al.: A 2-D Stockwell transform for AIRS gravity wave analysis 2565

Hertzog, A., Alexander, M. J., and Plougonven, R.: On the Intermittency of Gravity Wave Momentum Flux in the Stratosphere, J.Atmos. Sci., 69, 34333448, doi:http://dx.doi.org/10.1175/JAS-D-12-09.1

Web End =10.1175/JAS-D-12-09.1 http://dx.doi.org/10.1175/JAS-D-12-09.1

Web End = , 2012.

Hindley, N. P., Wright, C. J., Smith, N. D., and Mitchell, N. J.: The southern stratospheric gravity wave hot spot: individual waves and their momentum uxes measured by COSMIC GPS-RO, Atmos. Chem. Phys., 15, 7797-7818, doi:http://dx.doi.org/10.5194/acp-15-7797-2015

Web End =10.5194/acp-15-7797- http://dx.doi.org/10.5194/acp-15-7797-2015

Web End =2015 , 2015.

Hoffmann, L. and Alexander, M. J.: Retrieval of stratospheric temperatures from Atmospheric Infrared Sounder radiance measurements for gravity wave studies, J. Geophys. Res., 114, D07105, doi:http://dx.doi.org/10.1029/2008JD011241

Web End =10.1029/2008JD011241 http://dx.doi.org/10.1029/2008JD011241

Web End = , 2009.

Hoffmann, L., Xue, X., and Alexander, M. J.: A global view of stratospheric gravity wave hotspots located with Atmospheric Infrared Sounder observations, J. Geophys. Res., 118, 416434, doi:http://dx.doi.org/10.1029/2012JD018658

Web End =10.1029/2012JD018658 http://dx.doi.org/10.1029/2012JD018658

Web End = , 2013.

Hoffmann, L., Alexander, M. J., Clerbaux, C., Grimsdell, A. W., Meyer, C. I., Rler, T., and Tournier, B.: Intercomparison of stratospheric gravity wave observations with AIRS and IASI, Atmos. Meas. Tech., 7, 45174537, doi:http://dx.doi.org/10.5194/amt-7-4517-2014

Web End =10.5194/amt-7-4517-2014 http://dx.doi.org/10.5194/amt-7-4517-2014

Web End = , 2014.

Jiang, J. H., Wu, D. L., and Eckermann, S. D.: Upper Atmosphere Research Satellite (UARS) MLS observation of mountain waves over the Andes, J. Geophys. Res.-Atmos., 107, SOL 15-1SOL 15-10 , doi:http://dx.doi.org/10.1029/2002JD002091

Web End =10.1029/2002JD002091 http://dx.doi.org/10.1029/2002JD002091

Web End = , 2002.

Kocahan, O., Coskun, E., and Oezder, S.: Prole Measurement of Objects by Using Stockwell and Continuous Wavelet Transforms, 16th Signal Processing and Communications Applications Conference, Aydin, Turkey, 2022 April 2008, 804807, 2008.

Kuyuk, H. S.: On the use of Stockwell transform in structural dynamic analysis, Sadhana-Academy Proceedings in Engineering Sciences, 40, 295306, doi:http://dx.doi.org/10.1007/s12046-014-0301-2

Web End =10.1007/s12046-014-0301-2 http://dx.doi.org/10.1007/s12046-014-0301-2

Web End = , 2015.Liu, Y.: Localization Operators for Two-Dimensional Stockwell

Transforms, in: New Developments in Pseudo-Differential Operators, edited by: Rodino, L. and Wong, M. W., Operator Theory Advances and Applications, 189, 287296, 2009.

Liu, Y. and Wong, M. W.: Inversion formulas for two-dimensional Stockwell transforms, in: Pseudo-Differential Operators: Partial Differential Equations and time-frequency analysis, edited by: Rodino, L., Schulze, B. W., and Wong, M. W., Fields Institute Communications, 52, 323330, 2007.

Mansinha, L., Stockwell, R. G., and Lowe, R. P.: Pattern analysis with two-dimensional spectral localisation: Applications of two-dimensional S transforms, Physica A, 239, 286295, doi:http://dx.doi.org/10.1016/S0378-4371(96)00487-6

Web End =10.1016/S0378-4371(96)00487-6 http://dx.doi.org/10.1016/S0378-4371(96)00487-6

Web End = , 1997a.

Mansinha, L., Stockwell, R. G., Lowe, R. P., Eramian, M., and Schincariol, R. A.: Local S-spectrum analysis of 1-D and 2-D data, Phys. Earth Planet. In., 103, 329336, doi:http://dx.doi.org/10.1016/S0031-9201(97)00047-2

Web End =10.1016/S0031- http://dx.doi.org/10.1016/S0031-9201(97)00047-2

Web End =9201(97)00047-2 , 1997b.

McDonald, A. J.: Gravity wave occurrence statistics derived from paired COSMIC/FORMOSAT3 observations, J. Geophys. Res., 117, D15106, doi:http://dx.doi.org/10.1029/2011JD016715

Web End =10.1029/2011JD016715 http://dx.doi.org/10.1029/2011JD016715

Web End = , 2012.

Pinnegar, C. R. and Mansinha, L.: The S-transform with windows of arbitrary and varying shape, Geophysics, 68, 381385, doi:http://dx.doi.org/10.1190/1.1543223

Web End =10.1190/1.1543223 http://dx.doi.org/10.1190/1.1543223

Web End = , 2003.

Plougonven, R., Arsac, A., Hertzog, A., Guez, L., and Vial, F.: Sensitivity study for mesoscale simulations of gravity waves above

Antarctica during Vorcore, Q. J. Roy. Meteor. Soc., 136, 1371 1377, doi:http://dx.doi.org/10.1002/qj.639

Web End =10.1002/qj.639 http://dx.doi.org/10.1002/qj.639

Web End = , 2010.

Plougonven, R., Hertzog, A., and Guez, L.: Gravity waves over Antarctica and the Southern Ocean: consistent momentum uxes in mesoscale simulations and stratospheric balloon observations,Q. J. Roy. Meteor. Soc., 139, 101118, doi:http://dx.doi.org/10.1002/qj.1965

Web End =10.1002/qj.1965 http://dx.doi.org/10.1002/qj.1965

Web End = , 2013.

Sato, K., Tateno, S., Watanabe, S., and Kawatani, Y.: Gravity Wave Characteristics in the Southern Hemisphere Revealed by a High-Resolution Middle-Atmosphere General Circulation Model., J. Atmos. Sci., 69, 13781396, doi:http://dx.doi.org/10.1175/JAS-D-11-0101.1

Web End =10.1175/JAS-D-11- http://dx.doi.org/10.1175/JAS-D-11-0101.1

Web End =0101.1 , 2012.

Shutts, G. J. and Vosper, S. B.: Stratospheric gravity waves revealed in NWP model forecasts, Q. J. Roy. Meteor. Soc., 137, 303317, doi:http://dx.doi.org/10.1002/qj.763

Web End =10.1002/qj.763 http://dx.doi.org/10.1002/qj.763

Web End = , 2011.

Stockwell, R., Mansinha, L., and Lowe, R. P.: Localization of the complex spectrum: the S-Transform, IEEE Trans. Sig. Proc., 44, 9981001, 1996.

Stockwell, R. G.: S-Transform Analysis of Gravity Wave Activity from a Small Scale Network of Airglow Imagers, PhD thesis, University of Western Ontario, 1999.

Stockwell, R. G.: A basis for efcient representation of the S-transform, Digital Signal Processing, 17, 371393, doi:http://dx.doi.org/10.1016/j.dsp.2006.04.006

Web End =10.1016/j.dsp.2006.04.006 http://dx.doi.org/10.1016/j.dsp.2006.04.006

Web End = , 2007.

Wright, C. J.: Detection Of Stratospheric Gravity Waves Using HIRDLS Data, PhD thesis, Trinity College, University of Oxford, 2010.

Wright, C. J.: A one-year seasonal analysis of martian gravity waves using MCS data, Icarus, 219, 274282, doi:http://dx.doi.org/10.1016/j.icarus.2012.03.004

Web End =10.1016/j.icarus.2012.03.004 http://dx.doi.org/10.1016/j.icarus.2012.03.004

Web End = , 2012.

Wright, C. J. and Gille, J. C.: Detecting overlapping gravity waves using the S-Transform, Geophy. Res. Lett., 40, 18501855, doi:http://dx.doi.org/10.1002/grl.50378

Web End =10.1002/grl.50378 http://dx.doi.org/10.1002/grl.50378

Web End = , 2013.

Wright, C. J., Osprey, S. M., and Gille, J. C.: Global distributions of overlapping gravity waves in HIRDLS data, Atmos. Chem. Phys., 15, 8459-8477, doi:http://dx.doi.org/10.5194/acp-15-8459-2015

Web End =10.5194/acp-15-8459-2015 http://dx.doi.org/10.5194/acp-15-8459-2015

Web End = , 2015. Wright, C. J., Hindley, N. P., Moss, A. C., and Mitchell, N. J.: Multi-instrument gravity-wave measurements over Tierra del Fuego and the Drake Passage Part 1: Potential energies and vertical wavelengths from AIRS, COSMIC, HIRDLS, MLS-Aura, SAAMER, SABER and radiosondes, Atmos. Meas. Tech., 9, 877-908, doi:http://dx.doi.org/10.5194/amt-9-877-2016

Web End =10.5194/amt-9-877-2016 http://dx.doi.org/10.5194/amt-9-877-2016

Web End = , 2016a.

Wright, C. J., Hindley, N. P., and Mitchell, N. J.: Combining AIRS and MLS observations for three-dimensional gravity wave measurement, Geophy. Res. Lett., 43, 884893, doi:http://dx.doi.org/10.1002/2015GL067233

Web End =10.1002/2015GL067233 http://dx.doi.org/10.1002/2015GL067233

Web End = , 2015GL067233, 2016b.

Wu, D. L.: Mesoscale gravity wave variances from AMSUA radiances, Geophy. Res. Lett., 31, 19448007, doi:http://dx.doi.org/10.1029/2004GL019562

Web End =10.1029/2004GL019562 http://dx.doi.org/10.1029/2004GL019562

Web End = , l12114, 2004.

Yan, A., Zhou, W., Yuan, Q., Yuan, S., Wu, Q., Zhao, X., and Wang,J.: Automatic seizure detection using Stockwell transform and boosting algorithm for long-term EEG, Epilepsy Behav., 45, 8 14, doi:http://dx.doi.org/10.1016/j.yebeh.2015.02.012

Web End =10.1016/j.yebeh.2015.02.012 http://dx.doi.org/10.1016/j.yebeh.2015.02.012

Web End = , 2015.

www.atmos-meas-tech.net/9/2545/2016/ Atmos. Meas. Tech., 9, 25452565, 2016