Equatorial Wave Expansion of Instantaneous Flows for Diagnosis of Equatorial Waves from Data: Formulation and Illustration

Rossby waves, named after the seminal work of (Rossby, 1939), play a fundamental role in large-scale atmospheric and oceanic motions that affect weather and climate. An equally seminal work on Rossby waves is (Yeh, 1949), who discovered the downstream development mechanism of weather systems in the midlatitudes via the dispersion relation of Rossby waves. In the extratropics, more than 90% of wave fields at any given instance are in the form of Rossby waves. In the equatorial region, on the other hand, all types of waves can co-exist for a relatively longer time after their excitations owing to the much longer inertial time scale for horizontal motions near the equator. A theoretical as well as practical question is whether we can isolate these individual equatorial waves from the total wave flow field observed at a given instance, which is the subject of this study.

The study of equatorial waves has been extensive following the seminal analysis of free wave solutions to the linearized, equatorial shallow-water system on a β-plane by (Matsuno, 1966). (Longuet-Higgins, 1968) later showed that this is a limiting case of the more general set of Laplace's tidal equations on a sphere, with the vertical separation parameter (equivalent depth) sufficiently small. Thereafter, waves matching the characteristics of the theoretical solutions were discovered in observations (Yanai and Maruyama, 1966; Maruyama and Yanai, 1967; Wallace and Kousky, 1968a, 1968b), which was influential in the development of a theory for the mechanism behind the newly-discovered quasi-biennial oscillation (QBO) of the zonal flow in the equatorial stratosphere (Lindzen and Holton, 1968; Holton and Lindzen, 1972). Studies have shown that these waves are generated by, and coupled to, convection in the tropical lower atmosphere (Wheeler and Kiladis, 1999; Wheeler et al., 2000; Straub and Kiladis, 2003; Kiladis et al., 2009), and they are a dominant instrument of momentum and mass transport vertically into the stratosphere and meridionally into the sub- and extratropics (Baldwin and Tung, 1994; Choi et al., 2002; Giorgetta et al., 2002, 2006; Kawatani et al., 2010). In addition to being an integral part of the QBO evolution, equatorial Kelvin and Rossby waves are also thought to be equally important to the Madden-Julian oscillation (MJO), although the exact interaction is unclear (Zhang, 2005; Dias et al., 2013a).

Diagnostics of reanalysis datasets and outgoing longwave radiation measurements have reinforced the idea that equatorial waves in the real atmosphere closely match the waves that are free solutions to atmospheric models. This has in turn led to the development of several diagnostic techniques for identifying and studying waves in observations, mainly by searching for the associated equatorial wave structures in physical or spectral space in the data. As stated in the opening paragraph, we are most interested in diagnosing instantaneous wave fields so that we may study the temporal evolution of equatorial waves and their interactions with tropical convections. This is particularly useful for understanding interactions between waves and equatorial phenomena, such as the QBO and MJO, and the temporal evolution of equatorial waves, as they are continuously generated and propagate through changing background states. A method that decomposes a full wavenumber-frequency spectrum into all suitable wave classes and produces instantaneous wave fields without prescription of any a priori zonal or temporal scale could be a powerful tool for analysis of equatorial waves. Ideally, one could identify which waves are present, their locations, and their strengths, given only a snapshot of the atmospheric state.

The horizontal and vertical structures of equatorial waves are separable and may be modeled using solutions to the shallow-water equations (e.g., Matsuno, 1966) and Laplace's tidal equations (e.g., Lindzen, 1967; Longuet-Higgins, 1968), respectively, with the separation constant being the equivalent depth h_e. In general, only expansion in the horizontal direction is performed, despite knowing the vertical structure functions, because the available vertical resolution of most datasets providing dynamical variables is unsatisfactory relative to the vertical scale of the waves. As such, the shallow-water framework alone is primarily invoked for equatorial wave diagnosis, as will be done in this paper.

There are two popular varieties of methods that employ the shallow-water wave solutions to study equatorial waves. Spectral analysis techniques (e.g., Wheeler and Kiladis, 1999; Tindall et al., 2006a; Lott et al., 2014) are statistical diagnoses of the average wave properties in a dataset. This method does not explicitly identify instantaneous equatorial wave fields, except through application of a spectral filter to capture only the signal around a given dispersion relation. The evolution of wave fields diagnosed in this manner implicitly derive their spectral properties from the spatial and temporal scales chosen by the filter, which presupposes that the spectral characteristics closely match the dispersion curve for a static range of equivalent depths. However, the theoretical dispersion relation may be an inaccurate representation of the dispersion relation in the true atmosphere, due to Doppler shifting by non-zero mean zonal flow, heating, or other nonlinear effects (Yang et al., 2003). If waves additionally experienced an evolution in the local equivalent depth, their physical and spectral characteristics would be adjusted accordingly (Dias et al., 2013b).

Instead of assuming the dispersion relation and finding a physical field through spatiotemporal filtering, meridional mode decomposition (Yang et al., 2003) considers the horizontal field structure of individual equatorial waves following the analytic free solutions to the shallow-water system on an equatorial β-plane (Matsuno, 1966). By its nature, this type of diagnostic technique lends itself more to instantaneous identification of wave fields than the spectral analysis method. The frequency spectra and propagation characteristics are revealed simply from the time series of consecutive instantaneous fields, rather than from a prescribed timescale. These solutions are constructed using the basis set of parabolic cylinder functions (PCFs), although the shallow-water equatorial wave structures themselves are not necessarily orthogonal because they are linear combinations of PCFs, where each PCF projects onto multiple wave modes (Castanheira and Marques, 2015). This method relies on an initial spatiotemporal separation of the total field into eastward- and westward-propagating components, as well as a low-pass filter to remove inertia-gravity wave classes. Furthermore, EOF analysis (e.g., Fraedrich et al., 1997) is another method by which equatorial wave classes can be diagnosed in physical space, although the waves can only be named (if desired) by invoking the shallow-water framework regardless. Given the body of work on equatorial waves, we believe there is little question that the Matsuno modes reasonably represent the significant wave modes in the equatorial atmosphere.

Here, we propose a generalization of the meridional mode decomposition method by which we may isolate individual equatorial wave fields instantaneously given only a snapshot of the atmospheric state. As with current diagnostic techniques, we will employ the shallow-water framework, circumventing the non-orthogonality of the theoretical equatorial wave structures by synthesizing horizontal fields of each of the dynamical variables u, v, and φ, where u and v are zonal and meridional winds and φ is the geopotential field. This method endeavors to diagnose individual equatorial waves in a given instantaneous flow field using the equatorial wave structures themselves as a basis set. We begin by briefly describing the method by which the shallow-water framework may be utilized to diagnose the full equatorial wave spectrum in observations without spatiotemporal filtering of the total field. Following this, we illustrate the effectiveness of the equatorial wave expansion method by applying it to idealized data, shallow-water model output, and observation data.

To facilitate an easy presentation, all variables and equations in the remaining part of the paper have been non-dimensionalized by using \(\sqrt{1/(c\beta)}\) for the time scale and \(\sqrt{c/\beta}\) for the length scale, where β is the beta parameter evaluated at the equator and \(c=\sqrt{gh_e}\), where g is gravity and h_e is the equivalent depth of the shallow-water system. A convenient set of meridional basis functions in the equatorial domain is that of the PCFs, ψ_n(y), centered at the equator. Any horizontal instantaneous flow field on a discrete grid may be expanded as \begin{equation} \lambda'(x,y)=\sum_{k=-K}^K\sum_{n=0}^N\chi_{k,n}^\lambda \psi_n(y)e^{ikx} , \ \ (1)\end{equation} where Λ is a dummy flow variable, which can be zonal (u) or meridional (v) wind or geopotential height (φ), and the prime denotes the departure field from its zonal mean. Here, zonal wavenumber k runs from -K to K, the Nyquist sampling wavenumber. The meridional mode number n is a non-negative integer, and N is the largest meridional mode considered, either by choice or constraint due to meridional resolution. By projecting a flow field Λ'(x,y) onto ψ_n(y)e^ikx, the projection coefficient χ_k,n^Λ may be easily retrieved.

The same PCFs also appear in the construction of free solutions to the linearized equatorial β-plane shallow-water equations associated with the so-called "Matsuno modes" (Matsuno, 1966; Longuet-Higgins, 1968). Each meridional mode supports at most three distinct wave classes, whose frequencies are the roots of the equatorial dispersion relation for that mode number n. As a result, we can alternatively expand an instantaneous flow using the equatorial wave solutions, the Matsuno modes, as basis functions. The dispersion relation for equatorial waves is \begin{equation} \omega^2-k^2-\dfrac{k}{\omega}=2n+1 . \ \ (2)\end{equation} The dispersion relation yields three distinct real roots for each pair of zonal wavenumber k and meridional mode number n (n≥ 0), corresponding to a pair of eastward- and westward-propagating mixed Rossby-gravity waves for n=0 or a triplet of eastward- and westward-propagating inertia-gravity waves and a westward-propagating Rossby wave (n>0). Note that for n=0, the third root does not correspond to a physical solution since its eigenvector does not satisfy the boundary conditions. Also, Eq. (2) has one special root for n=-1; namely, ω=k, corresponding to a (eastward-propagating) Kelvin wave. As in (Matsuno, 1966), we use ω_k,n^(j) with the root index j=1,2,3 to denote these roots, following the convention j=1 for Kelvin (K), eastward-propagating mixed Rossby-gravity (EMRG) or eastward-propagating inertia-gravity (EIG) waves; j=2 for westward-propagating mixed Rossby-gravity (WMRG) or westward-propagating inertia-gravity (WIG) waves; and j=3 for westward-propagating Rossby (R) waves. Then, the Matsuno modes of (n,j,k) can be written as \begin{align} U_{k,n}^{(j)}(x,y)&=\dfrac{i}{2}[A_{k,n}^{(j)}\psi_{n+1}(y)-B_{k,n}^{(j)}\psi_{n-1}(y)]e^{ikx} ,\ \ (3a)\\ V_{k,n}^{(j)}(x,y)&=\psi_n(y)e^{ikx} ,\ \ (3b)\\ \Phi_{k,n}^{(j)}(x,y)&=\dfrac{i}{2}[A_{k,n}^{(j)}\psi_{n+1}(y)+B_{k,n}^{(j)}\psi_{n-1}(y)]e^{ikx} , \ \ (3c)\end{align} where $$ \left[ \begin{array}{l@{\quad}l} A_{k,n}^{(j)}=\dfrac{\sqrt{2(n+1)}}{\omega_{k,n}^{(j)}-k} & {\rm for}ñ\geq 0\\ A_{k,n}^{(1)}=1 & {\rm for}ñ=-1 \end{array} \right. $$ and \begin{equation} B_{k,n}^{(j)}=-\dfrac{\sqrt{2n}}{\omega_{k,n}^{(j)}+k}\quad {\rm for}ñ\geq 0 . \ \ (4)\end{equation} Note that in equations (3a), (3b), and (3c), U_k,n^(j)(x,y), V_k,n^(j)(x,y), and Φ_k,n^(j)(x,y) are, respectively, zonal velocity, meridional velocity, and geopotential height fields of a Matsuno mode (n, j, k). B_k,n^(j) do not exist (or are zero) for all n<1 and A_k,-1^(j) exists only for j=1, and all other A_k,-1^(j) are undefined (or zero) since there are no roots at n=-1 other than ω_k,-1⁽¹⁾=k for Kelvin waves. The fields defined in Eq. (3) for each (n,j) pair of a given k serve as a set of basis functions for their respective variables; namely, \begin{align} u'(x,y)&=\sum_{k=-K}^K\sum_{n=-1}^N\sum_{j=1}^3\alpha_{k,n}^{(j)}U_{k,n}^{(j)}(x,y) ,\ \ (5a)\\ v'(x,y)&=\sum_{k=-K}^K\sum_{n=-1}^N\sum_{j=1}^3\alpha_{k,n}^{(j)}V_{k,n}^{(j)}(x,y) ,\ \ (5b)\\ \phi'(x,y)&=\sum_{k=-K}^K\sum_{n=-1}^N\sum_{j=1}^3\alpha_{k,n}^{(j)}\Phi_{k,n}^{(j)}(x,y) , \ \ (5c)\end{align} where α_k,n^(j) is the coefficient of the (n,j) equatorial wave class for zonal wavenumber k. Because the PCFs, ψ_n(y), exist only for non-negative n, Eq. (3) indicates that the Matsuno mode for Kelvin waves is represented by U_k,-1⁽¹⁾ and Φ_k,-1⁽¹⁾ (i.e., v'=0 for Kelvin waves).

We next equate the mathematical construction of the total wave field, Eq. (1), where we use the PCFs as a basis, with the physical construction, Eq. (5), where we write the wave fields using the Matsuno modes as the basis set. After some laborious but otherwise straightforward manipulations, we have \begin{eqnarray} 2A_{k,-1}^{(1)}\alpha_{k,-1}^{(1)}&=&\chi_{k,0}^{\phi}+\chi_{k,0}^u , \ \ (6)\\ \left[ \begin{array}{c@{\quad}c} 2A_{k,1}^{(1)} & 2A_{k,1}^{(2)}\\ 1 & 1 \end{array} \right]\left[ \begin{array}{c} \alpha_{k,0}^{(1)}\\ \alpha_{k,0}^{(2)} \end{array} \right]&=&\left[ \begin{array}{c} \chi_{k,1}^{\phi}+\chi_{k,1}^u\\ \chi_{k,0}^v \end{array} \right],\ \ (7)\\ \left[ \begin{array}{c@{\quad}c@{\quad}c} 2A_{k,n}^{(1)} & 2A_{k,n}^{(2)} & 2A_{k,n}^{(3)} \\ 1 & 1 & 1 \\ 2B_{k,n}^{(1)} & 2B_{k,n}^{(2)} & 2B_{k,n}^{(3)} \end{array} \right] \left[ \begin{array}{c} \alpha_{k,n}^{(1)}\\ \alpha_{k,n}^{(2)}\\ \alpha_{k,n}^{(3)} \end{array} \right]&=& \left[ \begin{array}{c} \chi_{k,n+1}^{\phi}+\chi_{k,n+1}^u\\ \chi_{k,n}^v\\ \chi_{k,n-1}^{\phi}-\chi_{k,n-1}^u \end{array} \right],\ {\rm for}ñ\!\geq\! 1 . \ \ (8)\nonumber\\ \end{eqnarray} By solving Eqs. (6) to (8) for the vector containing the α_k,n^(j) coefficient(s) from the projections of (u',v',φ') on the PCFs, we obtain the coefficient of Kelvin waves from Eq. (6), the coefficients of the pair of eastward- and westward-propagating mixed Rossby-gravity waves from Eq. (7), and the coefficients of all inertia-gravity and Rossby wave triplets from Eq. (8). Specifically, the application of EWEIF to an instantaneous flow field is done with the following steps: (i) project u',v', and φ' onto the PCFs and perform a zonal FFT on those coefficients to retrieve χ_k,n^u,χ_k,n^v, and χ_k,n^φ for all k and n; (ii) calculate A_k,n^(j) and B_k,n^(j) for each j, n, and k; (iii) plug the results of (i) and (ii) into Eqs. (6) to (8) and then invert Eqs. (6) to (8) for each k and any desired n to solve the vector of α_k,n^(j) coefficient(s).

The system of Eqs. (6) to (8) describes a complete set of equations for expanding an instantaneous flow of a shallow-water equation model into its constituent equatorial wave class fields with a presumed equivalent depth uniquely. For this reason, we refer to it as the method of equatorial wave expansion of instantaneous flows (EWEIF). The coefficient α_k,n^(j) of each Matsuno mode of (n,j) can be calculated uniquely for a given instantaneous horizontal wave field, meaning that a full equatorial wave class field may be reconstructed using \begin{eqnarray} u'^{(j)}_n&=&\sum_{-K}^Ka_{k,n}^{(j)}U_{k,n}^{(j)} ,\nonumber\\ v'^{(j)}_n&=&\sum_{-K}^Ka_{k,n}^{(j)}V_{k,n}^{(j)} , \ \ (9)\\ \phi'^{(j)}_n&=&\sum_{-K}^Ka_{k,n}^{(j)}\Phi_{k,n}^{(j)}\nonumber \end{eqnarray} for each pair of (n,j), with the reminder that we only have j=1 for n=-1 (in which V_k,-1⁽¹⁾=0); j=1 and 2 for n=0; and j=1,2, and 3 for n≥ 1. The summation over all n and j will enable us to recover the original instantaneous horizontal wave field. In summary, one may follow the procedures below to obtain individual equatorial waves from instantaneous wave fields of (u',v',φ'):

(1) Obtain (χ_k,n^u,χ_k,n^v,χ_k,n^φ), by projecting (u',v',φ') on the PCFs in y and individual zonal waves in x according to Eq. (1). This yields all projection coefficients on the right-hand side of Eqs. (6) to (8);

(2) Invert the linear equations, Eqs. (6) to (8), to obtain α_k,n^(j), which are "unknowns" in these equations;

(3) Use Eq. (9) to obtain the fields (or maps) of individual waves (a single k, j, and n), or individual wave classes (e.g., a single j and n by summing up over k) or the total field (summing up over all k, j, and n).

Note that the temporal information of Matsuno modes or wave classes in the data is contained in the projection coefficients of (u',v',φ') on the PCFs because instantaneous fields at different times have different values of (χ_k,n^u,χ_k,n^v,χ_k,n^φ). By employing Eqs. (6) to (8) at each instantaneous time, we will obtain the temporal evolutions of all individual wave classes of equatorial waves from the data.

In the next three sections, we will validate and illustrate EWEIF for three scenarios. For the first scenario (section 3), we reconstruct a series of temporally evolving wave fields based on the exact (linear) solution of the equatorial β-plane shallow-water equation model, and then apply EWEIF to isolate the individual waves contained in the constructed wave fields. For the second scenario (section 4), we use a numerical model of a linear equatorial β-channel shallow-water equation model, starting at an arbitrarily specified local perturbation field, to generate temporally evolving wave fields. Then, we apply EWEIF to identify all individual waves from the output fields of the numerical model. For the final scenario (section 5), we briefly detail and present results of the application of the EWEIF method to the ERA-Interim dataset, in which we illustrate its utility in analyzing real data.

The main utility of EWEIF is to identify individual wave classes in a given instantaneous wave field. Therefore, the EWEIF method by itself does not directly provide the time scale information of individual waves. The temporal evolution of an individual wave can be obtained by applying the EWEIF analysis to a series of consecutive observations in a time period. In observations, there are no so-called "free waves" that propagate following the analytical dispersion relation derived from the linear model due to nonlinear interactions. Even if such nonlinearity is weak, one should not expect the temporal scale of observed waves to follow their dispersion curves exactly, because the same waves would get excited constantly and continuously throughout the tropics by diabatic heating events. We construct a synthetic dataset that contains waves freely propagating according to the dispersion relation for validating the EWEIF method's functionality in capturing the time scale of individual waves.

For this purpose, the synthetic wave field is generated as follows: Each wave class field is constructed individually using the spatial structure given in Eq. (3), with an assumed propagation speed as given by Eq. (2); namely, \begin{eqnarray} u'^{(j)}_n(x,y,t)&=&\sum_{k=-K \atop k\neq 0}^K a_{k,n}^{(j)}U_{k,n}^{(j)}(x,y)e^{-i\omega_{k,n}^{(j)}t} ,\nonumber\\ v'^{(j)}_n(x,y,t)&=&\sum_{k=-K \atop k\neq 0}^K a_{k,n}^{(j)}V_{k,n}^{(j)}(x,y)e^{-i\omega_{k,n}^{(j)}t} ,\ \ (10)\\ \phi'^{(j)}_n(x,y,t)&=&\sum_{k=-K \atop k\neq 0}^K a_{k,n}^{(j)}\Phi_{k,n}^{(j)}(x,y)e^{-i\omega_{k,n}^{(j)}t} ,\nonumber \end{eqnarray} where a_k,n^(j) is a complex constant for a given (n,j,k) and ω_k,n^(j) are the physical roots of the dispersion relation, Eq. (2). We wish to reiterate here that, for n=-1, there is only one root ω_k,-1⁽¹⁾=k with V_k,-1⁽¹⁾=0; there are only two physical roots j=1 and 2 for n=0; and j=1,2,3 for all n=1,2,… N.

Following the procedures outlined in section 2, we obtain the coefficients α_k,n^(j) for each of the 3 × N Matsuno modes (n,j) of zonal wave k at time t by inverting the linear algebraic equations given in Eqs. (6) to (8). Applying Eq. (9) enables us to obtain the wave fields of (u'^(j)_n,v'^(j)_n,φ'^(j)_n) from coefficients α_k,n^(j) for each of the 3\(\bot\) N Matsuno modes (n,j). The summation of (u'^(j)_n,v'^(j)_n,φ'^(j)_n) of these 3\(\bot\) N Matsuno modes (n,j) corresponds to the EWEIF-reconstructed total wave field at time t. Validations of the EWEIF method can be accomplished in several ways, one of which is to directly compare the original total wave field with the total wave field reconstructed by the EWEIF method at each time. For ease of reference, we refer to it as "total-field validation". We can also validate the EWEIF method at the wave level by comparing the map series of (u'^(j)_n,v'^(j)_n,φ'^(j)_n), given in Eq. (10), and the map series of (u'^(j)_n,v'^(j)_n,φ'^(j)_n), obtained by EWEIF, for all n and j (referred to as "wave-field validation"), and at the wave level by comparing the time series of a_k,n^(j)e^-iωk,n(j)t with that of α_k,n^(j)(t) obtained for all n, j, and k (referred to as "wave-evolution validation"). Obviously, both wave-field and wave-evolution validations only work in synthetic data in which we already know individual wave classes as well as individual waves before applying the EWEIF. The total-field validation can be used in observed data or model outputs.

As reviewed in the introduction, diagnostics of equatorial waves are typically based on techniques that require some temporal information of wave fields, such as spectral analysis techniques or wave decomposition of temporally filtered wave fields. The total-field validation is not applicable since the reconstructed wave fields are not directly from instantaneous wave fields. As a result, the validation is done by comparing the wavenumber-frequency spectra of the diagnosed waves against the dispersion curves of equatorial waves (e.g., Wheeler and Kiladis, 1999). Following (Wheeler and Kiladis, 1999), we perform a simple 2D FFT in zonal space and time, and then sum over all latitudes to obtain wavenumber-frequency spectra of the diagnosed waves. For the synthetic data, we can carry out the spectral validation using the total wave field or using individual wave class (or Matsuno mode) fields. If EWEIF successfully recovers the constructed equatorial wave fields, Eq. (10), the power spectra of the EWEIF-derived waves plotted in the wavenumber-frequency domain would be along the dispersion curves of these waves.

As a demonstration, we only use the first three PCFs (i.e., N=2) to construct the synthetic wave field. Therefore, we may include up to six Matsuno modes in the synthetic wave field: Kelvin (n=-1), EMRG (n=0), WMRG (n=0), EIG1 (n=1), WIG1 (n=1), and R1 (n=1). We set K=20 and T=1000. Without loss of generality, we first consider the complex amplitude a_k,n^(j)=(1,0) for all wave classes at each k, meaning that all waves in the synthetic data have equal amplitude. We carry out the total-field, wave-field, and wave-evolution validations. The nearly indistinguishable patterns between the fields obtained from EWEIF and the original synthetic wave fields at each time t already validate EWEIF's accuracy directly. Nevertheless, we still wish to report the validation of EWEIF's functionality in capturing the time scales of individual waves to connect EWEIF with the existing methods for diagnosing equatorial waves.

We wish to note that the results presented in this section may perhaps be considered trivial, given that the synthetic field is constructed in such a way that the EWEIF method would easily "invert" the field to get the wave class amplitudes. However, we must test that the EWEIF analysis is capable of recovering the dispersion curves of equatorial waves in even the simplest setting. Once we have established the accuracy of EWEIF for the simplest case, the later tests allow us to adjust the field in such a way that we may gauge the sensitivity of EWEIF to various features of a given input field, specifically regarding the choice of equivalent depth.

Shown by the shading in Fig. 1 is the wavenumber-frequency spectrum of the total wave field reconstructed using EWEIF on the synthetic data. Since the reconstructed total wave field is nearly indistinguishable from the original synthetic wave field, the wavenumber-frequency spectrum of the synthetic wave field overlaps (nearly) exactly with the shading shown in Fig. 1. It can be seen that the peaks of the power spectra of the EWEIF-derived waves are indeed along the dispersion curves of these waves (black lines). This validates the robustness of EWEIF in capturing the time scales of individual waves. However, there are two seemingly "major" discrepancies. The first is that the peaks of the wavenumber-frequency spectra are somewhat blurry and are not concentrated at their wavenumber and frequency as specified in the synthetic data, which are supposed to have pulse-like power at the exact wavenumber k and exact frequency ω_k,n^(j) with amplitude exactly equaling |a_k,n^(j)|². The second is that, despite having \(a_k,n^(j)\equiv (1,0)\) for all n, j, and k, the peak spectral power varies noticeably with n, j, and k. Because the same discrepancies also appear in the wavenumber-frequency spectra of the original synthetic wave field, they should not be regarded as evidence of inaccuracy of the EWEIF method. The first type of discrepancy is largely due to the sampling issue in discretizing the wave field in the x and t domains in producing the original synthetic data. In other words, if we discretize the wave field in a higher resolution in x and t in Eq. (10), the peaks of the wavenumber-frequency spectrum would become less blurry and more concentrated at their designated wavenumber and frequency, and the peak values would be closer to |a_k,n^(j)|². The sampling issue also causes spurious cancellations among waves belonging to different wave classes (or Matsuno modes), which leads to the second type of discrepancy. To confirm that the second type of discrepancy results from the spurious cancellations among waves belonging to different wave classes, we perform the wavenumber-frequency spectral analysis using the wave fields of the six wave classes individually. As indicated in Fig. 2, the second type of discrepancy has been alleviated greatly when performing the wavenumber-frequency spectral analysis using individual wave class fields without increasing the sampling resolution.

Figure 1.  Wavenumber-frequency spectral power of the total synthetic wave field (shading). The axes are non-dimensional wavenumber (abscissa) and frequency (ordinate). Dispersion curves (black lines) are drawn to show the signal corresponding to each wave class.

Figure 2.  Wavenumber-frequency spectral power of the geopotential height field (shading) of the individual wave classes as diagnosed by applying the EWEIF method to the total synthetic field. The axes are non-dimensional wavenumber (abscissa) and frequency (ordinate). Dispersion curves (black lines) are drawn to show the signal corresponding to each wave class.

Figure 3.  As in Fig. 2, except for the synthetic wave field with the Kelvin and R1 waves removed.

We also consider a non-uniform setting of a_k,n^(j) in constructing the synthetic wave field, to test if the EWEIF method would not produce spurious results for wave classes when their amplitudes are smaller than others or even do not exist in a given wave field. As an example of such tests, we remove all Kelvin and R1 waves in the synthetic field by setting α_k,-1⁽¹⁾=α_k,1⁽³⁾=0. As indicated in Fig. 3, there is no spectral power in the Kelvin and R1 fields as intended. The spectra of the remaining four wave classes match their own dispersion curves closely, as in the case shown in Fig. 2, and no anomalous spectral power exists as a result of the absence of some wave classes.

Figure 4 is the result of applying EWEIF to a field composed entirely of noise. In this way, we can confirm again that the methodology does not produce spurious results for wave class spectra when waves do not exist. As discussed before, a solution will be found for any field, and we must then verify that for a field composed of random noise, the resulting wave class spectra are also random noise. In doing so, we show that EWEIF does not introduce or preferentially partition any spectral power according to the dispersion curves, and that any spectra produced by the expansion comes strictly from the input data. The wavenumber-frequency spectrum of the results of the EWEIF analysis on a randomly-generated (or white noise) field in space and time (Fig. 4) shows that EWEIF does not produce spurious wave class fields that are not presented in the data.

Figure 4.  As in Fig. 2, except for a field with a white noise spectrum in both space and time without presence of any waves.

In reality, it is unlikely that h_e is known a priori and, in addition, multiple wave classes with differing depth scales could exist simultaneously. Here, we present further examination of EWEIF as it pertains to the choosing of h_e, and how this choice (correct or otherwise) may affect the result. To illustrate the nature of adverse effects due to incorrect equivalent depth, we repeat the otherwise identical EWEIF analysis on the same synthetic field in producing Fig. 2, except that the equivalent depth has been intentionally chosen incorrectly as one-half of the correct equivalent depth. In this case, the set of PCFs that the dataset is projected upon will have a narrower trapping latitude scale than that used in constructing the wave classes in the first place, leading to incorrect meridional basis structures of the waves diagnosed by the EWEIF method. The spectra in Fig. 5 show that some wave class fields have erroneous signals that are not along that wave's dispersion curves (black lines). For example, dispersion curves for the EIG1 and WIG1 waves show up erroneously as Kelvin and R1 waves (red lines), respectively. This indicates that power from the IG1 wave classes has leaked into the lower modes due to our choosing smaller equivalent depth. Note that no signal can leak from a symmetric mode into an anti-symmetric mode, and vice versa, because anti-symmetric and symmetric modes are orthogonal regardless of depth scale agreement. This explains that there is no signal leakage into (from) the mixed Rossby-gravity modes. Figure 5 further reiterates that EWEIF does not remove the signal or introduce a spurious signal simply by choosing an incorrect equivalent depth. It partitions all available spectral power within a single meridional mode projection to the wave classes and does so without altering the dispersion characteristics of the signals.

Lastly, we apply EWEIF to the same synthetic field as that used in Fig. 2, except for this test we add a WIG1 wave with an equivalent depth that is one-half of that for the rest of the waves. Figure 6 shows the individual wave class spectra diagnosed in this test, with the red curves indicating the dispersion relation for the new WIG1 wave. The EWEIF diagnosis was performed using the "correct" equivalent depth; that is, the equivalent depth used to create the original six wave classes for the test shown in Fig. 2. As in Fig. 2, EWEIF successfully separates the wave class fields with the chosen equivalent depth, and their spectra match the dispersion curves closely. We find that the additional signal due to the added WIG1 wave is present in both the WIG1 and R1 spectra. As in Fig. 5, this indicates the presence of a wave whose equivalent depth is not equal to the one chosen for the EWEIF analysis, which is an intentional result based on the design of this test. Because for a given equivalent depth there is no case of dispersion curves overlapping, we can safely assume that any signal that is present in the wave class spectra of multiple waves is not due to any wave of the chosen depth scale. Further analysis would be required on a case-by-case basis to determine the origin of such a signal. For example, in the case of Fig. 6, one could easily identify the signal as a WIG1 wave by observing that the red dispersion curve fits the anomalous signal closely.

Figure 5.  As in Fig. 2, except for the wave class fields diagnosed by assuming an "incorrect" equivalent depth for the PCF structures. The spectra still reveal the signal corresponding to the original, correct equivalent depth (black dispersion curves), even though some of the signal is misattributed. Red dispersion curves are drawn for the EIG1 and WIG1 wave signal leakage into the Kelvin and R1 wave spectra, respectively.

Figure 6.  As in Fig. 2, except for a field containing an additional WIG1 wave, which has a different equivalent depth than that used in Fig. 2. A red dispersion curve is drawn for the new WIG1 wave where its signal leaks into the WIG1 and R1 wave spectra.

In the previous section, the input data were designed so that only spatial and temporal structures adhering to the theoretical free solutions of the linear equatorial shallow-water system are allowed. It has been demonstrated that EWEIF is capable of recovering the original, fabricated equatorial waves (Fig. 2). In this section, we apply the EWEIF method to outputs of an equatorial β-channel shallow-water model to demonstrate that EWEIF is capable of recovering equatorial wave fields without knowing a priori exactly which waves exist in the data or their spatial and temporal scales.

We use a standard linear shallow-water model centered on the equator in a tropical latitude band on a β-plane with a periodic boundary in the zonal direction and rigid boundaries at 30^°S and 30^°N. The model employs standard values of a (Earth's radius), Ω (Earth's rotation rate), and g (the gravity parameter). Here, we illustrate EWEIF by considering the solution of an initial value problem obtained with h_e=0.3 m. There is no special purpose to using such a shallow equivalent depth for the illustration, except that it would favor the dominance of Rossby and Kelvin waves over inertia-gravity waves with a narrower meridional trapping latitude as well as slow propagation speeds of Kelvin waves and inertia-gravity waves. The specific solution used for the illustration below is generated by the model initialized at rest with a 20^°-wide hole centered at the equator and at 180^° longitude (Fig. 7a). Since it is a linear model, the depth of the initial hole (as well as its polarity) is inconsequential, which is set arbitrarily to 1 m. Obviously, there will be no mixed Rossby-gravity waves coming out of the initial hole due to its symmetry about the equator.

Figures 7b and c are height fields of the Kelvin and n=1 Rossby waves diagnosed by EWEIF for the initial time. Accompanied with Kelvin and Rossby waves are small-amplitude inertia-gravity waves (not shown). The EWEIF analysis reveals neither mixed Rossby-gravity waves nor any higher odd-mode waves, which have zero amplitude at the initial time as well as any time after, again confirming that EWEIF does not produce waves spuriously if they are not present in the dataset, as we demonstrated with synthetic wave fields in the previous section. The sum of all individual wave fields isolated by the EWEIF method is convergent to the total field displayed in the Fig. 7a.

Figure 7.  Maps of the geopotential height field (shading) of the initial condition (t=0): (a) total field; (b) Kelvin wave; (c) Rossby wave. The results in (b) and (c) are from the EWEIF method.

From the EWEIF analysis, we can clearly see that an isolated, motionless hole centered at the equator is made of Rossby waves of low pressure, Kelvin waves of low pressure, and eastward- and westward-propagating inertia-gravity waves that are all co-located at the center of the hole. Along the equator and in its immediate vicinity, the westerly wind associated with Rossby waves of low pressure is canceled out by the easterly wind associated with the Kelvin wave of low pressure. Away from the equator on the south and north rims of the hole, the zonal winds associated with Rossby waves of low pressure are canceled out by their counterparts associated with inertia-gravity waves. On the west and east rims of the hole, the meridional winds that circulate low pressure counter-clockwise in the Northern Hemisphere and clockwise in the Southern Hemisphere are canceled out by meridional winds associated with inertia-gravity waves. These cancellations in the velocity and accumulations of negative mass anomalies give rise to the isolated motionless hole. In other words, each wave is still moving underneath the hole, even at the initial time, except that the motions associated with these individual perturbations are canceled out, creating an apparent "calm" hole at that moment. As time goes by, these individual perturbations, together with their own mass anomalies and velocity fields, will propagate out of their initial positions to their own destinations dictated by the imbalance of physical forces associated with their motions and gradient of mass distribution. In this sense, EWEIF allows us to "predict" the portion of a (total) perturbation flow at a given instance that would propagate westwards as Rossby waves and the portion that would propagate eastwards in form of Kelvin waves, as well as the portion that would propagate out from the original locations as mixed Rossby-gravity and inertia-gravity waves.

The results shown in Fig. 8 confirm the accuracy of the EWEIF analysis made to the initial state in terms of the underlying dynamics. It is apparent that an exactly same shaped perturbation with the same amplitude identified as the Kelvin waves in Fig. 7b shows up in the longitude around 45^°W 100 days later. The EWEIF analysis performed on the field shown in Fig. 7a identifies the perturbation around 45^°W as a Kelvin wave (Fig. 8b), which is otherwise identical to that in Fig. 7b, except its location is shifted east by about 135^° in longitude 100 days later. In other words, one can predict the flow pattern shown in Fig. 8b from Fig. 7b using the phase speed of Kelvin waves (about 1.72 m s^-1 for h_e=0.3 m). The Rossby flow pattern shown in Fig. 7c appears at about 130^°E 100 days later in Fig. 8c, which has been confirmed as an R1 wave by the EWEIF analysis made to the total flow shown in Fig. 8a. Due to the dispersive nature of Rossby waves, the pair of low-pressure centers that are symmetric about the equator becomes zonally wider with smaller amplitude (Fig. 8 versus Fig. 7). Therefore, the flow shown Fig. 8c (together with its velocity fields) can be "predicted" from the flow pattern shown in Fig. 7c with the information of their ω_k,1⁽³⁾ used in the EWEIF analysis. Obviously, the "predictive" ability of the EWEIF analysis is not surprising at all, because it is built on the same dynamics framework as the numerical model, but it nevertheless validates the robustness of the EWEIF method in terms of physics, which is beyond the mathematical validation showing that the sum of the individual wave flows is convergent to the total wave flow.

Figure 8.  As in Fig. 7, except for the output of the numerical model at t=100 days.

Figure 9.  Hovmöller diagrams of geopotential height fields (shading) along the equator derived directly from the output of the numerical model (c), and the portions associated with (a) WIG1, (b) R1, (d) Kelvin, and (e) EIG1 wave classes, as diagnosed by the EWEIF method for the first 300 days of the model run.

Shown in Fig. 9 are Hovmöller diagrams of height fields along the equator derived directly from the output of the numerical model and the portions associated with WIG1, R1, Kelvin, and EIG1 wave classes, as diagnosed by the EWEIF method. It allows us to illustrate how the equatorial wave class fields diagnosed by the EWEIF method may be concatenated into time series to create the illusion of propagation, when only instantaneous wave fields are diagnosed. Besides cleanly isolating eastward-propagating Kelvin waves (Fig. 9d) and westward-propagating Rossby waves (Fig. 9b), the EWEIF method allows us to reveal in detail the spatiotemporal evolutions of inertia-gravity waves generated from the initial motionless hole. Specifically, the slowly eastward-drifting wide-stripe signals shown in Fig. 9e correspond to eastward-propagating inertia-gravity waves with the eastward drifting indicative of their eastward group velocity. As discussed in (Cane and Sarachik, 1976), westward-propagating equatorial long inertia-gravity waves can have both eastward and westward group velocity. This explains why the wide-stripe westward-propagating signals shown in Fig. 9a are much wider than their counterparts in Fig. 9e, because they correspond to a mix of the westward-propagating inertia-gravity waves with eastward group velocity (evident as the slow eastward drifting of the east edge of the stripe signals) and those with westward group velocity (evident as the slow westward drifting of the west edge of the stripe signals).

In addition to the wide-stripe signals, there exist (zonally) very narrow and westward drifting (temporally) stripe signals in the Hovmöller diagrams of both WIG1 (Fig. 9a) and EIG1 (Fig. 9e) waves. Their westward drifting speed is the same as the group velocity of R1 waves (note that the westward drifting speed shown in the Hovmöller diagram of R1 waves corresponds to group velocity of the R1 wave package, as it consists of many R1 waves with different zonal wavenumbers), which is faster than the westward group velocity of inertia-gravity waves (as indicated by the slower westward drifting of the west edge of the wide-stripe signals in the leftmost panel). This illustrates the ability of the EWEIF analysis to capture the continuous spinoffs of both westward and eastward inertia-gravity waves from Rossby waves as they propagate westwards dispersively. As elicited in the two-part series of papers by Cai and Huang (2013a, 2013b), the westward propagation of Rossby waves always involves crossing isobaric unbalanced flow. When the geostrophically balanced flow crosses the β-barrier (northerly and southerly winds), it excites crossing-isobaric unbalanced flows. The crossing-isobaric flow carries the air mass into/away from the original low-/high-pressure centers that are "enclosed" by geostrophically balanced paralleling-isobaric flow. Part of the air mass associated with "filling/leaking" processes will be surrounded by paralleling-isobaric flow again due to the Coriolis deflection of the crossing-isobaric unbalanced flows. This, together with the unfilled/unleaded part of the original low-/high-pressure centers, becomes the new geostrophically balanced low-/high-pressure centers in the next time step, which is responsible for changes in shape and amplitude of the original low-/high-pressure centers as they propagate westwards or the dispersive nature of Rossby waves (Cai and Huang, 2013b). The remaining part of the air mass associated with "filling/leaking" processes that has not been "enclosed" by geostrophically balanced paralleling-isobaric flow corresponds to the spin-offs of inertia-gravity waves. The discussion above explains why the spin-off signal itself propagates westwards with Rossby wave packages, whereas the signals within each spin-off both propagate westwards (Fig. 9a) and eastwards (Fig. 9e) at a speed even faster than Kelvin waves. It should be pointed out that such spin-offs of inertia-gravity waves from Rossby waves are neglected in quasi-geostrophic dynamics. According to (Cai and Huang, 2013a), the amplitude of the spin-off part is determined by the ratio of gravity wave frequency to inertia-gravity wave frequency, suggesting the spin-off is strongest near the equator.

In the previous sections, we detailed the ability of EWEIF to isolate shallow-water wave structures from both idealized and model data. We further wish to validate the methodology by applying it to reanalysis data in the equatorial stratosphere. Here, we wish to validate the EWEIF diagnosis of the evolution of stratospheric equatorial waves against their expected evolution as detailed in classical QBO theory (Lindzen and Holton, 1968; Holton and Lindzen, 1972; Baldwin et al., 2001). Specifically, we aim to examine if the stratospheric equatorial waves identified by the EWEIF method from instantaneous flow fields of the ERA-Interim dataset (Dee et al., 2011) would be composed of mainly eastward-propagating (westward-propagating) waves during the easterly (westerly) phase of the QBO.

We apply the EWEIF analysis to geopotential height, zonal wind, and meridional wind fields at the 10 hPa pressure surface derived from the four-times-daily ERA-Interim dataset covering an equatorial latitude band from 30^°S to 30^°N, with a 1.5^° longitudinal grid spacing, for the period 1 January 1979 through 31 December 2012. Following (Wheeler and Kiladis, 1999), we perform the spectral analysis using a sliding Fourier analysis technique with 120-day windowed segments and a 60-day overlap between consecutive segments. The mean and linear trend are removed from each segment individually. We also transform the data into a flow-relative coordinate so that the propagation characteristics revealed by the spectral analysis more closely represent intrinsic values with respect to a resting basic state [note that, in reality, the effect of a non-zero basic wind state extends beyond simple Doppler shifting, as discussed in (Yang et al., 2003), (Ern et al., 2008), and (Dias and Kiladis, 2014)]. A zonal FFT is then applied for every time and latitude, followed by a temporal FFT for each wavenumber and latitude. The wavenumber-frequency spectral power is obtained by squaring the complex spatiotemporal coefficients, averaging over all segments, and summing over all latitudes. The results of the wavenumber-frequency spectral power indicate that the spectral power peak matches the Kelvin wave dispersion curve with an equivalent depth of h_e=200 m (not shown here). We note that this value is nearly the same as that of the first baroclinic mode of the troposphere (Garcia and Salby, 1987; Salby and Garcia, 1987), which has been shown to project strongly upon dynamical variables (Wheeler and Kiladis, 1999). For this reason, we choose 200 m as the value of the equivalent depth, the only free parameter of the EWEIF analysis, to diagnose stratospheric equatorial waves in the ERA-Interim dataset.

Figure 10.  Hovmöller diagrams of the (a, c) Kelvin wave and (b, d) R1 wave geopotential height fields (shading) for the first 100 days of 1979 [(a, b) QBO locally easterly] and 1980 [(c, d) QBO locally westerly]. The diagrams have the same shading scheme to illustrate the change in amplitude across opposite QBO phases.

The EWEIF analysis yields 34 years (1979-2012) of four-times-daily maps of geopotential height, zonal wind, and meridional wind fields at the 10 hPa pressure surface for each of the following 12 equatorial wave classes: Kelvin waves, eastward- and westward-propagating mixed Rossby-gravity waves (EMRG and WMRG), eastward- and westward-propagating inertial-gravity waves with n=1,2, and 3 (EIG1-3 and WIG1-3), and Rossby waves with n=1,2 and 3 (R1-3). Following the same procedures outlined above, we perform the wavenumber-frequency spectral analysis for each of these 12 wave fields. We then compare the wavenumber-frequency spectra of individual wave classes with their dispersion curves obtained using an equivalent depth of h_e=200 m. Because the wavenumber-frequency spectral analysis only provides the time-mean properties of the dominant waves in each regime of the QBO cycle, we do not present the detailed results here of the wavenumber-frequency spectral analysis of these individual wave classes for illustrating the feasibility of applications of the EWEIF analysis in observations. Suffice it to say that we have confirmed (not shown here) that (i) the sum of the spectra of individual wave classes is very close to the raw spectra of the total geopotential height wave field, and (ii) the wavenumber-frequency spectra of individual wave classes capture the main characteristics of their corresponding dispersion curves, as has been shown in many other studies (e.g., Tindall et al., 2006a, 2006b; Alexander et al., 2008; Ern et al., 2008; Yang et al., 2011).

In the remaining part of this section, we turn our attention to illustrating how the EWEIF analysis might allow us to examine the temporally continuous evolution of individual equatorial waves in the stratosphere as the background QBO evolves with time. To the authors' knowledge, the continuous mutual evolution of individual equatorial waves and the QBO has not been examined from instantaneous fields, although the continuous evolutions of the wave-drag that drives the QBO (Ern and Preusse, 2009; Ern et al., 2014) and wave kinematics properties of the total field (Cai et al., 2014) have been analyzed. Specifically, we wish to validate the EWEIF method by illustrating that the continuous evolutions of individual wave fields diagnosed by EWEIF are consistent with the wave-mean flow interaction theory for the QBO (Lindzen and Holton, 1968), as has been well-documented (in terms of the time-mean properties of the dominant waves in each regime of the QBO cycle) in both observations (e.g., Tindall et al., 2006b; Alexander et al., 2008; Ern et al., 2008; Ern and Preusse, 2009; Yang et al., 2011) and simulations (e.g., Dunkerton, 1997; Giorgetta et al., 2006; Evan et al., 2012).

Figure 10 shows Hovmöller diagrams of the Kelvin and R1 waves over two otherwise arbitrarily selected periods, but one period for the easterly QBO phase and the other for the westerly phase. It can be seen that Rossby waves exhibit noticeably weaker amplitude in the easterly phase while Kelvin waves do so in the westerly phase, consistent with the wave selection mechanism via critical-level absorption of waves that propagate in the opposite direction of the mean zonal flow. Besides the evidence of wave amplitude modulations by the direction of the background zonal flow, we also see a lack of coherent continuous westward (eastward) propagation of Rossby (Kelvin) waves across the entire latitudinal circle during the easterly (westerly) phase of the QBO. This again is consistent with the critical-level absorption of waves that propagate in the opposite direction of the mean zonal flow.

Figure 11.  Equatorial wave class amplitude evolution synchronized with the background QBO cycle. On the abscissa, the QBO clock begins with a peak westerly phase (0), continues to a peak easterly phase (12), and returns to the peak westerly phase (24=0). The equatorial waves are grouped by propagation direction [(a) eastward; (b) westward]. Full domain geopotential height variance is used to represent wave class amplitude, with each curve normalized to its peak value.

We confirm that the amplitudes of wave classes diagnosed by EWEIF obey critical-level absorption effects throughout the entire time domain as the background QBO evolves. To succinctly summarize this feature, we present in Fig. 11 composite continuous evolutions of amplitude of individual wave classes with respect to the background QBO evolution. The amplitude of individual wave classes is defined as the total variance of the geopotential height field associated with individual wave classes diagnosed by the EWEIF in the entire latitudinal band between 30^°S and 30^°N. To show the evolutions of all wave classes in the same figure, we have normalized each curve by its own maximum value. Recalling that QBO theory predicts critical-level absorption, meaning that waves whose phase speed approaches the background mean zonal wind speed are dampened. Figure 11 vividly shows that those waves propagating to the east (west) have maximum amplitude in the easterly (westerly) phase and minimum amplitude in the westerly (easterly) phase. The ability to capture the temporally continuous evolutions of individual waves via EWEIF analysis allows us to uncover some new and detailed evidence for the underlying dynamic processes of the wave-mean flow interaction dynamics from observations. For example, it is apparent that the maximum and minimum amplitudes have a delay of 1-4 months from the peak phases of the QBO. Additionally, we find that wave classes are not entirely damped at the times when their propagation direction is along the direction of the background zonal wind. The multi-level analysis showing the vertical propagation of amplitude modulations for individual waves with different zonal wavenumbers (not shown here) suggests that the critical-level absorption acts as a wavenumber filter instead of a wave class filter, which dampens, individually and progressively, the waves that propagate along the zonal wind direction within their wave class field as they propagate towards their critical levels where \(\bar{u}=c(k)\). The results of the multi-level analysis of individual wavenumber selection mechanisms will be reported in a separate paper, as the main goal of this paper is to formulate and illustrate the EWEIF method.

This paper presents a formulation for a method of diagnosing individual equatorial waves from instantaneous total wave fields, derived either from observations or models, without prescribing a time or zonal length scale based on the equatorial dispersion relation in spectral space. By expanding flow fields into constituent equatorial waves using the free solutions to the shallow-water system as basis functions, we are able to diagnose instantaneous equatorial wave fields with spectral characteristics matching the dispersion relation. Of particular note is that the frequency spectrum can be diagnosed simply from time series of independently diagnosed instantaneous horizontal fields. In our analyses, we show that the frequency characteristics of equatorial waves are linked to their physical structure, and we demonstrate that EWEIF is able to illuminate this feature without consideration of the time dimension explicitly.

As with the spectral, meridional mode, and EOF decomposition methods for diagnosing equatorial waves, EWEIF relies on the shallow-water framework developed by (Matsuno, 1966) to identify the wave classes. In these methods, there is an a priori (or in the case of EOF analysis, a posteriori) assumption that the wave classes will adhere to some structure presumed by Matsuno's theory. The methods that use an assumption regarding the physical structure are beneficial to diagnosing the evolution of equatorial waves because the evolution is not predicated on a time filter, which inherently forces the evolution to reflect said filter. EWEIF further builds upon the meridional mode decomposition method (i.e., Yang et al., 2003) by separating signals between waves that share PCF projections.

The only free parameter in the EWEIF analysis is the equivalent depth for equatorial waves, which is also treated as a free parameter in all other methods that are based on equatorial wave dynamics. Choosing an equivalent depth is an integral part of the EWEIF method, because it controls the values of the A_k,n^(j) and B_k,n^(j) coefficients in Eq. (4), as well as the meridional trapping latitude through the frequency ω_k,n^(j). Analysis of results from testing EWEIF on synthetic data with a single vertical mode shows that the wave class fields are sensitive to the choice of equivalent depth only in that the signal may be attributed to the wrong wave classes if the wrong equivalent depth is chosen. Despite this, the signal maintains the correct structure in spectral space related to the true wave class field even when misattributed. Similarly, the presence of multiple equivalent depths in the dataset will be evident in the wave class spectra since the spectral power will match a dispersion curve other than one with the correct h_e. In both cases, the correct equivalent depth(s) may still be determined if the fields are not completely and correctly separated in spectral space by fitting the signal power to a dispersion curve. Any spectral power from waves with the assumed equivalent depth will be correctly diagnosed. It should be pointed out that, by design, the EWEIF method intrinsically supports the provision of using temporally evolving equivalent depth, but presently does not include an algorithm for the diagnosis or "best fit" of an equivalent depth for a given instantaneous flow field. Work is currently in progress to extend EWEIF's ability to also estimate the equivalent depth from instantaneous flow. In this case, we would be able to provide an in-line value of h_e, rather than choosing a best guess a priori using dynamical theory or the field spectra.

We also test the EWEIF method on ERA-Interim data and find that EWEIF successfully partitions the wave classes in spectral space from their physical fields. We validate the EWEIF method by confirming the mutual evolutions of equatorial waves with the background QBO in the stratosphere, as predicted by the wave-mean flow interaction theory for the QBO. The ability to capture the temporally continuous evolutions of individual waves via EWEIF analysis allows us to uncover some new and detailed evidence for the underlying dynamic processes of the wave-mean flow interaction dynamics from observations. For example, it is apparent that the maximum and minimum amplitudes have a delay of 1-4 months from the peak phases of the QBO. Additionally, we find that wave classes are not entirely damped at the times when their propagation direction is along the direction of the background zonal wind.

Because EWEIF only decomposes the instantaneous horizontal field of individual wave classes, we are able to produce time series of wave fields. In a similar fashion, we are also able to extend the domain of our study through multiple vertical levels. The fields at various levels would be diagnosed independently, so any features of the vertical structure or propagation characteristics of the wave classes would be diagnosed from information contained within the data, as opposed to the prescription of a modal structure or wavenumber spectrum. The goal of such a study would be to determine the primary vertical modal features and the evolution of equatorial waves in the vertical direction.

It is of importance to point out that the EWEIF method merely decomposes wave fields in data into individual equatorial waves based on the linearized shallow-water system with a motionless mean state that has a specific equivalent depth (or thermal structure). It is not expected or assumed that the actual evolutions of individual waves in data would follow linear equatorial waves. Nevertheless, the differences between the actual evolutions, as diagnosed from the instantaneous flows in the next time step (the next day, say) and those inferred from the linear dynamics of the individual waves diagnosed in the current time step, would be indicative of the combined effect of nonlinearity, basic flow, and diabatic heating, as well as the incorrect estimate of the equivalent depth used in the EWEIF analysis.

There is continued interest in tropospheric wave activity, as the tropical troposphere is the source of waves that drive the QBO, and is host to the MJO and the atmospheric component of ENSO. A challenge remains to study the interaction of waves with their environment within the noisy background of the lower atmosphere. Near genesis, spectral filters could excise meaningful signals representative of some transient state of the wave in its infancy. Being a full spectrum methodology, EWEIF could serve to illuminate the evolution of a wave from its point of genesis into the far field, as we have demonstrated in this paper. In particular, this method could be applied to a high-resolution model of a localized heating event to observe the temporal evolution of the wavenumber spectrum as the wave propagates away from its source, as well as the favored projection of some depth scale or selective initiation of some waves over others.

Use of EWEIF is not limited to theoretical models. The QBO is a feature that is important for long-range predictions but is under-represented in models and under-utilized in forecasting systems (Hamilton et al., 2015). Similarly, the MJO remains a challenge to accurately represent in GCMs (Deng et al., 2015). In both cases, insufficient parameterization of scale interactions is the likely culprit because convection and the resultant equatorial wave activity play a major role in the development and maintenance of these oscillations. To compare a GCM simulation of the QBO or MJO to reality, we may apply EWEIF to the model output and see how the individual waves are treated. A study of this fashion could advance parameterizations and lead to more accurate long-term forecasting for weather and climate.