the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Assessing atmospheric gravity wave spectra in the presence of observational gaps
Irina Strelnikova
Gerd Baumgarten
We present a thorough investigation into the accuracy and reliability of gravity wave (GW) spectral estimation methods when dealing with observational gaps. GWs have a significant impact on atmospheric dynamics, exerting influence over weather and climate patterns. However, empirical atmospheric measurements often suffer from data gaps caused by various factors, leading to biased estimations of the spectral powerlaw exponent (slope) β. This exponent describes how the energy of GWs changes with frequency over a defined range of GW scales. In this study, we meticulously evaluate three commonly employed estimation methods: the fast Fourier transform (FFT), generalized Lomb–Scargle periodogram (GLS), and Haar structure function (HSF). We assess their performance using time series of synthetic observational data with varying levels of complexity, ranging from a signal with one frequency to a number of superposed sinusoids with randomly distributed wave parameters. By providing a comprehensive analysis of the advantages and limitations of these methods, our aim is to provide a valuable roadmap for selecting the most suitable approach for accurate estimations of β from sparse observational datasets.
 Article
(3082 KB)  Fulltext XML
 BibTeX
 EndNote
Gravity waves (GWs) are ubiquitous phenomena that play a crucial role in the dynamics of the Earth's atmosphere, where they impact weather and climate patterns (Hines, 1960; Ern et al., 2018). Various sources, including convection, topography, and jet streams, generate these waves (Crowley and Williams, 1987; Fritts, 1989). As they propagate through the atmosphere, they can break and mix with the surrounding atmosphere, redistributing their energy and momentum. This leads to significant changes in the atmospheric thermodynamics and largescale circulation patterns of the atmosphere, including wind speeds and temperature gradients (Lindzen, 1981; Holton, 1983; Fritts and Alexander, 2003). Observations of these meteorological variables reveal that GWs exist for the most part in the form of a spectrum of superposed waves within a wave packet and occasionally as quasimonochromatic waves (Maekawa et al., 1984; Eckermann and Hocking, 1989). To understand the physical processes that govern the generation, propagation, and dissipation of these wave packets, it is often useful to examine their spectral properties such as their frequencies, amplitudes, and scales (Axford, 1971; Fritts and VanZandt, 1993).
On that note, VanZandt first introduced the concept of a “universal atmospheric GW spectrum” (VanZandt, 1982). This spectrum facilitated efficient parameterizations of how GWs affect the mean atmospheric state (Narendra Babu et al., 2008). For instance, the spectra of GWs are often used in model parameterization, including source spectra parameterization, Lagrangian spectral parameterization, and subgridscale parameterization, enabling the simulation of the dynamics of the middle and upper atmospheres (Beres et al., 2005; Song and Chun, 2008; Houchi et al., 2010). Overall, accurate predictions of GW activity can improve weather forecasting, while they contribute significantly to climate modelling in parameterizing physical processes like turbulence and mixing (Alexander et al., 2002; Smith, 2012; Liu et al., 2014).
This GW spectrum exhibits a power law scaled by an exponent (or slope) β, which describes the rate at which wave energy changes with its wavenumber (or frequency). The basis for this spectrum for atmospheric GWs is supported not only by a strong foundation in theoretical works (Dewan and Good, 1986; Weinstock, 1990; Hines, 1991; Dewan, 1994; Gardner, 1994), but also in observational studies (Smith et al., 1987; Fritts et al., 1988; Gardner et al., 1995; Nastrom et al., 1997; Zhang et al., 2006, 2017). These values of β depend on not only the type of spectra (e.g. temporal or horizontal, or vertical wavenumber) but also the geophysical variables measured (e.g. temperature, horizontal or vertical wind); see Table A1 for a summary. Thus, an accurate estimation of β is essential to validate different theoretical predictions of GW power spectral densities (PSDs) (Dewan and Grossbard, 2000) and improve climate models and weather forecasts (Lindgren et al., 2020).
Determining β from empirical atmospheric measurements is challenging due to various factors, such as the inevitable presence of data gaps, observational noise, and the finiteness of data length and resolution. Data gaps can occur for numerous reasons, including instrumental errors, data transmission issues (e.g. due to weather conditions like clouds in the case of lidar), and signal interference (in the case of radar). When gaps exist in multiscale time series, data points representing certain frequencies are lost, which distorts the spectra and introduces significant bias into the estimation of β (Brown and ChristensenDalsgaard, 1990; Rigling, 2012). To minimize the effect of these gaps on the spectra, datafilling schemes are often applied. Though linear interpolation is usually used to fill in these gaps (Meisel, 1978; Lepot et al., 2017), even adaptively implemented interpolators produce artefacts in the time series at low gap percentages (GPs), which contribute additional bias in the spectra (Schulz and Stattegger, 1997; Hall and Aso, 1999). Bias in spectral estimates can also be caused by other relevant sources, such as spectral leakage, steep spectra (β>2), and insignal components with larger periods than the observed time span T (Klis, 1994).
In this paper, we systematically quantify the advantages and limitations of estimation methods of GW spectra in handling these error sources. We also propose a procedure for selecting unambiguously suitable methods for β estimation. Even though this study is motivated by the analysis of atmospheric GWs, the conclusions and the methods can be generalized to different fields with similar timeseries characteristics in other branches of geophysics. Two commonly used methods are considered, namely the fast Fourier transform (FFT) (Cooley and Tukey, 1965) and the generalized Lomb–Scargle periodogram (GLS) (Zechmeister and Kürster, 2009), as well as the fairly recent Haar structure function (HSF) (Lovejoy and Schertzer, 2012). FFT is the standard method to analyse spectra of evenly sampled data. The Lomb–Scargle periodogram (LS) has been used in many studies as the main analysis method (or as a reference) of GW spectra (e.g. Hall and Aso, 1999; Zhang et al., 2006; Guharay and Sekar, 2011; Qing et al., 2014). As far as the authors are aware, HSF has never been used in atmospheric GW studies. Both GLS and HSF are specifically known to handle unevenly sampled data. In an effort to closely mimic real observations of GWs, we simulate timeseries data with varied levels of complexity, beginning with a signal with one frequency and increasing in complexity to a superposition of sinusoids with randomly distributed frequencies (and phases).
Previous studies have investigated these spectral methods and others for estimating powerlaw spectra and compared their performance using synthetic and observed data. For instance, Zhan et al. (1996) found that using FFT of linearly interpolated signals is the best approach to analyse radar wind data at 50 % GP (only for the case of $\mathit{\beta}=\mathrm{5}/\mathrm{3}$) compared to the correlogram and LS. However, a quantitative analysis of the effect of changing β or the GP was not conducted. Similarly, Munteanu et al. (2016) showed that FFT outperforms LS, Z transform, and discrete Fourier transform in estimating β from Venus' magnetic field data, although the effect of changing β was not considered either, since the powerlaw spectra were not simulated. In contrast, Hébert et al. found that HSF consistently surpassed other methods in estimating β, without the need to interpolate the gapped (simulated palaeoclimate) data for $\mathit{\beta}\in (\mathrm{0},\mathrm{3})$, except in the case of $\mathit{\beta}\in (\mathrm{1},\mathrm{0})$, where they concluded that LS would be the best option (Hébert et al., 2021). Nonetheless, the impact of altering the GP was not quantitatively presented; instead, the skewness of the gaps' (gamma) distribution was used as a parameter to refer to the irregularity of the time series. These three studies showed that the LS method suffers from significant leakage in the case of powerlaw spectra, which persists even when multitaper methods (MTMs) are used.
The rest of the paper is organized as follows: in Sect. 2, we describe the methods used in our study, including a description of FFT, GLS, and HSF. In Sect. 3 we introduce the data simulation procedures. In Sect. 4 we discuss data processing. In Sect. 5 we present the results of our simulations, comparing the performance of these methods in different scenarios. In Sect. 6 we discuss the implications of our findings and provide recommendations for spectral analysis of GW time series with data gaps. Finally, in Sect. 7, we present a summary of our relevant results and conclusions.
2.1 The fast Fourier transform
FFT is the most commonly used method for estimating frequency spectra of evenly sampled data (Cooley and Tukey, 1965). It enables the approximation of a time series sampled from a continuous distribution over discrete time steps, through a series of complex sine and cosine waves with varying frequencies. Under the assumption of a unit sample interval, (forward) FFT converts a time series z_{n} of length N from its original domain (time or space) into a set of coefficients Z_{k} in the (temporal or spatial) frequency domain by employing the following relation:
In our work, FFT will serve as the benchmark spectral estimation method. The expected Fourier transform of a discretized signal is given by the convolution of the true transform and the transform of a Dirac comb window function designating those measurement times (Vanderplas, 2018). In the case of gapped data, the symmetry in the Dirac comb is destroyed, causing the resulting transform to be noisy with incorrect peak positions and heights. Consequently, the true transform of gapped data will not be recoverable. This disadvantage can be bypassed by applying data reconstruction methods such as interpolation, sparse approximation, etc., to approximate the true Fourier transform (Babu and Stoica, 2010). Unfortunately, these reconstruction methods can introduce artefacts into the signal, which depend on the distribution of the gaps and their sizes (Munteanu et al., 2016).
2.2 The generalized Lomb–Scargle periodogram
The GLS periodogram developed by Zechmeister and Kürster (2009) offers a method for estimation of the PSD of unevenly sampled time series. It is a generalization of Lomb's leastsquares approach (Lomb, 1976), which is equivalent to the modified Schuster's periodogram (Schuster, 1898; Scargle, 1982) (based on FFT) in the case of evenly sampled data. GLS produces a spectrum by leastsquares fitting a model of a weighted sinusoid given by
to the time series at each sampled frequency ω. The offset c compensates for the assumption that the mean of the time series $\stackrel{\mathrm{\u203e}}{z}$ is equal to the mean of the fit $\stackrel{\mathrm{\u203e}}{y}$. This floatingmean approach is advantageous, considering that the mean of a periodic signal may change statistically, especially for small N (FerrazMello, 1981). Furthermore, the purpose of using weighted sums is to account for the observational noise for which the original LS does not.
The LS method has often been used to seek dominant periodic frequencies or cycles (Zhang et al., 1993; Pichon et al., 2015; Rao et al., 2017), analyse seasonal changes in significant modulations of GW fields (Beldon and Mitchell, 2010), and estimate the spectral indices β and amplitudes of GW powerlaw spectra (Hall and Aso, 1999; Zhang et al., 2006; Guharay and Sekar, 2011; Qing et al., 2014). In addition, LS is known as the most efficient method for estimating the variance in both gapped and nongapped stationary time series with a single periodicity, without the need to fill in missing data (Marinna et al., 2019). In contrast, Vio et al. (2010) found that LS is reliable for analysing neither semiperiodic nor aperiodic signals with nonstationary noise or signals made of more than one wave, without additional steps.
2.3 The Haar structure function
HSF is a mathematical tool used in conducting scaling analysis of signals (Lovejoy and Schertzer, 2012), which is based on the Haar wavelet (Haar, 1910). It is a simple yet powerful method for decomposing a signal x(t) whose power spectral density exhibits a power law, i.e. PSD∝τ^{H} with a scaling (Hurst) exponent H, over a scale (lag) $\mathit{\tau}=\mathrm{1}/f$, into fluctuations $\mathrm{\Delta}x=x(t+\mathit{\tau})x\left(t\right)$. The firstorder Haar fluctuations ℋ_{τ} at a lag τ are defined by the following relation:
The qthorder structure function S_{q}(τ) is then obtained as an approximation by ensembleaveraging these fluctuations as follows:
In the quasiGaussian case, the moment scaling function is K(q)≈0, so for the firstorder structure function (q=1) only H determines the scaling of mean fluctuations.
A power spectrum which follows a powerlaw PSD∝τ^{β} is related to the Hurst exponent by $\mathit{\beta}=\mathrm{1}+qHK\left(q\right)$; since the power spectral density is a secondorder moment, we take q=2. Thus, under our quasiGaussian approximation, we rescaled HSF to a comparable scale to the PSD using the following relation:
Despite the term K(2) being fairly small, it is nontrivial in the nonGaussian case and for higher moments q, and HSF allows for its calculation (Lovejoy and Schertzer, 2012), but this is beyond the scope of this paper. HSF is particularly suitable for estimating the scaling exponent of time series with $H\in (\mathrm{1},\mathrm{1})$ or $\mathit{\beta}\in (\mathrm{1},\mathrm{3})$. This range of β values covers the vast majority of atmospheric processes from weather (where τ<10 d and $\mathrm{1}<\mathit{\beta}<\mathrm{3}$) to macroweather (where $\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{d}<\mathit{\tau}<\mathrm{10}$–30 years and $\mathrm{1}<\mathit{\beta}<\mathrm{1}$) systems. HSF also possesses the advantage of handling unevenly sampled data, which is a consequence of the fact that it is computed by taking the mean of absolute fluctuations (Lovejoy, 2014). Nevertheless, HSF is not employed to estimate the frequency of a wave or its amplitude, since it only measures how much frequency components contribute to the total variance. The Python code implementation of HSF is readily accessible (Mossad, 2023) and was derived from the R code originally developed by Raphaël Hébert (Hébert, 2021).
In this section, we present the simulation procedures used to generate time series similar to actual GWs measurements. In measurements, GWs can exhibit various behaviours, ranging from superposed waves within wave packets with multiple frequencies, amplitudes, and phases to more coherent quasimonochromatic waves (Maekawa et al., 1984; Eckermann and Hocking, 1989; Sica and Russell, 1999). In Sect. 3.1, we simulate signals with each having one frequency to mimic quasimonochromatic waves, by which the goal is to accurately estimate the correct frequency and amplitude of each signal. In Sect. 3.2, however, we simulate time series composed of a superposition of waves with random frequencies and powerlaw amplitudes. The exponents/slopes β of their powerlaw spectra are used to assess the bias in the spectralanalysis methods. While our simulation adopts a simplified linear saturation theory approach through the superposition of sine waves (Dewan and Good, 1986; Smith et al., 1987), we acknowledge that other explanations for the spectral character of GWs exist, including “nonlinear damping” (Weinstock, 1982, 1990; Gardner, 1994), “Doppler spreading” (Hines, 1991), “saturatedcascade similitude” (Dewan, 1994), and highReynoldsnumber “stratified turbulence” (Pinel and Lovejoy, 2014); see Table A1 for more details.
By analysing both simulations, we can determine the accuracy of each of the methods at different levels of signal complexity and identify potential limitations and sources of error in the analysis of GWs spectra. Random gaps are then introduced to resemble observational gaps for both simulations. The units and values of the variables used in this simulation have been chosen to represent average values or ranges characteristic of typical GW time series.
3.1 Singlefrequency signal simulation
Quasimonochromatic GWs can be observed under specific conditions where a single frequency dominates other components (Muraoka et al., 1988; Swenson et al., 1999). This kind of GW can be approximated as an evenly sampled single sinusoid $x\left(t\right)=A\mathrm{sin}(\mathrm{2}\mathit{\pi}ft+\mathit{\phi})$, with a known frequency f, phase shift φ, and amplitude A, at time t. For both simulations, the time resolution Δt is 5 min with a total span of T=6 h, as this resolution and this duration align with the average values of lidar measurements commonly used in atmospheric studies (Gardner et al., 1995; Gerding et al., 2008). As a result, the number of points N in each signal is equal to $\frac{T}{\mathrm{\Delta}t}=\mathrm{72}$. Each simulated sinusoid has an amplitude of A=4 K and a randomly chosen frequency f from the set $\mathrm{1}/\mathit{\{}\mathrm{6},\mathrm{3},\mathrm{1.5},\mathrm{1},\mathrm{0.5},\mathrm{1}/\mathrm{3}\mathit{\}}$ [h^{−1}]. Changing the frequency serves as a test to examine whether the bias in the methods is frequencydependent. The phase shift φ is also randomly chosen from a uniform distribution within the interval [0,2π].
The amplitude of the time series is equivalent to the estimated height of the main peak in the spectrum. It is computed from the FFT coefficients (Eq. 1) as ${A}_{\text{FFT}}={max}_{k}\left(\frac{\mathrm{2}\left{Z}_{k}\right}{N}\right)$ and from the GLS fit coefficients (Eq. 2) as ${A}_{\text{GLS}}={max}_{k}\left(\sqrt{{a}_{k}^{\mathrm{2}}+{b}_{k}^{\mathrm{2}}}\right)$. The frequency of that peak f_{k} corresponds to the estimated frequency of the signal. As a metric for the accuracy of estimation of the true values of frequencies and amplitudes, we used the relative bias given by
Since real data are susceptible to observational noise, it is crucial to consider the case where white noise is added to the simulated signal as a random variable r(t) from a standard normal distribution. Here, the signaltonoise ratio is defined by $\text{SNR}={A}^{\mathrm{2}}/\mathrm{2}{\mathit{\sigma}}_{r}^{\mathrm{2}}$, where ${\mathit{\sigma}}_{r}^{\mathrm{2}}$ is the noise variance (Horne and Baliunas, 1986). To strike an appropriate balance between capturing meaningful noise characteristics and minimizing low SNR bias, an average SNR value of 8 is chosen for this simulation.
3.2 Spectral powerlaw simulation
As reported before (see Table A1), the spectra of GWs are characterized by a power law; i.e. $\text{PSD}\propto \mathrm{1}/{f}^{\mathit{\beta}}$ (in the case of spatial data, $\text{PSD}\propto \mathrm{1}/(k,l,m{)}^{\mathit{\beta}}$, where k, l, and m are horizontal and vertical wavenumbers). On that account, we are interested in estimating the value of β of simulated time series whose spectra would have $\mathit{\beta}\in \mathit{\{}\mathrm{1},\mathrm{0},\mathrm{1},\mathrm{5}/\mathrm{3},\mathrm{2},\mathrm{2.5},\mathrm{3}\mathit{\}}$ and comparing it with the true value. The simulated evenly sampled time series x(t) consists of a sum of M sinusoids, each with frequency f_{i} and powerlaw amplitudes ${A}_{i}={f}_{i}^{\mathit{\beta}}$, as follows (Rice, 1944; Billah and Shinozuka, 1990; Kirchner, 2005):
We used for our simulation M=35 for each time series. It is driven by the objective of reconciling the inclusion of a minimum of 20 waves based on observational (Sica and Russell, 1999) and modelling (Dewan, 1994; Hamilton, 1997) suggestions while simultaneously incorporating a sufficiently large number of waves to mitigate random spectral errors (Keisler and Rhyne, 1976; Shinozuka, 2005). This approach enables the demonstration of the powerlaw spectrum without the need for excessive averaging. The phase shifts φ_{i} are also randomly chosen from a uniform distribution within the interval [0,2π]. The frequencies f_{i} are statistically independent, uniformly distributed random values, selected within the range $\left[\frac{\mathrm{1}}{T},\frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta}t}\right]=\left[\frac{\mathrm{1}}{\mathrm{6}\phantom{\rule{0.125em}{0ex}}\mathrm{h}},\frac{\mathrm{1}}{\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{min}}\right]$. Hence, the time series is composed of nonharmonic components and there are no favoured frequencies, which is a better approximation of atmospheric GWs than an idealistic case where frequencies are only integer multiples of a fundamental frequency. Here, x(t) is proportional to the square root of the amplitudes A_{i}, since β is estimated from PSDnormalized spectra, which are the squared modulus of the amplitudes.
The PSD is obtained by FFT using the relation $\frac{\mathrm{2}{Z}_{k}{Z}_{k}^{*}\mathrm{\Delta}t}{N}$, where ${Z}_{k}^{*}$ is the complex conjugate of Z_{k}. This definition is equivalent to the GLS spectrum $\frac{N\left{a}_{k}^{\mathrm{2}}+{b}_{k}^{\mathrm{2}}\right\mathrm{\Delta}t}{\mathrm{2}}$ of evenly sampled data. HSF is however normalized according to Eq. (5c) to estimate a comparable scale to the PSD. The bias in β estimation is defined for this simulation as
3.3 Gap simulation
After creating a time series with the desired spectral properties, gaps are introduced by randomly removing data points (except both endpoints), assuming that all data points are equally likely to be removed (i.e. a uniform distribution). Based on the simulated GP p in the data, an integer number of random points ${N}_{G}=\frac{N\cdot p}{\mathrm{100}}$ is removed. Thus, a 0 % GP means that no points were removed, while a 50 % GP means that 36 points were randomly removed, since N=72. To assess the dependence of bias in spectralanalysis methods on the gaps, we conducted simulation runs spanning GPs ranging from 0 % to 90 % in increments of 10 % for each time series analysed. Each of these simulation runs was repeated 1000 times at each GP increment to ensure the statistical significance of our results, since the frequencies, phases, and gaps are randomized. Then we computed the average values of the estimated amplitude and period (for the singlefrequency signal simulation) and β (for the spectral powerlaw simulation).
Before applying spectral methods to the generated time series from the simulations in Sect. 3, the following steps are taken:

The time series is first interpolated using the original time step of 5 min; this is only necessary for FFT.

The mean of the signal $\stackrel{\mathrm{\u203e}}{z}$ is subtracted to account for the zerofrequency component of the Fourier transform Z_{0}.
When computing the spectrum for a time step Δt, the frequency grids of all methods are defined as follows:

The frequency range spans from $\frac{\mathrm{1}}{T}$ to $\frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta}t}$.

The frequency spacing is given by $\mathrm{\Delta}f=\frac{\mathrm{1}}{T}$.
In the presence of gaps where Δt is not constant, the Nyquist frequency is then defined as follows:

${f}_{\text{Ny}}=\frac{\mathrm{1}}{\mathrm{2}p}$, where p is the largest value that allows ${t}_{i}={t}_{\mathrm{0}}+{n}_{i}p$ to be possible for all t_{i} values, and n_{i} comprises integers (Eyer and Bartholdi, 1999). This value corresponds to the same Nyquist frequency as that of our nongapped data, which is 0.1 min^{−1}.

This approach is more appropriate for GLS and HSF than applying a “pseudoNyquist” limit based on an average or a minimum value of Δt (Scargle, 1982; Vanderplas, 2018).
To estimate β, the spectra are fitted by taking the following steps:

A maximum likelihood estimator (MLE) is employed to determine the fit parameter β (Duvall and Harvey, 1986).

The MLE fit involves minimizing the negative loglikelihood function −ln L(O) of observations O_{i} at frequency f_{i} using the equation
$$\begin{array}{}\text{(9)}& \mathrm{ln}L\left(O\right)=\sum _{i=\mathrm{1}}^{n}\mathrm{ln}\langle {O}_{i}\rangle +{\displaystyle \frac{{O}_{i}}{\langle {O}_{i}\rangle}}\phantom{\rule{0.33em}{0ex}},\end{array}$$where 〈O_{i}〉 refers to the powerlaw model $c\left(\frac{\mathrm{1}}{{f}^{\mathit{\beta}}}\right)$ being fitted, with c as a normalization coefficient.

The MLE fit is recommended over leastsquares regression because the latter assumes a Gaussian distribution of periodogram residuals, leading to a biased estimate of β (Clauset et al., 2009).
5.1 Singlefrequency signal
First, we show an example of the time series generated by the simulation described in Sect. 3.1, which consists of a 0.5 h wave in a 6 h time series. As can be seen in Fig. 1a, in the absence of gaps, an accurate estimation of the 4 K amplitude and the 0.5 h period of this signal is acquired by both GLS and FFT. In addition, the spectra obtained by FFT and GLS are (as expected) equivalent in the case of evenly sampled data (Scargle, 1982). When the random gaps replace 50 % of the data points, a significant difference between the amplitude spectra is observed; see Fig. 1b. Both methods still provide an accurate estimate of the signal's period. Linear interpolation of the 50 % gapped signal preserves the structure of the wave but loses some of the highfrequency components, which leads to a significant underestimation of the amplitude by 43 % in the FFT spectrum. In contrast, the amplitude of the highest peak in the GLS spectrum is not affected by the gaps and has not changed from the expected value.
When comparing average relative period bias for different simulation periods, Fig. 2a shows that GLS demonstrates no period bias below 80 % GP and a negligible bias beyond (within ±20 % deviation interval). Figure 2b shows, contrarily, that the smaller the simulated period of the signal, the more FFT overestimates it at GPs larger than 40 %. This is due to linear interpolation replacing the removed highfrequency components with lowerfrequency ones, which eventually dominate as the GP gradually increases. To put these results in perspective, at 70 % GP, FFT inaccurately determines the period of a 0.5 h wave as 2.93 h (i.e. overestimates it by 486 %). Given that there are quite few data points left in the time series at 70 % GP, it is expected that FFT is not able to recover the correct period; however, GLS is still capable of obtaining the exact value of the period (i.e. an error of 0 %). Not only is GLS a much better estimator of the period on average, but also since its standard deviation remains trivially small until 80 % GP, it is a more reliable choice on a casebycase basis.
Similarly, FFT's amplitude bias experiences clear dependency on the frequency of the signal, while GLS demonstrates a negligible amplitude bias at GPs below 80 %; see Fig. 3. In particular, the results show that the mean estimated amplitude from the FFT spectrum deviates exponentially from the expected value as the frequency of the signal increases, even at gap percentages below 50 %. In addition, the standard deviation of amplitude estimation by GLS increases as the GP increases, although it remains within a ±10 % deviation interval up to 80 % GP. On the contrary, the standard deviation of amplitude estimation by FFT significantly increases as both the frequency of the signal and the GP increase, implying that FFT is more inconsistent and highly sensitive to missing data, especially for highfrequency signals.
Overall, GLS provides a more robust estimation of the period and amplitude of gapped time series, while FFT's performance is simultaneously dependent on the GP and frequency of the signal.
5.2 Spectral powerlaw signal
A timeseries example is shown in Fig. 4 to illustrate the complexity of a signal produced according to Sect. 3.2 for β=2, showcasing the signal before and after the introduction of gaps. As the percentage of gaps in the data increases, the impact on the spectral components varies depending on their frequency range. Highfrequency components (rapid fluctuations over short periods) are most vulnerable when data are being removed. Subsequently, the lowerfrequency components (slow variations over longer periods) follow, exhibiting greater resilience to data gaps. In essence, a considerably greater number of gaps are required to significantly affect the estimation of the latter. A distorted spectrum of this time series can be seen, which is a result of the limited sample length and resolution (Roberts et al., 1987). In addition, as the signal comprises numerous waves with random frequencies, the presence of closely spaced frequencies leads to the emergence of complex and broad peaks in the spectrum (Horne and Baliunas, 1986; Dewan and Grossbard, 2000).
In the absence of gaps, the true value of β is accurately estimated from the spectra of all methods (see Fig. 4a). However, after removing 50 % of the data points (Fig. 4b), the estimated spectrum by HSF remains relatively unaltered, while the GLS and FFT spectra diverge. The overestimation of β (bias=0.6) by FFT is due to the amplitudes of highfrequency components being underestimated, which is a result of the interpolation (Schulz and Stattegger, 1997; Hall and Aso, 1999). This overestimation of β can also be explained by the fact that these interpolated (highfrequency) components contribute locally by β=3 to the overall slope of the spectrum, which results in positive bias when the true β is less than 3 (Lovejoy, 2014). In the case of HSF, where no interpolation of the signal is done, an even smaller bias (0.1) results. In contrast, the true power law can be seen for the first few low frequencies in the GLS spectrum; then it starts to flatten at intermediate and high frequencies with a substantial bias of −1.41. This occurs because the lack of data points constrains the leastsquares fit by GLS, which leads to the interpolation of power at these frequencies (i.e. a flat line).
For a statistically significant picture, we show the distributions of estimated β values from 1000 simulation runs in Fig. 5. In the nongapped case (see Fig. 5a), the distributions of all methods overlap within a small standarddeviation range of ±0.2 around β=2. In contrast, we see that gaps cause the estimated β values from the FFT spectra to shift their distribution to higher values around β=2.3 in Fig. 5b, while the distribution of GLS becomes narrower and is diverted far below the expected value of 2. The mean of estimated β values from HSF spectra is clearly the closest to both the true value and its mean in the absence of gaps, which shows consistency and less sensitivity to the gaps. It is worth noting that the distributions in Fig. 5a show that even in the absence of gaps, estimated β lies mostly within the range of [1.5,2.5] and not exactly at 2, as single spectra are distorted without averaging. It is also important to mention that the results of the bias are almost identical whether β is estimated from averaging powerlaw exponents of single spectra or it is estimated from fitting an averaged spectrum.
To further explore the behaviour of the bias under different conditions, we evaluated the effect of changing the simulated value of β on the estimation bias (see Fig. 6). At 0 % GP, all methods show no bias, except for β>2 and $\mathit{\beta}=\mathrm{1}$, where there is an apparent deviation. The first deviation is expected because the more β is larger than 2, the more the spectrum suffers from “lowfrequency leakage” due to the finite length of the time series (Klis, 1994; Schulz and Mudelsee, 2002). The other case of small deviation takes place on the opposite end of the spectral slope range ($\mathit{\beta}=\mathrm{1}$), where high frequencies dominate. The power at these frequencies is quite easily aliased, hence underestimated, as they are the closest to the Nyquist limit. These deviations do not mean that spectra of $\mathit{\beta}=\mathrm{1}$ or β>2 cannot be obtained but that on average, they are very likely to be misestimated.
In the instance of gapped time series, our results show that as the GP increases, the biases in the estimated exponent also became more pronounced for all three methods. Similarly to the nongapped situation, GLS demonstrates an exceptional efficiency in estimating flat spectra where β=0 with a negligible bias. This is a consequence of the absence of frequency dependency in a flat spectrum, which renders the gaps irrelevant in terms of introducing bias, since the GLS spectrum is already flat. As β increases (indicating a steeper decline in power with increasing frequency) and the percentage of gaps in the data increases, the bias in the GLS spectrum becomes more prominent. For instance, in the case of β transitioning from 1→3 where power is skewed towards low frequencies, gaps cause GLS to mistakenly assign excessive power at the missing high frequencies, ultimately resulting in a steady underestimation of β. In contrast, when $\mathit{\beta}=\mathrm{1}$, highfrequency components dominate lowfrequency ones. Consequently, there is a loss of power at these high frequencies as the gaps disrupt their sampling. This loss of power causes the GLS method to overestimate β, mirroring the GLS bias observed when β=1.
In similar fashion, both FFT and HSF demonstrate a relatively constant bias for β=3 of approximately −0.3 at all GPs. However, as β decreases from $\mathrm{2}\to \mathrm{1}$, their biases monotonically increase as the GP increases. Nevertheless, HSF shows substantially less bias than FFT when the GP exceeds 10 %. The FFT bias is attributed to the established interpolation effects. Therefore, as more data points are interpolated, the FFT spectrum progressively underestimates the amplitude of the high frequencies. This underestimation results in the bias being positive for all β values, except β=3, where leakage causes FFT to overestimate these frequencies. Overall, averaging β values from single spectra is a good measure of the expected value because of their low standard deviation except at very high GPs.
In light of the aforementioned considerations, it can be argued that the FFT technique demonstrates competence in generating accurate spectral estimations for nongapped time series. Nevertheless, it encounters challenges as the data incorporate an increasing number of gaps, necessitating interpolation techniques which introduce inherent biases. Meanwhile, HSF is demonstrated to be a particularly reliable approach for analysing GW time series with spectral powerlaw exponents $\mathit{\beta}\in \mathit{\{}\mathrm{1},\mathrm{5}/\mathrm{3},\mathrm{2},\mathrm{2.5},\mathrm{3}\mathit{\}}$ and in between. Its performance, however, exhibits limitations primarily in cases where the spectrum of a time series possesses a powerlaw exponent β<1. Notably, such occurrences have only been observed and predicted within measurements of vertical wind time series, as indicated in Table A1. One can reasonably anticipate that spectra falling within the range of β values between 1 and 3 will be prevalent across the majority of atmospheric time series.
Conversely, the GLS method yields similarly favourable outcomes, particularly for time series whose spectra are flat and highfrequencydominated, so that it even surpasses the accuracy of HSF when β possesses values between −1 and 0. Nonetheless, the GLS method exhibits an increasing bias as the value of β increases beyond 0, rendering it a less certain choice for timescales extending beyond a few hours, which commonly occur in the context of atmospheric gravity waves. Clearly, the consistent and overall impeccable GLS performance in estimating periods and amplitudes of signals with single frequencies does not seem to translate to universally resolving the level of superposition of many random periodicities with powerlaw amplitudes.
6.1 Lowfrequency leakage
The problem of power leakage from low frequencies into higher ones arises as a result of the constrained frequency range, which itself is limited by the observed time span T. This leakage takes place not only in the case of spectra with β>2, but also in the spectra of time series with periods longer than T. GWs can often have these kinds of periods longer than the simulated time span of 6 h, and normally these periods are not resolved. Thus, we also quantitatively tested the effect of these long periods on the estimated spectra by adding three extra waves of 8, 10, and 12 h periods into the simulated time series. Two cases were examined, one where each of the three waves has an amplitude equivalent to the lowestfrequency component in each simulation. In another case, we scaled their amplitudes by the same powerlaw exponent β as all frequencies in the simulation.
In both cases, longerthanT periods produced quite similar effects on the spectra, with a substantial positive bias observed for all signals with β<2 and a significant negative bias for β>2, even in the absence of gaps. For instance, Fig. 7a and b show an example of this leakage in averaged spectra for β=3, without and with the extra waves. A common feature of both cases is the spectral power being excessively concentrated at the lowest frequencies. When the extra long waves are added, the leakage becomes more drastic for GLS and FFT, and their absolute biases of β increase. In contrast, HSF is less affected by these long periods compared to FFT and GLS. This effect of longerthanT waves in the time series resembles that of trends, which contribute to the spectral shape with a powerlaw exponent β=2 (Klis, 1994).
While a weighted fit of the spectra can reduce the bias, it does not fully rectify the problem of leakage; it also requires a smoothed spectrum and may be confounded by other biases from observational noise, gaps, or method inefficiencies. One approach to counteract this leakage from unresolved lowfrequency power is based on “prewhitening” the time series and then “postdarkening” the spectra (Blackman and Tukey, 1958). Prewhitening is a technique to decorrelate the time series (reduce its autocorrelation close to that of a white noise signal) before calculating the PSD (Keisler and Rhyne, 1976; Dewan and Grossbard, 2000; Guharay and Sekar, 2011). This step serves to decrease the rate of change in PSD(f) with frequency f; hence, the distribution of the power is more even. The prewhitened time series is obtained by the means of first differencing, where each data point is subtracted from its successive one; i.e. $\mathrm{\Delta}x\left(t\right)=x\left(t\right)x(t\mathrm{\Delta}t)$. By doing this, the contribution of the mean becomes negligible and trends are transformed into constants (Bieber et al., 1993). The spectrum of this prewhitened series Δx(t) is related to the original spectrum of x(t) by
This factor 2(1−cos (2πf_{n}Δt)) is derived from the Fourier transform of the autocorrelation function of Δx(t) and is used to compensate for the prewhitening process (Houbolt et al., 1964). This compensation postdarkens (recolours) the resulting spectrum. Note that the spectrum of the prewhitened series Δx(t) should be smoothed first using a Hann window to suppress the random fluctuations before postdarkening.
In Fig. 7c, we present the postdarkened spectra after prewhitening the time series for β=3. This approach completely cancels out the bias in all methods for both cases. This confirms the effectiveness of the prewhitening and postdarkening method in correcting the leakage problem. However, this approach is not a perfect solution, since it may introduce additional bias for less steep spectra (where β<1) which do not suffer from leakage. When the spectrum is constant in the lowfrequency part, the prewhitening process might cause it to be distorted as well.
6.2 Method selection procedure
The spectral analysis of GW timeseries data is a complex task that requires careful consideration of various factors. Based on our simulation results, we propose a flowchart (see Fig. 8) that outlines a practical guide for selecting appropriate spectral estimation methods for GW studies, taking into account the characteristics of the observed data such as their complexity and percentage of gaps in them. From the flowchart, it is clear that there are recognizable differences between the patterns of the time series of superposed waves and signals with one frequency. Even within the former classification, the antipersistent time series whose spectra have $\mathit{\beta}\in [\mathrm{0},\mathrm{1}]$ are still differentiable from those longrangedependent ones where $\mathit{\beta}\in [\mathrm{1},\mathrm{3}]$. It is also safe to say that theoretical predictions of GW spectra with $\mathit{\beta}\in [\mathrm{0},\mathrm{1}]$ exist only for measurements of vertical wind time series (see Table A1). Otherwise, $\mathit{\beta}\in [\mathrm{1},\mathrm{3}]$ spectra should be expected for the vast majority of atmospheric time series. Note that, even in the absence of gaps in the signal, caution must be taken if the estimated β is approximately equal to or less than 2. If it is, then it can be an accurate estimation or caused by one of the systematic errors such as interpolation or longerthanT variations.
The domain of observational analyses, especially in atmospheric physics, is vast and complex. The intrinsic nature of these measurements requires precision, reliability, and adaptability of analytical methods to extract meaningful insights. Observational gaps due to instrument failures or adverse conditions are also common in atmospheric physics. Our research on synthetic atmospheric gravity wave (GW) time series offers a comprehensive overview of different spectral estimation methods, shedding light on their strengths and limitations. The methods compared in this study are the fast Fourier transform (FFT) (which requires interpolation of gaps) and the generalized Lomb–Scargle method (GLS) and the Haar structure function (HSF) (both which can handle unevenly sampled data without interpolation). Their performance is assessed by evaluating whether the output spectra of simulated time series (with known a priori features and gap distributions) match the input parameters of the data, including the frequency, amplitude, and spectral slope. Building on our findings, we propose the following recommendations:
Determination of the timeseries nature. Initially, the series is assessed to determine if it represents a pattern of a signal with one frequency. If it does, the analysis proceeds directly to using the GLS method, especially if the time series contains gaps and/or demonstrates rapid variations. Rather than interpolating these gaps to use FFT, which alters the data's characteristics, GLS offers far more accurate spectral estimates.
Evaluation of data gaps. For more complex series, particularly those characterized by a superposition of waves or the powerlaw spectrum, the data are first assessed for observational gaps. The presence of such gaps profoundly affects the accuracy of traditional spectralanalysis methods like FFT. If variations are rapid (i.e. the slope of that powerlaw spectrum β is in the range $[\mathrm{1},\mathrm{0}]$), the series is again subjected to the GLS method. However, for slower variations such as the vast majority of atmospheric measurements, represented by a β value in the range of [1,3], the HSF method is consistently the least biased.
Mitigate lowfrequency leakage. When the spectra of the time series with slow variations are analysed (even in the absence of gaps), it should be investigated whether lowfrequency leakage takes place. Finite time series, very steep spectra, and waves with periods longer than the observed time span cause this kind of leakage. All tested methods, to varying extents, struggled with this leakage from unresolved lowfrequency power, leading to biases in estimating the spectral slope β. Thus we recommend cautiously applying the approach of prewhitening the time series and postdarkening the spectrum to address this problem. Prewhitening can be actualized by firstdifferencing the series, which renders its spectrum close to that of white noise. Subsequently, postdarkening the spectrum is required to compensate for the prewhitening process via division by a frequencydependent factor. Nevertheless, depending on the characteristics of the analysed time series, further steps (e.g. windowing) might be necessary to help counter other causes of leakage.
These sequential decisions highlight the importance of tailoring analysis methods to the characteristics of the data at hand. Employing a onesizefitsall method can result in biased spectra, which are especially critical in GW parameterizations in all atmospheric models. Our recommendations, grounded in rigorous research on synthetic gravity wave time series, aim to guide researchers in making informed decisions, ensuring the accuracy of spectral results and advancing our understanding of the dynamics of the atmosphere.
MM wrote the codes to conduct these simulations, analysed their results, and drafted the manuscript. IS, RW, and GB provided supervision and scientific insight and edited the text of the manuscript.
At least one of the (co)authors is a member of the editorial board of Atmospheric Measurement Techniques. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This paper is a contribution to the project W1 (Gravity Wave Parameterization for the Atmosphere) of the Collaborative Research Centre TRR 181 “Energy Transfers in Atmosphere and Ocean” and the project Analyzing the Motion of the Middle Atmosphere Using Nighttime RMRlidar Observations at the Midlatitude Station Kühlungsborn (AMUN).
This research has been supported by the Deutsche Forschungsgemeinschaft (grant nos. 274762653 and 445400792).
The publication of this article was funded by the Open Access Fund of the Leibniz Association.
This paper was edited by Lars Hoffmann and reviewed by two anonymous referees.
Alexander, M., Tsuda, T., and Vincent, R. N.: Latitudinal Variations Observed in Gravity Waves with Short Vertical Wavelengths, J. Atmos. Sci., 59, 1394–1404, https://doi.org/10.1175/15200469(2002)059<1394:LVOIGW>2.0.CO;2, 2002. a
Axford, D.: Spectral analysis of an aircraft observation of gravity waves, Q. J. Roy. Meteor. Soc., 97, 313–321, https://doi.org/10.1002/qj.49709741306, 1971. a
Babu, P. and Stoica, P.: Spectral analysis of nonuniformly sampled data – a review, Digit. Signal Process., 20, 359–378, https://doi.org/10.1016/j.dsp.2009.06.019, 2010. a
Beldon, C. L. and Mitchell, N. J.: Gravity wave–tidal interactions in the mesosphere and lower thermosphere over Rothera, Antarctica (68^{∘} S, 68^{∘} W), J. Geophys. Res., 115, D18101, https://doi.org/10.1029/2009jd013617, 2010. a
Beres, J., Garcia, R. R., Boville, B. A., and Sassi, F.: Implementation of a gravity wave source spectrum parameterization dependent on the properties of convection in the Whole Atmosphere Community Climate Model (WACCM), J. Geophys. Res., 110, D10108, https://doi.org/10.1029/2004jd005504, 2005. a
Bieber, J. W., Chen, J., Matthaeus, W. H., Smith, C. W., and Pomerantz, M. A.: Longterm variations of interplanetary magnetic field spectra with implications for cosmic ray modulation, J. Geophys. Res., 98, 3585–3603, https://doi.org/10.1029/92ja02566, 1993. a
Billah, K. Y. R. and Shinozuka, M.: Numerical method for colorednoise generation and its application to a bistable system, Phys. Rev. A, 42, 7492–7495, https://doi.org/10.1103/PhysRevA.42.7492, 1990. a
Blackman, R. B. and Tukey, J. W.: The measurement of power spectra from the point of view of communications engineering – Part I, Bell Syst. Tech. J., 37, 185–282, https://doi.org/10.1002/j.15387305.1958.tb03874.x, 1958. a
Brown, T. M. and ChristensenDalsgaard, J.: A technique for estimating complicated power spectra from time series with gaps, Astrophys. J., 349, 667, https://doi.org/10.1086/168354, 1990. a
Clauset, A., Shalizi, C. R., and Newman, M.: Powerlaw distributions in empirical data, Siam Rev., 51, 661–703, https://doi.org/10.1137/070710111, 2009. a
Cooley, J. and Tukey, J. W.: An algorithm for the machine calculation of complex Fourier series, Math. Comput., 19, 297–301, https://doi.org/10.1090/s00255718196501785861, 1965. a, b
Crowley, G. and Williams, P.: Observations of the source and propagation of atmospheric gravity waves, Nature, 328, 231–233, https://doi.org/10.1038/328231a0, 1987. a
Dewan, E.: The saturatedcascade model for atmospheric gravity wave spectra, and the wavelengthperiod (WP) relations, Geophys. Res. Lett., 21, 817–820, https://doi.org/10.1029/94gl00702, 1994. a, b, c, d, e, f, g
Dewan, E. M. and Good, R. E.: Saturation and the “universal” spectrum for vertical profiles of horizontal scalar winds in the atmosphere, J. Geophys. Res., 91, 2742, https://doi.org/10.1029/jd091id02p02742, 1986. a, b, c
Dewan, E. M. and Grossbard, N.: Power spectral artifacts in published balloon data and implications regarding saturated gravity wave theories, J. Geophys. Res., 105, 4667–4683, https://doi.org/10.1029/1999jd901108, 2000. a, b, c
Duvall, T. L. and Harvey, J. T.: Solar Doppler shifts: sources of continuous spectra, Springer, https://doi.org/10.1007/9789400946088_11, 1986. a
Eckermann, S. D. and Hocking, W. K.: Effect of superposition on measurements of atmospheric gravity waves: A cautionary note and some reinterpretations, J. Geophys. Res., 94, 6333, https://doi.org/10.1029/jd094id05p06333, 1989. a, b
Ern, M., Trinh, Q. T., Preusse, P., Gille, J. C., Mlynczak, M. G., Russell III, J. M., and Riese, M.: GRACILE: a comprehensive climatology of atmospheric gravity wave parameters based on satellite limb soundings, Earth Syst. Sci. Data, 10, 857–892, https://doi.org/10.5194/essd108572018, 2018. a
Eyer, L. and Bartholdi, P.: Variable stars: Which Nyquist frequency?, Astron. Astrophys. Sup., 135, 1–3, https://doi.org/10.1051/aas:1999102, 1999. a
FerrazMello, S.: Estimation of Periods from Unequally Spaced Observations, Astron. J., 86, 619, https://doi.org/10.1086/112924, 1981. a
Fritts, D. C.: A review of gravity wave saturation processes, effects, and variability in the middle atmosphere, Pure Appl. Geophys., 130, 343–371, https://doi.org/10.1007/bf00874464, 1989. a
Fritts, D. C. and Alexander, M.: Gravity wave dynamics and effects in the middle atmosphere, Rev. Geophys., 41, 1003, https://doi.org/10.1029/2001rg000106, 2003. a
Fritts, D. C. and VanZandt, T. E.: Spectral Estimates of gravity wave energy and momentum fluxes. Part I: Energy dissipation, acceleration, and constraints, J. Atmos. Sci., 50, 3685–3694, https://doi.org/10.1175/15200469(1993)050<3685:SEOGWE>2.0.CO;2, 1993. a
Fritts, D. C., Tsuda, T., Sato, T., Fukao, S., and Kato, S.: Observational evidence of a saturated gravity wave spectrum in the troposphere and lower stratosphere, J. Atmos. Sci., 45, 1741–1759, https://doi.org/10.1175/15200469(1988)045<1741:OEOASG>2.0.CO;2, 1988. a
Gardner, C. S.: Diffusive filtering theory of gravity wave spectra in the atmosphere, J. Geophys. Res., 99, 20601, https://doi.org/10.1029/94jd00819, 1994. a, b, c, d, e, f
Gardner, C. S., Tao, X., and Papen, G. C.: Simultaneous lidar observations of vertical wind, temperature, and density profiles in the upper mesophere: Evidence for nonseparability of atmospheric perturbation spectra, Geophys. Res. Lett., 22, 2877–2880, https://doi.org/10.1029/95gl02783, 1995. a, b, c, d, e, f, g, h
Gerding, M., Höffner, J., Lautenbach, J., Rauthe, M., and Lübken, F.J.: Seasonal variation of nocturnal temperatures between 1 and 105 km altitude at 54^{∘} N observed by lidar, Atmos. Chem. Phys., 8, 7465–7482, https://doi.org/10.5194/acp874652008, 2008. a
Guharay, A. and Sekar, R.: Seasonal characteristics of gravity waves in the middle atmosphere over Gadanki using Rayleigh lidar observations, J. Atmos. Sol.Terr. Phy., 73, 1762–1770, https://doi.org/10.1016/j.jastp.2011.04.013, 2011. a, b, c
Haar, A.: On the theory of orthogonal function systems, Math. Ann., 69, 331–371, 1910. a
Hall, C. M. and Aso, T.: Mesospheric velocities and buoyancy subrange spectral slopes determined over Svalbard by ESR, Geophys. Res. Lett., 26, 1685–1688, https://doi.org/10.1029/1999gl900340, 1999. a, b, c, d
Hamilton, K.: The role of parameterized drag in a tropospherestratospheremesosphere general circulation model, Springer, https://doi.org/10.1007/9783642606540_23, 1997. a
He, Y., Sheng, Z., and He, M.: Spectral Analysis of Gravity Waves from Near Space HighResolution Balloon Data in Northwest China, Atmosphere, 11, 133, https://doi.org/10.3390/atmos11020133, 2020. a
Hébert, R.: RScaling (1.0.0), Zenodo, https://doi.org/10.5281/zenodo.5037581, 2021. a
Hébert, R., Rehfeld, K., and Laepple, T.: Comparing estimation techniques for temporal scaling in palaeoclimate time series, Nonlin. Processes Geophys., 28, 311–328, https://doi.org/10.5194/npg283112021, 2021. a
Hertzog, A. and Vial, F.: A study of the dynamics of the equatorial lower stratosphere by use of ultralongduration balloons: 2. Gravity waves, J. Geophys. Res., 106, 22745–22761, https://doi.org/10.1029/2000jd000242, 2001. a, b
Hines, C. O.: Internal atmospheric gravity waves at ionospheric heights, Can. J. Phys., 38, 1441–1481, https://doi.org/10.1139/p60150, 1960. a
Hines, C. O.: The Saturation of Gravity Waves in the Middle Atmosphere. Part II: Development Of DopplerSpread Theory, J. Atmos. Sci., 48, 1361–1379, https://doi.org/10.1175/15200469(1991)048<1361:TSOGWI>2.0.CO;2, 1991. a, b, c, d
Holton, J. M.: The influence of gravity wave breaking on the general circulation of the middle atmosphere, J. Atmos. Sci., 40, 2497–2507, https://doi.org/10.1175/15200469(1983)040<2497:TIOGWB>2.0.CO;2, 1983. a
Horne, J. H. and Baliunas, S. L.: A prescription for period analysis of unevenly sampled time series, Astrophys. J., 302, 757, https://doi.org/10.1086/164037, 1986. a, b
Houbolt, J. C., Steiner, R., and Pratt, K. G.: Dynamic response of airplanes to atmospheric turbulence including flight data on input and response, vol. 199, National Aeronautics and Space Administration, https://hdl.handle.net/2027/uiug.30112106585620 (last access: 25 January 2024), 1964. a
Houchi, K., Stoffelen, A., Marseille, G.J., and De Kloe, J.: Comparison of wind and wind shear climatologies derived from highresolution radiosondes and the ECMWF model, J. Geophys. Res., 115, D22123, https://doi.org/10.1029/2009jd013196, 2010. a
Keisler, S. R. and Rhyne, R. H.: An assessment of prewhitening in estimating power spectra of atmospheric turbulence at long wavelengths, Tech. rep., NASA Langley Research Center, https://ntrs.nasa.gov/api/citations/19770005667/downloads/19770005667.pdf (last access: 25 January 2024), 1976. a, b
Kirchner, J. W.: Aliasing in $\mathrm{1}/{f}^{\mathit{\alpha}}$ noise spectra: Origins, consequences, and remedies, Phys. Rev. E, 71, 066110, https://doi.org/10.1103/PhysRevE.71.066110, 2005. a
Klis, M.: Rapid variability in xray binaries–towards a unified description, NATO Science Series C, 1995 edn., Springer, Dordrecht, Netherlands, https://hdl.handle.net/11245/1.421015 (last access: 25 January 2024), 1994. a, b, c
Lepot, M., Aubin, J.B., and Clemens, F.: Interpolation in Time Series: An introductive overview of existing methods, their performance criteria and uncertainty assessment, Water, 9, 796, https://doi.org/10.3390/w9100796, 2017. a
Lindgren, E., Sheshadri, A., Podglajen, A., and Carver, R. W.: Seasonal and latitudinal variability of the gravity wave spectrum in the lower stratosphere, J. Geophys. Res.Atmos., 125, e2020JD03285, https://doi.org/10.1029/2020jd032850, 2020. a
Lindzen, R. S.: Turbulence and stress owing to gravity wave and tidal breakdown, J. Geophys. Res., 86, 9707, https://doi.org/10.1029/jc086ic10p09707, 1981. a
Liu, H., McInerney, J., Santos, S. A., Lauritzen, P. H., Taylor, M. J., and Pedatella, N.: Gravity waves simulated by highresolution Whole Atmosphere Community Climate Model, Geophys. Res. Lett., 41, 9106–9112, https://doi.org/10.1002/2014gl062468, 2014. a
Lomb, N. R.: Leastsquares frequency analysis of unequally spaced data, Astrophys. Space Sci., 39, 447–462, https://doi.org/10.1007/bf00648343, 1976. a
Lovejoy, S.: A voyage through scales, a missing quadrillion and why the climate is not what you expect, Clim. Dynam., 44, 3187–3210, https://doi.org/10.1007/s0038201423240, 2014. a, b
Lovejoy, S. and Schertzer, D.: Haar wavelets, fluctuations and structure functions: convenient choices for geophysics, Nonlin. Processes Geophys., 19, 513–527, https://doi.org/10.5194/npg195132012, 2012. a, b, c
Maekawa, Y., Fukao, S., Sato, T., Kato, S., and Woodman, R. F.: Internal inertia–gravity waves in the tropical lower stratosphere observed by the Arecibo Radar, J. Atmos. Sci., 41, 2359–2367, https://doi.org/10.1175/15200469(1984)041<2359:IIWITT>2.0.CO;2, 1984. a, b
Marinna, A. M., Alfredo, L. A., and Christopher, R. S.: Using the LombScargle method for wave statistics from gappy time series, IEEE Conference Proceedings, 2019, 1–9, https://doi.org/10.1109/CWTM43797.2019.8955285, 2019. a
Meisel, D. D.: Fourier transforms of data sampled at unequal observational intervals, Astron. J., 83, 538, https://doi.org/10.1086/112233, 1978. a
Mossad, M.: GapWaveSpectra (1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.8136556, 2023. a, b
Munteanu, C., Negrea, C., Echim, M., and Mursula, K.: Effect of data gaps: comparison of different spectral analysis methods, Ann. Geophys., 34, 437–449, https://doi.org/10.5194/angeo344372016, 2016. a, b
Muraoka, Y., Sugiyama, T., Kawahira, K., Sato, T., Tsuda, T., Fukao, S., and Kato, S.: Cause of a monochromatic inertiagravity wave breaking observed by the MU radar, Geophys. Res. Lett., 15, 1349–1352, https://doi.org/10.1029/gl015i012p01349, 1988. a
Narendra Babu, A., Kishore Kumar, K., Kishore Kumar, G., Venkat Ratnam, M., Vijaya Bhaskara Rao, S., and Narayana Rao, D.: Longterm MST radar observations of vertical wave number spectra of gravity waves in the tropical troposphere over Gadanki (13.5^{∘} N, 79.2^{∘} E): comparison with model spectra, Ann. Geophys., 26, 1671–1680, https://doi.org/10.5194/angeo2616712008, 2008. a
Nastrom, G. D., Van Zandt, T. E., and Warnock, J. M.: Vertical wavenumber spectra of wind and temperature from highresolution balloon soundings over Illinois, J. Geophys. Res., 102, 6685–6701, https://doi.org/10.1029/96jd03784, 1997. a
Pichon, A. L., Assink, J., Heinrich, P., Blanc, E., CharltonPerez, A., Lee, C. S., Keckhut, P., Hauchecorne, A., Rüfenacht, R., Kämpfer, N., Drob, D. P., Smets, P., Evers, L., Ceranna, L., Pilger, C., Ross, O. A., and Claud, C.: Comparison of colocated independent groundbased middle atmospheric wind and temperature measurements with numerical weather prediction models, J. Geophys. Res.Atmos., 120, 8318–8331, https://doi.org/10.1002/2015jd023273, 2015. a
Pinel, J. and Lovejoy, S.: Atmospheric waves as scaling, turbulent phenomena, Atmos. Chem. Phys., 14, 3195–3210, https://doi.org/10.5194/acp1431952014, 2014. a
Podglajen, A., Hertzog, A., Plougonven, R., and Legras, B.: Lagrangian temperature and vertical velocity fluctuations due to gravity waves in the lower stratosphere, Geophys. Res. Lett., 43, 3543–3553, https://doi.org/10.1002/2016gl068148, 2016. a
Qing, H., Zhou, C., Zhao, Z., Chen, G., Ni, B., Gu, X., Yang, G., and Zhang, Y.: A statistical study of inertia gravity waves in the troposphere based on the measurements of Wuhan Atmosphere Radio Exploration (WARE) radar, J. Geophys. Res.Atmos., 119, 3701–3714, https://doi.org/10.1002/2013jd020684, 2014. a, b
Rao, N. K., Ratnam, M. V., Vedavathi, C., Tsuda, T., Murthy, B. V. K., Sathishkumar, S., Gurubaran, S., Kumar, K. S., Subrahmanyam, K. V., and Rao, S. V.: Seasonal, interannual and solar cycle variability of the quasi two day wave in the lowlatitude mesosphere and lower thermosphere, J. Atmos. Sol.Terr. Phy., 152–153, 20–29, https://doi.org/10.1016/j.jastp.2016.11.005, 2017. a
Rice, S. O.: Mathematical analysis of random noise, Bell Syst. Tech. J., 23, 282–332, https://doi.org/10.1002/j.15387305.1944.tb00874.x, 1944. a
Rigling, B. D.: Application of temporal gap filling to natural power law spectrums, IEEE Geosci. Remote S., 9, 624–628, https://doi.org/10.1109/lgrs.2011.2177062, 2012. a
Roberts, D. A., Lehar, J., and Dreher, J.: Time Series Analysis with Clean – Part One – Derivation of a Spectrum, Astron. J., 93, 968, https://doi.org/10.1086/114383, 1987. a
Scargle, J. D.: Studies in astronomical time series analysis. II – Statistical aspects of spectral analysis of unevenly spaced data, Astrophys. J., 263, 835, https://doi.org/10.1086/160554, 1982. a, b, c
Schulz, M. and Mudelsee, M.: REDFIT: estimating rednoise spectra directly from unevenly spaced paleoclimatic time series, Comput. Geosci., 28, 421–426, https://doi.org/10.1016/s00983004(01)000449, 2002. a
Schulz, M. and Stattegger, K.: Spectrum: spectral analysis of unevenly spaced paleoclimatic time series, Comput. Geosci., 23, 929–945, https://doi.org/10.1016/s00983004(97)000873, 1997. a, b
Schuster, A.: On the investigation of hidden periodicities with application to a supposed 26 d period of meteorological phenomena, J. Geophys. Res., 3, 13, https://doi.org/10.1029/tm003i001p00013, 1898. a
Shibata, T., Ichimori, S., Narikiyo, T., and Maeda, M.: Spectral analysis of vertical temperature profiles observed by a lidar in the upper stratosphere and the lower mesosphere, J. Meteorol. Soc. Jpn. Ser. II, 66, 1001–1005, https://doi.org/10.2151/jmsj1965.66.6_1001, 1988. a
Shinozuka, M.: Simulation of Multivariate and Multidimensional Random Processes, J. Acoust. Soc. Am., 49, 357–368, https://doi.org/10.1121/1.1912338, 2005. a
Sica, R. J. and Russell, A. G.: How many waves are in the gravity wave spectrum?, Geophys. Res. Lett., 26, 3617–3620, https://doi.org/10.1029/1999gl003683, 1999. a, b
Smith, A. M.: Global Dynamics of the MLT, Surv. Geophys., 33, 1177–1230, https://doi.org/10.1007/s1071201291969, 2012. a
Smith, S. M., Fritts, D. C., and VanZandt, T. E.: Evidence for a saturated spectrum of atmospheric gravity waves, J. Atmos. Sci., 44, 1404–1410, https://doi.org/10.1175/15200469(1987)044<1404:EFASSO>2.0.CO;2, 1987. a, b, c
Song, I.S. and Chun, H.Y.: A Lagrangian Spectral Parameterization of Gravity Wave Drag Induced by Cumulus Convection, J. Atmos. Sci., 65, 1204–1224, https://doi.org/10.1175/2007jas2369.1, 2008. a
Swenson, G. R., Haque, R., Yang, W., and Gardner, C. S.: Momentum and energy fluxes of monochromatic gravity waves observed by an OH imager at Starfire Optical Range, New Mexico, J. Geophys. Res., 104, 6067–6080, https://doi.org/10.1029/1998jd200080, 1999. a
Vanderplas, J.: Understanding the lomb–scargle periodogram, Astrophys. J. Suppl. S., 236, 16, https://doi.org/10.3847/15384365/aab766, 2018. a, b
VanZandt, T.: A universal spectrum of buoyancy waves in the atmosphere, Geophys. Res. Lett., 9, 575–578, https://doi.org/10.1029/gl009i005p00575, 1982. a, b
Vio, R., Andreani, P., and Biggs, A. D.: Unevenlysampled signals: a general formalism for the LombScargle periodogram, Astron. Astrophys., 519, A85, https://doi.org/10.1051/00046361/201014079, 2010. a
Weinstock, J.: Nonlinear Theory of Gravity Waves: Momentum Deposition, Generalized Rayleigh Friction, and Diffusion, J. Atmos. Sci., 39, 1698–1710, https://doi.org/10.1175/15200469(1982)039<1698:NTOGWM>2.0.CO;2, 1982. a
Weinstock, J.: Saturated and unsaturated spectra of gravity waves and scaledependent diffusion, J. Atmos. Sci., 47, 2211–2226, https://doi.org/10.1175/15200469(1990)047<2211:SAUSOG>2.0.CO;2, 1990. a, b, c, d
Zechmeister, M. and Kürster, M.: The generalised LombScargle periodogram. A new formalism for the floatingmean and Keplerian periodograms, Astron. Astrophys., 496, 577–584, 2009. a, b
Zhan, Q., Manson, A. H., and Meek, C.: The impact of gaps and spectral methods on the spectral slope of the middle atmospheric wind, J. Atmos. Terr. Phys., 58, 1329–1336, https://doi.org/10.1016/00219169(95)00159x, 1996. a
Zhang, S., Huang, C., Huang, K., Gong, Y., Chen, G., Gan, Q., and Zhang, Y.: Latitudinal and seasonal variations of vertical wave number spectra of threedimensional winds revealed by radiosonde observations, J. Geophys. Res.Atmos., 122, 13174–13190, https://doi.org/10.1002/2017jd027602, 2017. a, b, c, d
Zhang, S. D., Huang, C., and Yi, F.: Radiosonde observations of vertical wave number spectra for gravity waves in the lower atmosphere over Central China, Ann. Geophys., 24, 3257–3265, https://doi.org/10.5194/angeo2432572006, 2006. a, b, c
Zhang, S.N., Peterson, R. N., Wiens, R. H., and Shepherd, G. G.: Gravity waves from O_{2} nightglow during the AIDA '89 campaign I: emission rate/temperature observations, J. Atmos. Terr. Phys., 55, 355–375, https://doi.org/10.1016/00219169(93)900749, 1993. a
Zhu, X.: A new theory of the saturated gravity wave spectrum for the middle atmosphere, J. Atmos. Sci., 51, 3615–3626, https://doi.org/10.1175/15200469(1994)051<3615:ANTOTS>2.0.CO;2, 1994. a, b