2025 年 103 巻 3 号 p. 371-388
The earliest attempts to study the global normal mode oscillations of the atmosphere used time series of barometric in situ observations, but such approach is limited by the spatial and temporal inhomogeneity of meteorological station data. A major advance on the subject was recently made by applying a zonal-time spectral analysis to the surface pressure field in hourly gridded ERA5 reanalysis data, which disclosed an array of spectral peaks at theoretically predicted zonal wavenumber-frequency pairs, including many peaks with periods between 2 hours and 12 hours. However, this result relies on adequate representation of the modes in ERA5, which (i) ingests data sources that cannot explicitly resolve high frequency modes (e.g., radiosondes and polar satellite observations), and (ii) employs a numerical forward model that potentially introduces spurious effects. The present study provides “ground truth” for the reanalysis by a simple analysis of hourly barometric observations taken at ∼ 3800 stations over the globe. For each putative global mode, a time series of its index is computed by filtering the hourly ERA5 pressure fields. The station data is then regressed onto this index, revealing, for each mode, a characteristic, globally coherent spatial pattern of regression coefficients. The meridional structures of the regression patterns agree fairly well with the corresponding Hough functions, not only for low-frequency Rossby and Rossby-gravity modes, but also for high-frequency modes such as Kelvin and inertia-gravity modes. Even the Pekeris resonance is identified for a couple of Kelvin modes. These findings both solidify the evidence for a rich spectrum of global normal modes in the real atmosphere and also lend credence to their representation in ERA5. It is impressive that ERA5, by combining a numerical model with scattered meteorological observations, even reproduces the tiny (∼ 0.1 – 1 Pa amplitude) pressure signals of the high-frequency global normal modes.
This study aims at confirming the presence of the recently reported array of normal modes (Sakazaki and Hamilton 2020; hereinafter SH20) through a simple analysis of raw barometric observations taken at many individual stations over the globe.
In common with other natural systems, the global atmosphere displays normal mode (also called resonant or free) oscillations, occurring at discrete frequencies and each associated with its own horizontal and vertical structure. Resonant mode solutions are predicted by classical tidal theory that considers the inviscid primitive equations for an atmosphere above a rotating, smooth, spherical Earth. These equations are linearized about a motionless mean state with temperature assumed to be a function only of the vertical coordinate (e.g., Chapman and Lindzen 1970). Within the approximations of classical tidal theory, the governing equations can be combined into a single second-order partial differential equation for, e.g., the geopotential perturbation. That equation is separable into ordinary differential equations for the time, zonal, meridional, and vertical domains. The solutions in time and the zonal direction are simple Fourier harmonics, while the meridional equation is the Laplace tidal equation (LTE) with known Hough function solutions (Longuet-Higgins 1968; Kasahara 1976). The vertical structure equation (VSE) is a second-order equation that requires boundary conditions at the Earth’s surface and a specified “top of the atmosphere.” The normal mode oscillations correspond to the homogeneous solutions in this system (i.e., the eigen solutions in case of no external forcing) and the geopotential perturbation for each mode is expressed as
 |
where z is altitude, θ is latitude, k is zonal wavenumber, λ is longitude and t is time. Here, each vertical mode (m: vertical mode index) is characterized by the equivalent depth (hm) and the vertical structure function [Zm (z)] that are obtained as the eigenvalue and eigenfunction of the VSE, respectively. For each vertical solution [hm, Zm (z)], there is an associated complete set of horizontal modes (i.e., LTE solutions), describing linear shallow-water waves with a specific zonal wavenumber (k) component in a fluid with depth hm on a rotating sphere. The LTE solutions are associated with infinite sets of eigen-frequencies
and meridional structure functions [Hough function 
].
Because the atmosphere is unbounded at its upper limits, only one or two eigen solutions(s) are theoretically predicted for the VSE (i.e., m = 1, 2), at least for a realistic vertical mean temperature profile (Salby 1979; Ishioka 2023; Ishizaki et al. 2023). The more robust solution is the so-called “Lamb resonance” (m = 1, h1 ∼ 10 km), characterized by a Lamb wave structure in Z and energy trapped near the surface (Lamb 1911). The other solution is “Pekeris resonance” (m = 2, h2 ∼ 6.5 km) with its energy trapped both around the stratopause and the surface (Pekeris 1937). The Lamb resonance solution satisfies the physically reasonable “top of the atmosphere” boundary condition and is robust to changes in the assumed mean temperature profile. By contrast, the existence of Pekeris resonance depends on the detailed mean temperature profile assumed. Ishioka (2023; see his Fig. 2) showed that the Pekeris solution does not exactly satisfy the upper and lower boundary conditions and so may correspond to a mode that would actually leak some energy rather than act as a perfect resonance. The relevance of the Pekeris mode in real atmospheric flow has been doubted (e.g., Salby 1980) and clear observational evidence for the presence of this mode has only recently been obtained (see below).
Figure 1 shows the theoretical dispersion curves for Lamb (h = 10 km; closed circles) and Pekeris (h = 6.5 km; open circles) resonance. The modes are further classified into Rossby (blue), Rossby-gravity (orange), Kelvin (red), and inertia-gravity (magenta) modes (Matsuno 1966). Note that eastward components in Rossby-gravity modes are sometimes classified into inertia-gravity modes (“n = 0 eastward inertia-gravity modes,” e.g., Kiladis et al. 2009), but this study refers to both westward and eastward components as Rossby-gravity modes.
Most previous attempts to find observational evidence for various modes in the real atmosphere considered low-frequency modes, mostly Rossby or Rossby-gravity modes (see SH20 and references therein). Because of the “red” nature of the power spectrum, these modes have relatively large amplitudes and thus are easier to detect. In their pioneering work on the westward propagating 5-day wave, Madden and Julian (1972, 1973) found a westward propagating signal with a period of ∼ 5 days by performing a cross-spectral analysis of surface pressure data from world-wide stations. Along with the composite analysis, they showed that its meridional structure agreed closely with the Hough function of the corresponding normal mode. Later studies often used geopotential height data from global data assimilation products (i.e., analyses or reanalyses) and satellite measurement data to identify several Rossby and Rossby-gravity modes, such as those referred to as the 4-day wave, 10-day wave, and 16-day wave (e.g., Ahlquist 1982; Hirota and Hirooka 1984; Madden 2007; Sassi et al. 2012; Madden 2019; Sekido et al. 2024).
For high-frequency modes, Matsuno (1980) found a signal of k = 1 Kelvin mode using a cross-spectral analysis of barometric observations at a few tropical stations. This was later confirmed by extended analysis of surface pressure data by Hamilton (1984) and Matthews and Madden (2000), and it is now known as the “33-hr Kelvin wave.” Hamilton and Garcia (1986) applied spectrum analysis to an exceptionally long record of raw barometric data taken at Batavia (6°S). Several peaks were found in the high-frequency band (≲ 1 day), which the authors tentatively identified as theoretically predicted high-frequency modes, including Kelvin and inertia-gravity modes. A few inertiagravity modes have been also tentatively identified more recently by Shved et al. (2015) and Ermolenko et al. (2018).
A much ampler array of free modes was discovered by SH20 based on zonal wavenumber-frequency spectral analysis of near equatorial surface pressure data from ERA5 at hourly temporal resolution. The SH20 result for the ratio of spectral power to the background noise for equatorially symmetric and anti-symmetric components is revisited in Fig. 1 (color shading), using gridded surface pressure data between 20°S and 20°N during 1980 to 2021 (See SH20 for the definition of background noise and other details). It is apparent that there are numerous isolated spectral peaks whose frequencies match the expected (slightly Doppler-shifted) normal mode frequencies of Lamb resonance (h1 = 10 km; Fig. 1). SH20 found that the meridional and vertical structures also matched well those expected from classical theory (i.e., Hough function, Θ (θ), and vertical structure function, Z (z), respectively).
All noted modes were identified as Lamb resonances (h1 = 10 km; closed circles in Fig. 1). Recently, by analyzing the pressure pulse forced by the volcano eruption at Tonga Hap’pai in 2022 and by reexamining the SH20 spectral results, Watanabe et al. (2022) discovered the other type of normal mode resonance, namely, Pekeris resonance (h2 = 6.5 km). Indeed, one can see small spectral peaks of k = 1 and k = 2 Kelvin modes for the Pekeris resonance in Fig. 1 (red open circles; see also Fig. 3).
The results by SH20 and Watanabe et al. (2022) appear to have provided solid evidence of the atmosphere “ringing” regularly at many resonant zonal wavenumber-frequency pairs. One might wonder, however, how accurately ERA5 (or, more generally, any atmospheric reanalysis) can represent such free oscillations, given potential structural errors in the numerical forecast model or a shortage of pertinent high-frequency observations in the data assimilation. One aspect of interest in this regard is the treatment of the upper boundary in the dynamical model, which has to be somewhat arbitrary. Any model boundary condition that acts to reflect energy back downward from the top model level could potentially introduce spurious free oscillations (Lindzen et al. 1968; Kasahara and Shigehisa 1983). This concern might be most relevant for the large-scale, high-frequency Kelvin, Rossby-gravity, and inertia-gravity modes identified by SH20, as these should be least affected by various sources of mechanical damping in the model. Moreover, in regard to assimilation, two key data sources, namely twice-daily balloon-borne radiosondes and polar orbiting satellites, cannot by themselves adequately resolve waves with periods less than 12 hours, potentially leading to distortions of the high-frequency modes identified by SH20. Since free oscillations could be excited by various internal processes of the atmosphere such as cumulus convection and baroclinic activity (e.g., Miyoshi and Hirooka 1999; Zurita-Gotor and Held 2021), some of the normal mode signals in the ERA5 might be self-generated in the model, independent of the signals in the real atmosphere.
Whereas the accuracy of reanalyses is relatively clear in the case of low-frequency modes (e.g., Rossby and westward Rossby-gravity modes; see Sakazaki 2021), SH20 also made an argument for ERA5 to accurately depict high-frequency, global-scale atmospheric variability. In particular, SH20 noted that the principal lunar semidiurnal (L2) oscillation, as seen in sub-daily barometric observations at many individual stations throughout the world (e.g., Haurwitz and Cowley 1969), has been previously shown to be well captured in the predecessor of ERA5 (Kohyama and Wallace 2014; Schindelegger and Dobslaw 2016, hereafter referred to SD16). Since the gravitational forcing is not included in the underlying dynamical model, any analyzed L2 signal must be introduced through the assimilation of real data. This inference in turn implies that at least some portion of high-frequency normal modes (those with amplitudes of > 0.1 hPa and periods of < 12 hours) would be also realistically represented and constrained by observations.
Although encouraging, the SD16 “ground truth” study applies directly only to L2. In the present paper we report on a somewhat similar effort to obtain direct evidence of normal modes in raw measurements, specifically globally-distributed, hourly pressure observations. However, there are difficulties in extracting the normal mode signal exclusively from raw barometric data at individual stations, even when decades-long records are available (Hamilton and Garcia 1986). One problem is that the mode amplitudes are generally small compared to the background noise, while also the frequencies of some modes are too close together to be separated well enough (cf. Fig. 3 of Hamilton and Garcia 1986). In this regard, cross-spectral analysis or zonal wavenumber-frequency spectrum, as was done by SH20, would be desirable; however, such a procedure is challenging due to the inhomogeneity of raw barometric data, both in time and space (see Fig. 2).
Here, we pursue a different approach to obtain ground-based evidence for the wide variety of global normal modes. Specifically, we create a “pacemaker index” of normal mode signals from ERA5 surface pressures. By regressing raw barometric data onto this putative mode index time series, we can condense these data to any fluctuations synchronized with the normal modes in ERA5. A key advantage of our method is that, at individual stations, we can calculate the regression coefficient even for irregularly sampled observations. By examining the horizontal distribution of the regression coefficients, it is also possible to delineate the mode structure. Of course, this approach would work only if the normal mode variations in ERA5 are realistic, but it will be demonstrated that the extracted signals have a clear global wave structure that closely matches theoretical predictions.
The remainder of the manuscript is organized as follows. Section 2 describes the datasets, while Section 3 explains the analysis procedure, underpinned by step-by-step examples. Section 4 presents the meridional structure for various normal modes. Finally, Section 5 summarizes the main findings and discusses some implications.
Following SD16, we mainly analyze data from the International Surface Database (ISD, Smith 2011). This dataset contains surface meteorological observations taken at nearly 20,000 stations over more than 100 years (1900 to present). For the work at hand, we consider sea level pressure data over 42 years, between 1980 and 2021 (i.e., restricted to the period in which satellite data are routinely assimilated into atmospheric reanalyses). Sampling times and intervals vary by station, and sometimes even change at a given station over the observation period. As an approach to quality control, we focus on the lunar semidiurnal tide (L2). Compared to synoptic observations, L2 signals are quite small, about 5 Pa in amplitude, and hence their representation serves as an indication for the data quality at a given station. We therefore subset the ISD to ∼ 4,100 stations that were considered to feature realistic L2 signals in the analysis of SD16.
Zonal wavenumber-frequency spectrum for equatorially (a, b) symmetric and (c, d) anti-symmetric components calculated with ERA5 pressure data for 20°S–20°N during 1980–2021. The ratio of spectra to the background spectra are presented, with the color bar shown on the right. Panels (a) and (c) span the frequency range up to 12 cycles day−1 (cpd), while panels (b) and (d) do so up to 2.5 cpd. Solid circles and open circles [only for panels (b) and (d)] are the theoretical dispersion curves for h1 = 10 km and h2 = 6.5 km, respectively, with their colors representing the wave type: Blue, red, orange, and magenta are for Rossby, Kelvin, Rossby-gravity and inertiagravity modes, respectively (note that eastward components in Rossby-gravity modes are sometimes classified into inertia-gravity modes (“n = 0 eastward inertia-gravity modes”). For Rossby and inertia-gravity modes, the first three gravest modes are shown.
Distribution of (blue circles) ISD, (magenta triangles) buoy and (light blue circles) ISPD stations used for normal mode analysis.
Since the equatorial region is rather sparsely covered by the ISD stations, meteorological buoy data with a sampling of 10 minutes from the Global Tropical Moored Buoy Array (GTMBA) Program at 43 stations are also analyzed. The GTMBA program consists of three subsets: (i) the Tropical Atmosphere Ocean/ Triangle Trans-Ocean Buoy Network (TAO/TRITON) over the Pacific (McPhaden et al. 1998), (ii) the Prediction and Research Moore Array (PIRATA) over the Atlantic (Bourles et al. 2008), and (iii) the Research Moored Array for African–Asian–Australian Monsoon Analysis and Prediction (RAMA) over the Indian (McPhaden et al. 2009). Note that portions of these datasets were previously used to detect small L2 signals over the ocean (SD16; Sakazaki and Hamilton 2018).
We additionally consider data at 43 stations during 1980–2015 from the International Surface Pressure Data bank (ISPD v4) (Compo et al. 2019) that were again found to have a good representation of L2 signals by SD16.
For further quality control, we only analyze data from stations at altitudes < 1,000 m and with records including at least 10,000 individual observations (we confirmed that eliminating the criterion for the altitude does not significantly change the results). In total, 3,734 stations from ISD, 42 buoys and 21 ISPD stations are used; see Fig. 2 for their location across the globe. As mentioned above, the data distribution is far from uniform in either latitude or longitude: ISD/ISPD stations are densely concentrated in Europe, Eastern Asia, North America, and Australia, while there are few over the tropical regions of the African and South American continents or over the ocean. Tropical buoys are distributed somewhat uniformly in the zonal direction, filling gaps in the ISD/ISPD spatial coverage over the tropical ocean. In any case, the inhomogeneity in distribution should be borne in mind as a possible limitation of our analysis, especially for the Southern Hemisphere and for high-order normal modes characterized by shorter length scales. Hereafter, we refer to our final compilation as “ISD/Buoy” data. Note that most of these data are assimilated in ERA5 (Hersbach et al. 2020).
2.2 ERA5ERA5 (Hersbach et al. 2020), the latest atmospheric reanalysis by the European Centre for Medium-Range Weather Forecasts (ECMWF), is used for creating the normal mode index, as well as for examining the horizontal mode structure that will be compared to that deduced from ISD/Buoy data. Similar to other atmospheric reanalysis, ERA5 was constructed by adjusting the forward-integrated state of a numerical weather model to agree, within specified uncertainties, with available in situ and remote-sensing observations. We analyze hourly ERA5 diagnostics for surface pressure and mean sea-level pressure at a grid spacing of 1° over the period 1980–2021 (Hersbach et al. 2023).
The working hypothesis of our analysis is that ERA5 represents normal mode oscillations in a realistic manner. Using the ERA5 surface pressure data, we create single time series that depict the magnitude and phase of each mode. In short, individual time series of ISD/Buoy data are then regressed onto this ERA5-based index. If the global pattern of regression coefficients agrees with the theoretical mode structure, we can conclude that our working hypothesis is true and that the raw barometric (ISD/Buoy) data contain global normal-mode signals that are consistent with those evident in ERA5.
The regression analysis is repeated for two versions of ERA5 sea level pressures: (1) a subset of the ERA5 data, sampled at times and locations of the ISD/Buoy stations (this is hereafter referred to as “Obs-sampled-ERA5”), and (2) the entire gridded, hourly dataset over 1980–2021 (referred to as “All-ERA5”). By comparing the two regression branches, we can evaluate how the results are affected by the inhomogeneous global distribution of the ISD/Buoy data.
In the following subsections, we explain the analysis strategy in more detail by showing two examples: Kelvin modes with k = 1 and 5 (these are referred to as KL1 and KL5, respectively). Note that the former has been identified as the “33-hour Kelvin wave” in station barometric data (Matsuno 1980; Hamilton 1984), while the latter has only been identified so far through analysis of the ERA5 data (SH20).
Frequency spectrum for the equatorially symmetric, eastward-propagating zonal (black) wavenumber 1 (E1) and (gray) 2 (E2) components as deduced with surface pressure data between 20°S and 20°N in ERA5. The spectrum is calculated every year from 1980 to 2021 and the results over the 42 years are averaged. The E2 spectrum is multiplied by a factor of 10−1 for better presentation. The spectral peaks marked by “L” and “P” denote the Lamb (m = 1) and Pekeris (m = 2) resonance, respectively. Red and blue horizontal lines denote the frequency range used for the filtering for Lamb and Pekeris peaks, respectively: For the Lamb peaks, the range is determined objectively based on the fitting to the Lorentzian function (green solid curves; the fitting parameters are adopted from Table 1 of SH20), while it is determined subjectively for the Pekeris peaks; see the main text.
Following SH20 and Sakazaki (2021), we create the normal mode index (time series) from meridionally averaged ERA5 surface pressures, adopting the 20°S–20°N latitude band for the equatorially symmetric component, and the difference between 0°–20°N and 20°S–0° (divided by a factor of 2) for the equatorially anti-symmetric component (note that the results do not change significantly when sea level pressure data, instead of surface pressure data, are used for producing the index). Figure 3 exemplarily shows the frequency spectrum of the equatorially symmetric, eastwardpropagating k = 1 and 2 components. The spectral peaks denoted by “L” correspond to those for Kelvin modes of Lamb resonance and are well approximated by the Lorentzian function; the red horizontal lines are estimates of the width of each resonance based on an objective Lorentzian fit (see below). The filtering is based on the two-dimensional Fourier transform, mapping space to the zonal wavenumber-frequency domain: Only spectral coefficients with the corresponding zonal wavenumber k and within the frequency range (f0 − 3d, f0 + 3d) are retained (the other Fourier components are simply set to zero) and are Fourier transformed back to the physical space. Here f0 and d are the central frequency and spectral width, respectively, empirically determined by the fitting to the Lorentzian function (see Table 1 of SH20). Note that the Lorentzian function is defined such that it takes its maximum at f0 and decreases by a factor of 10−1 at f0 ± 3d.
For m = 2 vertical modes (Pekeris resonance), two Kelvin modes (k = 1 and 2) are examined, as in Watanabe et al. (2022). This limited selection is understandable from Fig. 1; the peaks for Pekeris modes (denoted by “P”) are well separated from those for Lamb modes (denoted by “L”) only for high-frequency modes (Kelvin or inertia-gravity modes). On the other hand, the peak amplitude is one order smaller than the Lamb series (m = 1) (cf. Fig. 3), making it rather difficult to identify high-frequency modes such as Kelvin modes with large k or inertia-gravity modes. Such a trade-off results in our study being confined to two isolated Pekeris resonance peaks. For the filtering, we set (f0, d) = (0.58, 0.02) for k = 1 and (f0, d) = (1.165, 0.038) for k = 2 (see blue horizontal lines denoted with “P” in Fig. 3; these are subjectively determined instead of using Lorentzian fitting).
The actual index time series is calculated every year, using the 365 (or 366 for leap years) ± 20 day data. The resultant pressure anomalies are a function of longitude and time. From these data, four time series at 0°E lagged by lπ/4 (l = 0, 1, 2, 3) are produced and used as normal mode indices after normalized by their standard deviation over time. Figures 4a and 4b illustrate the index time series for KL1 and KL5, respectively, over a week in early October 2010. It is evident that the index indeed oscillates with the characteristic periods (∼ 32 hours for KL1 and ∼ 7 hours for KL5; see Fig. 1), while amplitudes change over a much longer time scale.
Illustration of the procedure to extract normal mode signals from ISD/Buoy data for (left columns) KL1 and (right columns) KL5. (a, b) Normal mode index created with the filtered, normalized tropical surface pressure data from ERA5 at 0°E with different lags: (blue) 0 and (orange) π/2. (c–f) Regression coefficients on the normal index time series with lag = 0 as deduced from (c, d) ISD/Buoy data and (e, f) All-ERA5 data (unit: Pa). Color bars are shown at the bottom of panels (e, f). (g–j) Zonal distribution of regression coefficients (unit: Pa) obtained for the index with (g, h) lag = 0 and (i, j) lag = π/2, with the upper and bottom panels showing the results for 35°N (32.5–37.5°N) and 0°N (2.5°S–2.5°N), respectively. Blue, closed circles are for ISD and ISPD data, magenta, triangles are for buoy data, gray curves are for gridded ERA5 data, and red solid curves are the harmonic fitting for (left) k = 1 and (right) k = 5 for observation (ISD/Buoy/ISPD). See text for details.
We regress the time series of raw barometric data from ISD/Buoy onto the normalized normal mode index, after first removing the annual and semiannual harmonics. For the actual calculation, the index is linearly interpolated to the times of observation at individual stations. Figures 4c and 4d show the global distribution of the regression coefficient (R) on the normal index time series with l = 0 (lag zero) for KL1 and KL5 modes, respectively, as deduced from ISD/ Buoy data, and corresponding coefficients from All-ERA5 data are depicted in Figs. 4e and 4f. We see that the regression pattern of ERA5 (Figs. 4e, f) exhibits a Kelvin mode structure with the expected wavenumbers [i.e., k = 1 for (e), and k = 5 for (f)]. Similar inferences can be made for ISD/Buoy (Figs. 4c, d), although data coverage in the tropics is too sparse to clearly delineate every modal trough and high for KL5. We emphasize that the results in Fig. 4 are obtained without any a priori assumption for the modes’ horizontal structure.
Next, a zonal wavenumber (harmonic) fitting is performed for R to obtain the meridional structure in amplitude and phase of the target zonal wavenumber components [again, k = 1 (5) for KL1 (KL5)]. The fitting is done using data binned in latitude bands with 5° width. The calculation is made only when there are at least 20 valid data points in longitude, after removing values exceeding the 3-σ level in each latitude band (σ: standard deviation in zonal direction). Figures 4g – j, for instance, show how the fitting works for latitude bands of 2.5°S–2.5°N and 32.5–37.5°N, with the results based on two lagged indices (cf. Section 3.1, Figs. 4a, b). Again, the zonal distribution of R is clearly represented by a single harmonic, albeit modulated by noticeable short spatial scale variability in the extratropics.
A special treatment is necessary for the k = 0 component for which R (ideally) takes the same value at all longitudes, since it is the zonally symmetric component (in this case, zonal wavenumber fitting cannot determine the amplitude and phase). As shown in Fig. 5a for k = 0 Rossby-gravity wave mode, the zonally averaged R changes with the lag for the index (l). We thus determine its amplitude (A) and phase (δ) by harmonic fitting to the zonal-mean R (at each latitude belt) as a function of lag, i.e., A ‹ 
as shown by red curve in Fig. 5b.
(a) Meridional structure of the zonally averaged regression coefficient (R) (unit: Pa) calculated for normal mode index time series with the four different lags (x-axis) for Rossby-gravity wave mode k = 0, as derived with ISD/Buoy data. (b) Zonally averaged R for 40°N (37.5°–42.5°N band) for the four different lags (black circles). Red curve shows the harmonic fit, wherein the gray, dashed vertical line denotes the phase determined with this fit (it is nearly zero in this case).
The 95 % confidence intervals for the calculated amplitude and phase are estimated using a bootstrap method. We iterate the calculation of amplitude and phase values for the zonal harmonic of interest 1,000 times with resampled datasets, where each resample is generated from the original data with random replacement. The 97.5 and 2.5 percentile values are obtained from the resulting distribution and taken as the upper and lower confidence bounds (see error bars in, e.g., Fig. 6).
Meridional structure of (top) amplitude and (bottom) phase for Kelvin modes of Lamb resonance (m = 1) for (from left to right) k = 1 to 5. Red circles denote the results from ISD/Buoy, with their vertical bars showing the 95 % confidence level. The number in the parentheses at the top of each panel denotes the wave frequency (unit: CPD). Gray curves (for amplitude) and open circles (for phase) show the results for All-ERA5 data, while green curves for amplitude represent the Hough function fitted to the results for All-ERA5 (gray curves). The phase represents the lag from the index time series and is drawn only if the amplitude is > 0.1 Pa. “Not shown” indicates either that the spectral peaks are too close to diurnal harmonics or that the fitting to the Lorentzian function failed. See Supplementary for the comparison with the results from Obs-sampled-ERA5.
Figure 6 shows the meridional structures for amplitude and phase of R for Kelvin modes, as obtained with ISD/Buoy (red circles) and All-ERA5 (gray curves for amplitude and gray open circles for phase). See Supplementary for the results including Obs-sampled-ERA5 and also Table 1 for the difference between ISD/Buoy and Obs-sampled-ERA5. The results for k = 1 and k = 5 are of course identical to KL1 and KL5 as discussed in Section 3. Panels reading “Not shown” indicate that either that the spectral peaks are too close to diurnal harmonics or that the fitting to the Lorentzian function failed (see Fig. 9 of SH20 for details). Here, the amplitude has been multiplied by a factor of 
so that it represents a typical wave amplitude (recall that the normal mode index is normalized by its standard deviation).
Moreover, the depicted amplitudes and phases are the average of the four results obtained with the four lagged indices. Because the four indices and the resultant zonal distributions in R are lagged by a quarter cycle (e.g., Figs. 4g, i), R obtained from the index with lag of lπ/4 is shifted by 901/k °E in zonal direction so that the four results are in phase; after that, the average for harmonic coefficients and, then their amplitude and phase, are calculated. The phase shown in Fig. 6 thus simply represents the lag from the variation at the equator (positive and negative values indicate the variation precedes and follows that at the equator, respectively). Additionally, as in SH20, the corresponding Hough functions are fitted to the All-ERA5 results to see how close the meridional structure is to the theoretical mode structure (green curves in upper panels).
We find that both ISD/Buoy and ERA5 show a global Kelvin mode structure with maximum amplitudes at the equator and the phase being constant at almost all latitudes (note again that we did not assume any meridional structure). It is also discernible that the meridional extent in amplitude decreases with increasing zonal wavenumber. These structures are largely consistent with the theoretical Hough functions, clearly indicating that the normal-mode signals in ERA5 do exist in the raw observation data, i.e., ISD/Buoy in our case.
The extracted amplitude for KL1 is ∼ 10 Pa, commensurate with the value reported in previous studies that spectrally analyzed barometric measurements at tropical stations (Hamilton 1984; Matthews and Madden 2000). On the other hand, the components of k ≥ 3 have never been conclusively identified based on raw observation data, meaning that the present results are the first, ground-based evidence for these modes. Notably, their amplitude is fairly small (e.g., ∼ 1 Pa for KL5), but they are well defined in that the amplitude and phase have meridionally systematic structures. It is impressive that ERA5, synthesizing a necessarily imperfect numerical model with noisy observations, captures such small-amplitude, high-frequency pressure signals.
Despite the generally close correspondence between the two tested datasets, there are some differences, especially for larger zonal wavenumbers. In particular, the KL5 amplitudes deduced from ISD/Buoy are slightly smaller than those from ERA5 in the tropics. Considering the agreement between Obs-sampled-ERA5 and All-ERA5 (Fig. S1), this small discrepancy may not be attributable to the sampling inhomogeneity of the barometric network, but rather imperfections in ERA5.
b. Inertia-gravity modesFigures 7 and 8 (see also Figs. S2, S3) show the results for 1st gravest, equatorially symmetric and antisymmetric inertia-gravity modes, respectively (note that these are sometimes called “n = 1 inertia-gravity mode” and “n = 2 inertia-gravity mode,” respectively; Kiladis et al. 2009). These modes have very small amplitudes (∼ 1 Pa) and short periods, ranging from 12 hours (k = 1) down to 4 hours (k = 5). Nevertheless, our analysis method is well capable of extracting even these types of oscillations. For the 1st gravest symmetric modes (Fig. 7), for instance, the meridional structures in amplitudes and phases for the modes with k = −4, −3, −1, 1 show quantitatively good agreement between ISD/Buoy and ERA5 except for high latitudes. For the k = 4 mode, while the amplitude in ISD/Buoy is smaller than that in ERA5 especially in the tropics, the phase structure agrees fairly well. As was the case for k = 5 Kelvin mode (Fig. 6), the difference in amplitude is likely not due to the sampling inhomogeneity given the significant difference between ISD/Buoy and Obs-sampled-ERA5 (Fig. S2). The anti-symmetric mode (Fig. 8) reveals good agreement for k = −3, −2, −1, 0, 1, but for k = 3, −5 the ISD/Buoy amplitudes are smaller in comparison to ERA5 by more than 50 %.
As in Fig. 6 but for the first gravest, equatorially symmetric inertia-gravity modes for m = 1.
As in Fig. 6 but for the first gravest, equatorially anti-symmetric inertia-gravity modes with m = 1.
Figures 9 and 10 (see also Figs. S4, S5) show the 2nd gravest, equatorially symmetric and anti-symmetric modes (note that these are sometimes called “n = 3 inertia-gravity mode” and “n = 4 inertia-gravity mode”, respectively). Amplitudes decrease to about 0.7 Pa or less, and the periods are ∼ 1 hour shorter compared to the 1st gravest modes with the corresponding wavenumbers (cf. Fig. 1); the meridional structures are more complicated with more nodes in latitude (see green curves in Figs. 9, 10). Nevertheless, even these higher modes are detected in ISD/Buoy data, at least for the small zonal wavenumber components. For symmetric modes, we observe relatively tight agreement in both amplitude and phase for k = −2, −1, 0, while amplitudes in ISD/Buoy are somewhat smaller than those in ERA5 for k = −4, 2 (Fig. 9). Similarly, for anti-symmetric modes, good agreement is observed for k = −1 while for k = −3, 1, 4 amplitudes in ISD/buoy are too small, yet phases are in good agreement (Fig. 10).
As in Fig. 6 but for the second gravest, equatorially symmetric inertia-gravity modes with m = 1.
As in Fig. 6 but for the second gravest, equatorially anti-symmetric inertia-gravity modes with m = 1.
In summary, the signals extracted from ISD/buoy have a meridionally coherent structure that is consistent with the corresponding Hough function, although amplitudes tend to be smaller than those seen in ERA5 for some high zonal wavenumber (k) components. We reiterate that no previous study has obtained robust inertia-gravity mode signals including their horizontal structures based on ground-based data. The present work thus provides solid evidence for the existence of high-frequency, global inertia-gravity modes.
c. Rossby and Rossby-gravity modesFigures 11 and 12 (see also Figs. 6, 7) show the results for the symmetric, gravest Rossby mode and antisymmetric Rossby-gravity mode (note again that the eastward components of the latter is sometimes called as “n = 0 eastward inertia-gravity mode”). As noted in the Introduction, Rossby and westward Rossby-gravity modes have relatively low frequencies (cf. Fig. 1), and thus they can attain relatively large amplitudes. We indeed see that these signals are clearly extracted from ISD/Buoy data and their amplitudes and phases closely correspond to those in Obs-sampled-ERA5. For the k = −1 Rossby mode (“5-day wave”), the maximum amplitude (∼ 70 Pa) is consistent with that reported in previous studies (e.g., Madden and Julian 1973). Eastward Rossby gravity modes are also clearly detected, though again for the weak modes (e.g., k = 5) amplitudes in ISD/Buoy appear somewhat smaller than those in ERA5.
As in Fig. 6 but for the gravest, equatorially symmetric Rossby mode with m = 1.
As in Fig. 6 but for the Rossby-gravity mode with m = 1. Note that the eastward component of the latter is sometimes called as “n = 0 eastward inertia-gravity mode.”
It is worth mentioning that the meridional structure as deduced from ISD/Buoy and ERA5 is slightly different from the theoretical Hough function solutions for large, low-frequency k components. For example, the k = −3 and −4 Rossby modes and k = −5 Rossby-gravity modes exhibit a more compact structure in meridional direction (i.e., confined to equatorial region) in ISD/Buoy and ERA5, compared to the corresponding Hough function. These modes are more susceptible to background winds because of their slow phase speed (cf. Fig. 1) and thus the meridional structure may deviate from the theoretical prediction obtained under the assumption of no background winds. Recently, Ishizaki et al. (2025) generalized the classical linear theory to treat mean states with prescribed height-latitude distributions of zonal wind and temperature. By solving the resulting eigenvalue problem, they obtained the predicted frequencies and vertical and meridional structures of the normal mode oscillations. They showed that for a realistic mean state, the meridional structures of the Rossby and Rossby-gravity mode solutions differed somewhat from the corresponding classical theory solutions. Furthermore Ishizaki et al. determined that the deviation from classical theory for these modes is mainly attributable to the vertical mean flow shear (rather than the horizontal shear considered by e.g., Kasahara 1980). The Ishizaki et al. solutions allow the meridional structure of the modes to vary in the vertical and indeed the modes become more confined near the equator, being consistent with the present findings (Figs. 11, 12). This is likely related to the fact that the mode is confined near the equator at heights where Doppler-shifted frequency is reduced (in agreement with Lindzen 1970, and other earlier studies of equatorial waves).
4.2 Pekeris resonance (m = 2)Figure 13 (see also Fig. S8) shows the extracted meridional structures for k = 1 and 2 Kelvin modes for the Pekeris resonance (m = 2). We observe clear signals in ISD/buoy data that follow a Kelvin mode structure and peak at equatorial amplitudes of ∼ 5 Pa and ∼ 3 Pa for k = 1 and 2, about 50 % smaller than their counterparts with Lamb resonance (Fig. 6). This constitutes the first robust evidence for the existence of Pekeris normal modes in the real atmosphere relying purely on raw ground-based measurements. Notably, unlike the Lamb series (m = 1), such internal modes (m = 2) might be expected to be more susceptible to the effects of top boundary of the model used in the reanalysis procedure. The close agreement in amplitude between ISD/Buoy and ERA5 shown in Fig. 13, however, indicates that such effects may be negligible at least for ERA5 data. Žagar et al. (2022), using the Hough-mode expansion technique, showed that the energy in Kelvin mode takes a maximum for the components having an equivalent depth of ∼ 7 km. This may correspond in some part to the Pekeris resonance detected in this study, although Žagar et al. (2022) calculated the equivalent depths for a vertically bounded atmosphere and so their physical meaning was slightly different from our study.
As in Fig. 6 but for the Kelvin modes with m = 2 (Pekeris resonance).
Previous attempts to detect the global normal mode oscillations of the atmosphere only using station observations have had to cope with the limitations, particularly the irregular space-time distribution, of the available in situ data. Here we have adopted a simple regression analysis to extract globally coherent signals from barometric observations. The key point of our approach consists in regressing the raw pressure measurements (ISD/Buoy) onto a single modal index time series created by analysis of gridded reanalysis data. The method does not require homogeneity in temporal coverage and sampling among different stations, allowing records to differ in length and measurement frequency.
As a result, we have successfully identified the normal mode signals in the ISD/Buoy dataset not only for low-frequency modes such as Rossby and Rossby-gravity modes, but also for high-frequency modes such as Kelvin and inertia-gravity modes (even down to 2nd gravest modes). Despite not assuming any a priori horizontal dependence, a globally coherent, characteristic mode structure emerges from our computed regression coefficients, and the meridional structures obtained agree fairly well with the corresponding Hough functions. In addition, the same analysis was repeated for ERA5 itself (global data), yielding modal amplitudes and phases consistent with those in ISD/ Buoy data.
These findings corroborate the evidence (SH20) for a spectrum of global normal modes in the real atmosphere that match the expectations from classical tidal theory. Also, the results of the regression analysis indicate that the signals in ERA5 are generally in phase with, and of the same magnitude as, those revealed by analysis of raw barometric data. The agreement underscores the validity and usefulness of ERA5 for normal mode detection and investigation. Although the results are generally consistent between ISD/Buoy and ERA5, we find a few minor differences. In particular, amplitudes derived from ISD/Buoy are smaller than those of ERA5 for some high-frequency, high zonal wavenumber modes (note that the phase mostly display good agreement, though). In addition, the meridional structure deduced from ISD/Buoy and ERA5 deviates somewhat from the theoretical Hough function for high zonal wavenumber Rossby and westward Rossby-gravity modes. It is likely that background zonal winds alter the mode structure given that higher wavenumber modes are associated with smaller phase speeds (cf. Fig. 1).
Taken together, our results can be viewed as a unique evaluation of reanalysis data. As noted in the Introduction, it has been known that at least some reanalyses contain realistic depictions of the principal lunar air tide (∼ 10 Pa). The present study shows that ERA5 captures even much smaller global wave signals, e.g., ∼ 1 Pa for some inertia-gravity modes. Accurate representation of the pressure signal associated with the Pekeris internal resonance is particularly compelling and helps allay concerns about possible impacts of the upper boundary in the ERA5 product. Our simple technique that uses reanalyses to extract the signals from spatiotemporally inhomogeneous observations could be applied to any meteorological disturbances associated with identifiable spectral peaks such as equatorial waves. This may be a powerful tool to validate the detailed representation of other smallmagnitude fluctuations in atmospheric reanalyses.
ISD data are obtained by National Centers for Environmental Information though the website: https://www.ncei.noaa.gov/data/global-hourly/. Tropical buoy data are provided by the GTMBA Project Office of NOAA/PMEL through the web site: https://www.pmel.noaa.gov/tao/drupal/disdel/. ISPD data are obtained by the NSF National Center for Atmospheric Research (NCAR) Research Data Archive (RDA) through the website https://rda.ucar.edu/datasets/d132002/. ERA5 data (Hersbach et al. 2023) were downloaded from the Climate Data Store (CDS) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview).
Supplementary file contaions figures (Figs. S1–8) showing the meridional structure of each normal mode as derived from (1) ISD/Buoy, (2) Obs-sampled-ERA5, and (3) All-ERA5 data, as well as the Hough mode fitting for (3). Note that the corresponding figures in the main text (Figs. 6–13) do not include the results for (2).
This study was supported in part by JSPS Grant-in-Aid for Scientific Research (21K03661, 24K00706). We are grateful to Prof. Kevin Hamilton, Prof. Keiichi Ishioka, and Mr. Hideaki Ishizaki for their valuable comments and discussion on the manuscript. T. S. also thanks Prof. Keiichi Ishioka for providing a numerical code for solving Laplace’s tidal equation. We thank Dr. Nedjeljka Žagar and another anonymous reviewer for their constructive reviews which greatly helped to improve the manuscript.
