A Brief Overview of Gravity Wave Retrieval Techniques From Observations

Atmospheric gravity waves (GWs) are important in driving the middle and upper atmosphere dynamics on Earth. Here, we provide a brief review of the most common techniques of retrieving gravity wave activity from observations. Retrieval of gravity wave activity from observations is a multi-step process. First, the background fields have to be removed as the retrieval of the wave activity highly depends on this. Second, since a broad spectrum of internal waves contributes to atmospheric fluctuations, the contribution of GWs has to be extracted carefully. We briefly discuss the strengths, limitations/barriers, and applications of each technique. We also outline some future research questions to improve the treatment of these wave extraction methods.


INTRODUCTION
Atmospheric gravity waves (GWs) are present in all stably stratified atmospheres and play a pivotal role in their mixing, energy budget, and momentum budget (Medvedev and Yiğit, 2019). The dynamical and thermodynamical role of GWs in the vertical coupling between the lower and upper atmosphere is being increasingly studied in a global sense in a number of studies (Yiğit et al., 2016, and see the references therein). Numerical models of varying flavors are used to study GW processes in the atmosphere (e.g., Yiğit et al., 2009;Hickey et al., 2011;Heale et al., 2014;Gavrilov and Kshevetskii, 2015;Yiğit et al., 2021). Coarse-grid general circulation models require incorporation of GW models in order to account for subgrid-scale GW and to reproduce the variable and mean structure of the middle and upper atmosphere. GW observations are required to validate these models.
Space-based and ground-based observations are both used to determine GW climatology (Vincent and Fritts, 1987). Orbiting satellites can provide global coverage, while ground-based measurements (e.g., lidars, radars (MST (Mesosphere-Stratosphere-Troposphere), MU (Middle and Upper atmosphere), Meteor radar, etc.), rocket sounding, or radiosondes) provide local observations but with different height and temporal resolution. Weather balloon-driven radiosondes were extensively used to determine GW parameters in the lower atmosphere (~40 km) due to excellent height resolution (Wang et al., 2005). Middle and upper atmospheres are mostly observed by satellites (Fritts and Alexander, 2003).
Here, we provide a brief review of techniques (Section 2) of GW retrieval from different observations. We also show the applicability of these techniques to the Martian atmosphere. The final section provides a summary with concluding remarks.

GW ANALYSIS TECHNIQUES
Most observational retrieval of GW activity involves separation of the instantaneous observation of an observable into a background mean, which is assumed to be a slowly varying part and fluctuation (wave) part.
Considering temperature as a representative field variable, we can write where T is the background (mean) temperature, and T′ denotes the wave-induced term, in principle including contributions from a broad spectrum of waves, for example, tides, planetary waves, and GWs. One challenge in analyzing atmospheric gravity waves is the effective separation of these large scale waves from gravity waves. Sometimes, depending on the nature of the observations, the fluctuating part is directly attributed to GWs.

Polynomial Fit
Separation into a background and a wave part can be done either vertically or horizontally. In the vertical direction, it is common to use a polynomial fit by using the least square method. Subtraction of this fitted polynomial from observations provides a proxy for gravity wave activity.
Polynomial fitting is a form of regression analysis in which the relation between a dependent variable and an independent variable is modeled into an n th degree polynomial of the dependent variable using the least square method. The general equation of the polynomials can be written as: f x ( ) n n 0 a n x n + ϵ a 0 + a 1 x 1 + a 2 x 2 + a 3 x 3 + / + a n x n + ϵ, where n = 0, 1, 2, 3,. . ., that is, any real valued integer, ϵ is the random error conditioned on the independent variable x of the highest degree of polynomial n, and a's are the set of polynomial coefficients. Different orders of polynomials are decided based on the height range of the analysis. For example, if L is the total height range, then by fitting second order polynomials, one can retrieve waves of maximum wavelength L/2. In other words, perturbations with vertical wavelengths greater than 2L may be removed by the use of a second order polynomial. Likewise, a third order and fourth order polynomial fit may remove perturbations with vertical wavelengths greater than L and 2L/3, respectively. Any polynomials of order greater than four may remove the signals of inertia-gravity waves from the analysis (Guest et al., 2000). To calculate GW source spectrum to use in numerical models by using radiosonde observation in the lower atmosphere of the Southern Hemisphere, Pfenninger et al. (1999) used different order polynomials to be fitted as the background for different atmospheric layers, such as lower order polynomials (order of 2 or 3) as background for the troposphere and comparatively higher order polynomial fit (above 3) in the stratosphere. A number of studies have used a second order polynomial fit to observations in order to determine the mean field (e.g., Allen and Vincent, 1995;Vincent and Joan Alexander, 2000;Wang and Geller, 2003;Yi, 2005, 2007;Zhang et al., 2010).
It is also possible to use higher order polynomials in GW analysis. A third order polynomial fit has been used in different geolocational radiosonde observations to define atmospheric background. For example, Koushik et al. (2019) reported inertia gravity waves in the lower stratosphere using radiosonde from the measurement of vertical profiles of temperature, ozone, zonal wind, and meridional wind. A wavelet analysis (Section 2.6) was applied to the retrieved perturbed quantities in order to characterize GWs.
To investigate the applicability of limited angle tomography in GW research by using GLORIA limb sounder, Krisch et al. (2017) used a one-dimensional Savitzky-Golay filter (Savitzky and Golay, 1964) to determine the atmospheric background. This technique applies a cubic polynomial fit in a small subset of data and continues this process to the next subsets to smooth the data and eventually create the background planetary scale wave shapes. They applied this filter onto 25 vertical and 60 horizontal neighboring points in all three spatial directions. In this technique, they have been able to separate the waves having horizontal wavelengths of~750 km and vertical wavelengths of~3 km. The use of this filter has also been seen in the work by He et al. (2020). They stated the complete gravity wave spectrum characteristics that were obtained from three consecutive stages of balloon movement for the first time.
A fourth order polynomial fit was used to examine the properties of inertia gravity waves in the lower stratosphere using high resolution Omega-sonde ozone soundings (Guest et al., 2000). They fitted the polynomial to zonal wind and temperature measurements over the intervals of 0-9 km and 12-30 km. The difference between soundings and the background profiles defines the wave field which can be filtered out by using a Fourier transform to remove the very high or very low frequency oscillations. Dutta et al. (2017) reported the use of different order polynomial fits (2-9) as the background on observed instantaneous wind profiles by global positioning system radiosondes over Hyderabad (17.4°N, 78.5°E). They observed the difference of wave parameters calculated by varying order of polynomials, and found reasonably reduced difference when a "Butterworth" filter was applied to extract inertia gravity wave components.

Use of Polynomial Fit in Martian GW Analysis
In a similar fashion, polynomial fits to determine the background state have been used in the studies of GWs in other planetary atmospheres such as Mars. Nakagawa et al. (2020) reported vertical wave propagation in the middle atmosphere of Mars by using MAVEN/IUVS observations. They presented a study on wave perturbation in a temperature profile from 20 to 140 km. They used second order polynomial fit as background temperature. A least square third order fit has been used as background temperature by Creasey et al. (2006) to analyze gravity wave activity in Mars Global Surveyor (MGS) radio occultation experiments. Ando et al. (2012) calculated vertical wavenumber spectra from temperature profiles acquired by the same MGS radio occultation experiments in an altitude range between 3 and 32 km. Authors in the latter two studies fitted a third order polynomial fit as the shape of the background temperature. An increased order (seventh order) polynomial fit has been used in the work by Yiğit et al. (2015) in the analysis of the Martian thermospheric gravity wave using the NGIMS instrument on board MAVEN. They demonstrated that the seventh order polynomial fit is an adequate estimate of the background atmospheric density profile fitted on a logarithmic scale of CO 2 number density. On Mars, it is demonstrated that almost identical features of GWs have been found from Ar and CO 2 density perturbations if a seventh order polynomial fit is used (England et al., 2017).

Spectral Filtering or High/Low/Band Pass Filtering
Another common method to separate gravity wave contributions from the large scale waves is to apply a high/low-pass filter of any cutoff wavelength depending on the properties of the filter. Generally for GW analysis, high-pass filters are used in the time domain to remove low frequency components of other large scale waves, and in the height domain to remove contribution from tides (Hirota, 1984;Hirota and Niki, 1985;Eckermann et al., 1995). Low-pass filter is used in the height domain to remove lower wavelengths. Tsuda et al. (1994) delineated vertical and horizontal propagation characteristics by analyzing radiosonde measurement over tropics using a high-pass filter to determine the fluctuation. Hodograph analysis works best for a monochromatic wave or dominant gravity wave energy packet; therefore, Tsuda et al. (1994) used high-pass filtering before applying hodograph analysis to determine the propagation characteristics of the wave. Using these techniques, they found that most GWs were generated in the middle of the troposphere and can propagate upward into the stratosphere. The typical observed vertical wavelength range was 2-2.5 km, while horizontal phase velocities were 5-7 ms −1 at the troposphere and 12 ms −1 at the stratosphere. Eckermann et al. (1995) used this technique in isolating fluctuations from vertical profiles of temperature and wind in the altitude range of 20-60 km at 15 sites measured by meteorological rocket sounding. Later, Tsuda et al. (2000) applied high-pass filter with a cutoff of 10 km wavelength to remove contributions from tides and planetary waves. He also used a low-pass filter with a cutoff wavelength of~200 m in the troposphere and2 -3 km in the stratosphere to extract GW temperature perturbations from GPS/MET temperature profiles aiming to calculate potential energy and Brunt-Väisäla frequency, which are an indicator of atmospheric stability. However, appearance of equatorial Kelvin waves (vertical wavelengths of~3-5) may also be found in the calculation of GW perturbations by using this method (Holton et al., 2002). Also, this technique removes GWs having vertical wavelengths lower than the cutoff, which may be important in the tropical GW observations. De la Torre et al. (2006) identified horizontal wavelengths of two different modes of waves over Andes (mountain wave~40 km and IGW~2 00 km) by using the band-pass filtered temperature perturbation retrieved from CHAMP (Challenging Minisatellite Payload), SAC-C (Satélite dé Aplicaciónes Cientificás-Ć), and Low Earth Orbit (LEO) satellites, between May 2001 and February 2006. A similar filter called empirical mode decomposition (EMD), which is a part of Hilbert Huang Transformations (Huang et al., 1998), has been used in the work by Nayak and Yiğit (2019) in analyzing SABER vertical profiles of temperature to calculate the relative temperature fluctuations. GWs calculated by this method effectively rule out the contribution of PWs or tides and have predominant vertical wavelengths of~5-15 km.

Least Square Fitting
This method is similar to a polynomial fit. It is often used in space-based observation and can be used as a horizontal filter. In this technique, atmospheric background is estimated by a simple least square fit imposed on to instantaneous profiles of global measurement to remove wave contributions of larger scale waves having zonal wavenumbers typically less than 7. Binning subtracted data in a specific latitude-longitude, longitude-altitude, or latitude-altitude grid yields the global map of GWs of particular interest Kumar, 2012, 2013).
John and Kumar (2013) discussed two specific methods in the context of global maps of gravity wave potential energy from temperature profiles observed by a single Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument on board the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite. Here, the first method is described as the LSF method and the second one is the spectral filtering (Section 2.2). In this LSF method, the observation of the instantaneous profiles first need to be binned in specific grid points. John and Kumar (2013) binned the instantaneous temperature profile over the globe in a 5°× 20°( latitude and longitude) grid first. Then, they estimated 0-6 zonal wavenumber components from the gridded data by employing a simple least square fit aiming to construct planetary scale wave amplitudes and phases. This constructed background is then subtracted from the instantaneous profiles falling into the respective grids to obtain perturbations due to gravity waves only. This technique is applied in analyzing gravity wave variation in the stratosphere and mesosphere (Yamashita et al., 2013;Liu et al., 2014). A spectral analysis or wavelet analysis is required afterward to remove noise and isolate wave-like features on each vertical temperature profile to determine dominant vertical wavelengths (section 2.6).

Kalman Filtering
Kalman filter is an algorithm which carries an advantage of getting prediction by using observation. The use of Kalman filter has been seen in multidisciplinary research fields including gravity wave observations (Fetzer and Gille, 1994;Preusse et al., 2002;Ern et al., 2004Ern et al., , 2006Krebsbach and Preusse, 2007). Fetzer and Gille (1994) provided a description of the variability of the variance observed in temperature profiles of LIMS (Limb Infrared Monitor of the Stratosphere). Gravity wave Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 824875 profiles were identified as the difference between instantaneous and background profiles. The background was constructed from daily mapped zonal planetary wave 0-6 temperature coefficients estimated by a Kalman filter which, in general, supplies information about the larger scales.
In a detailed analysis of stratospheric gravity waves originated by mountain waves, Preusse et al. (2002) described an estimation of sensitivity to GWs from emission of limb-sounding CRISTA-1 (Cryogenic Infrared Spectrometers and Telescopes for the Atmosphere satellite). To account for planetary scale wave structures including zonal mean (zonal wavenumber 0) to higher scale planetary waves having zonal wavenumbers up to 1-6, a "Kalman" filter of 0-6 wavenumber is used to obtain an estimation of background atmosphere [following Rodgers (1977); Ern (1993)], as planetary wave related atmospheric zonal modes in the stratosphere are expected to have zonal wavenumbers up to 4. Ern et al. (2004) calculated the total vertical flux of horizontal momentum due to GWs at 25 km altitude following the same technique. Similarly, a comparison study of gravity wave momentum flux (GWMF) between observation and model analysis (CRISTA-1 and CRISTA-2 with GWMF calculated using Warner and McIntyre parameterization scheme) was presented by Ern et al. (2006), where the same 0-6 zonal wavenumber "Kalman" filter had been applied to construct the atmospheric background from instantaneous observations. A detailed description of validation of "Kalman" filtering can be found in appendix A of that study. Similar use of this filter had been done by Krebsbach and Preusse (2007) in a spectral analysis of intra-annual, semiannual, and inter-annual GW activity from SABER temperature observation.

Longitude-Time 2D Fourier Decomposition
Two-dimensional Fourier analysis has successfully been applied to the investigations of global-scale waves in a number of studies (e.g., Lait and Stanford, 1988a,b;Garcia et al., 2005;Ern et al., 2008Ern et al., , 2009aErn et al., ,b, 2015Ern et al., , 2016Meyer et al., 2018). Ern et al. (2011) provided a detailed description of this technique in order to retrieve global properties of GW absolute momentum flux (GWMF). A two-dimensional Fourier decomposition in terms of longitude and time had been performed to obtain shorter period planetary scale waves with periods as short as 1-2 days. In this investigation, authors used a high-resolution data set of HIRDLS (High Resolution Dynamics Limb Sounder)/Aura temperatures that cover the upper troposphere and the stratosphere along with SABER temperature observations that cover both the stratosphere and mesosphere. A comparison with CRISTA has also been reported in terms of GWMF with HIRDLS and SABER. An example of the use of this method has been shown in Figure 1. Zonal averages of the GW squared temperature amplitudes, vertical wavelengths, reciprocal horizontal wavelengths, and absolute momentum fluxes have been retrieved from the SABER data by using this technique. Lower atmospheric (below 40 km) GW activity had been detected during the month of July 2006 ( Figure 1I).

Wavelet Analysis
A well-known method to retrieve gravity wave parameters from observations after background removal is the wavelet analysis (Zink and Vincent, 2001;Moldovan et al., 2002;De la Torre et al., 2006;Kramer et al., 2016;Koushik et al., 2019;Colligan et al., 2020). The procedure of the wavelet analysis technique is fully described in the work by Torrence and Compo (1998) and in the book of Nappo (2013). A detailed description of its application in the context of gravity waves is described by Zink and Vincent (2001). Regions in high wind variance in height-wavenumber spectrum are identified as gravity wave packet in their analyses.
The first step of doing a wavelet analysis is the choice of the mother wavelet (MW). Then, the Fourier transform of the MW and Fourier transform of time series data are needed to calculate a minimum scale and to choose all other scales. For each scale, computing daughter wavelet at that scale and normalizing daughter wavelet by the square root of the total wavelet variance are required. Then, multiplying by the Fourier transform of the time series and converting the transform back to the real space can visualize the dominant waves present in the observations. Generally, for analysis, one may use a Morlet wavelet because horizontal perturbations of any wind profiles caused by gravity wave packets are primarily amplitude modulated sine functions (Zink and Vincent, 2001). The wavelet transform of an arbitrary signal at space f(x) can be defined as (Moldovan et al., 2002): where ψ(x) is the basis wavelet, p is called the dilation or compression scale factor, and q represents the translation at space. Morlet wavelet can be defined as: where ω 0 denotes angular frequency of the wavelet. When q is known in this equation, p can be inferred from the local wavelength and the wavelet transform would have the maximum amplitude (Torresani, 1995).

Stockwell or S-Transform
This is a spectral analysis technique, similar to the continuous Gaussian wavelet analysis. Alexander et al. (2008) and Wright et al. (2010) applied the S-transform to temperature perturbations in analyzing the background of the observations by HIRDLS/ Aura. Sinusoidal functions modulated with a Gaussian width inversely proportional to the wavelengths form the basis of an S-transform. Fourier transform may be obtained by the spatial integral of the S-transform. Gaussian is translated in a way that it can provide localized information of spectra (Alexander et al., 2008). Wang et al. (2006) provided detailed mathematical analysis and used a combination of S-transform and hodograph analysis for estimating gravity wave horizontal propagation directions.

Maximum Entropy Method/Harmonic Analysis
The MEM/HA method is also applicable only after the background removal and has extensively been used in a number of studies by Eckermann and Preusse (1999); Preusse et al. (2000Preusse et al. ( , 2001a; and Ern et al. (2004Ern et al. ( , 2011Ern et al. ( , 2015. MEM/HA is a combination of maximum entropy method [more details of MEM can be found in the work by Press and Teukolsky (1992)] and harmonic analysis, which contain certain advantages compared with the Fourier decomposition and the standard deviation, and is well suited for studying a single monochromatic wave of arbitrary wavelength from retrieved temperature (Ern et al., 2016). This technique is useful to obtain the height profiles of the amplitudes, phases, and vertical wavelengths of the two largest oscillations present in any given pair profiles of temperature (Preusse et al., 2001a). A detailed study of the method can be found in the work by Preusse et al. (2002). By following this technique and using CRISTA retrieved Kalman detrended background temperature, vertical wavelengths of 3-5 km had been detected and large GW amplitudes were found over the southernmost part of South America. This method obeys the dispersion relation for the linear 2D mountain waves. We show an example of MEM/HA analysis in Figure 2 adapted from this study. A dominant peak of vertical wavelength was found at 5 km (right panel-bottom to top vertical lines) and at 10 km.

Stokes Parameter Analysis
Historically, Vincent and Fritts (1987) utilized Stokes parameters for the first time in the context of atmospheric gravity waves. By utilizing vector wind radar observations over the region of Adelaide during the period of November 1983 to December 1984, they calculated GW motions and all four Stokes parameters directly or indirectly. Stokes parameter analysis for GWs are generally done in wavenumber domain, not in height domain, because all the parameters can be extended in the Fourier spectral domain (Eckermann and Vincent, 1989). Horizontal velocity fluctuations should be polarized for gravity waves. When the analysis is done in Fourier spectral domain, one may decompose the waves at some extent. Additionally, cross spectral relationship between the zonal component and potential temperature may be established by Stokes parameter analysis (Cho et al., 1999). Stokes parameters were developed originally to explain the state of polarization of plain electromagnetic waves (EM waves) by G. G. Stokes in 1852 (Jackson, 1999). The gravity wave field may be partially polarized which triggers information retrieval by Stokes parameter analysis analogously as GWs are considered plain transverse waves, moving horizontally or vertically in nature (Colligan et al., 2020). The use of Stokes parameters can provide the polarization relations of GWs similar to EM waves and calculate the horizontal propagation direction of waves which is related to wind perturbation hodograph axial ratio (Vincent and Fritts, 1987;Eckermann and Hocking, 1989;Zink and Vincent, 2001;Wang et al., 2005;Colligan et al., 2020). There might be one way to check the presence of coherent waves in the calculation by the measurement of degree of polarization which can be calculated by using Stokes parameter analysis (Vincent and Joan Alexander, 2000).

Other Methods
There are other methods including the harmonic analysis, bispectral analysis, rotary spectral analysis, hodograph analysis, and radars. The harmonic analysis is a branch of mathematics where the representation of superposition of basic waves is discussed by the employment of Fourier series and Fourier transforms. As in the real atmosphere, GWs exist as a superposition with other large scale waves; the harmonic analysis technique is useful in tidal analysis or gravity wave analysis to characterize wave frequency. The bispectral analysis is helpful to explain nonlinear resonant wave-wave interactions (Wüst and Bittner, 2006). The rotary spectral technique (Vincent, 1984) or Hodograph analysis (Hirota and Niki, 1985) is popular to characterize the inertia gravity wave intrinsic frequency and propagation direction, though it demands the presence of a single coherent wave source and does not yield good results at the mixture of different frequencies present in the perturbation profiles (Dutta et al., 2017). In general, in a plot of fluctuation components of winds, an ellipse would be fit and the rotational direction of the ellipse would decide the propagation direction of waves. Clockwise ellipse with height in the Northern Hemisphere and anticlockwise in the Southern Hemisphere indicates upward propagation. Nevertheless, in the presence of multiple wave events, the gravity wave parameters calculated by using this method might be inaccurate as similar features could arise by superposition of linearly polarized waves in a broad spectrum (Eckermann and Hocking, 1989).
Radars have been used extensively in the stratosphere, mesosphere, and the lower thermosphere (Fritts and Alexander, 2003). Combined observation of a lidar and an MST radar at Aberystwyth (52.4°N, 4.1°W) to investigate gravity waves between 2-and 50-km revealed a dominance of long period quasi- Frontiers in Astronomy and Space Sciences | www.frontiersin.org February 2022 | Volume 9 | Article 824875 monochromatic motions in the stratospheric gravity wave field, although this picture is not the same for the troposphere (Mitchell et al., 1994). Fritts et al. (1992) used the Jicamarca MST radar near Lima, Peru (12°S, 77°W), to understand the dynamics of equatorial mesosphere during June and August 1987, and found a very active wave field with considerable gravity wave and tidal motions with periods of~5 min to many hours (less than 1 day). Wang and Fritts (1990) used an MST radar at Poker Flat, Alaska (65°N, 147°W), to measure mesospheric GW momentum fluxes by analyzing the velocity data during the winter, summer, and fall of 1986. They showed that the wave field in the high-latitude mesosphere is highly variable. A radar at Arecibo, Puerto Rico (18°N, 67°W), and an MU radar at Shigaraki, Japan (35°N, 136°E), were used by Tsuda et al. (1989) to measure the saturated gravity wave spectrum. Meteor radar is a relatively inexpensive tool to conduct research in the upper mesosphere and lower thermosphere neutral wind fields by utilizing the reflection of electromagnetic waves from meteor trails in the range of 80-100 km altitude (Avery, 1990

CONCLUSION
The objective of observations is to provide an accurate picture of the real world, which is needed to validate and improve global models of the atmosphere and gravity waves. On the other hand, improved models can help with the interpretation of observations. Also, data assimilation heavily relies on observations of atmospheric field variables. A number of techniques are used to extract GWs from observations. Briefly, two main challenges of GW retrievals are noteworthy. First, the mean field has to be adequately determined so that the fluctuating component of the atmosphere is evaluated. Second, since a broad spectrum of internal waves contributes to atmospheric fluctuations, the contribution of GWs has to be extracted carefully. It is well known from theory and observations that not only gravity waves but also acoustic (infrasonic) waves exist in the atmosphere. These listed waves differ significantly in the limiting cases of low-frequency and high-frequency waves. However, in the intermediate range, it is difficult to distinguish between them. Overall, there are different ways of retrieving GW activity from observations, depending on the nature of observations. More systematic analyses of the intercomparison of the observational retrieval techniques are needed in order to adequately evaluate the weaknesses and strengths of the different approaches. This mini review has presented some of the existing techniques of extracting GWs from atmospheric observations. We have not attempted to provide a detailed mathematical analysis; however, we have reviewed various examples of the applications of the different techniques in order to provide an exposure to the subject of GW analysis.

AUTHOR CONTRIBUTIONS
MS and EY planned and outlined the manuscript; MS researched and wrote the manuscript. EY checked and edited the entire paper content and contributed significantly to all the sections of the paper.

FUNDING
MS is supported by NSF award 1849014, EY is partially funded by NASA grant 80NSSC20K0941.