## Abstract

We used newly available ERA5 hourly global data to examine the variations of atmospheric circulation on global scales and high frequencies. The space–time spectrum of surface pressure displays a typical red background spectrum but also a striking number of isolated peaks. Some peaks represent astronomically forced tides, but we show that most peaks are manifestations of the ringing of randomly excited global-scale resonant modes, reminiscent of the tones in a spectrum of a vibrating musical instrument. A few such modes have been tentatively identified in earlier observational investigations, but we demonstrate the existence of a large array of normal mode oscillations with periods as short as 2 h. This is a powerful and uniquely detailed confirmation of the predictions of the theory of global oscillations that has its roots in the work of Laplace two centuries ago. The delineation of the properties of the modes provides valuable diagnostic information about the atmospheric circulation. Notably the amplitudes and widths of the normal mode spectral peaks contain information on the forcing mechanisms and energy dissipation for the modes, and the simulation of these properties for each of the many modes we have identified can serve as tests for global climate and weather prediction models.

## 1. Introduction

Studies of atmospheric dynamics typically deal with phenomena within some limited range of space and time scales, and the focus has generally been on those scales most directly connected to major weather phenomena, notably the synoptic scales (several days and horizontal wavelengths ~1000–10 000 km) and mesoscales (minutes to several hours and wavelengths ~1–1000 km) as shown in Fig. 1 (Orlanski 1975; Williams et al. 2017). In this paper we explore a much less studied spectral range of motions: those with high frequencies (periods ~2 h–2 days) and very large horizontal wavelengths (>5000 km). This region includes the familiar astronomically forced lunar and solar atmospheric tides and, as we explain further below, is expected to feature an array of normal mode oscillations of the global atmosphere, which will be the main focus of our investigation.

Like other mechanical systems such as the solid Earth (Nawa et al. 1998; Nishida et al. 2000), the atmosphere displays multiple quasi-monochromatic normal mode oscillations spanning a wide range of frequencies. There is a well-developed “classical” theory that treats the dynamics of an inviscid atmosphere on the rotating spherical Earth, but for an idealized case of small-amplitude perturbations about a motionless mean state characterized only by a mean temperature profile as a function of height, *T*(*z*). In this case the governing equations are separable in time and the zonal, meridional, and vertical directions, and so the solutions for, e.g., the perturbation geopotential can be written as the sum of individual modes with each represented as

where *λ* is longitude, *θ* is latitude, *z* is altitude (in log-pressure coordinate), *k* is the zonal wavenumber, *ω* is the angular frequency and *t* is time (e.g., Chapman and Lindzen 1970). Here, the meridional structures, or so-called Hough modes Θ(*θ*) are solutions to Laplace’s tidal equation (LTE) and *Z*(*z*) are solutions of a vertical structure equation (VSE). This formalism can be used to solve for the tidal response to a forcing with a fixed frequency [as described in detail in Chapman and Lindzen (1970)] but here we consider the eigensolutions to the unforced system, which should describe the possible free normal mode oscillations of the global atmosphere. In the following we generally refer to the free normal modes (i.e., modes in unforced system, or, homogeneous solutions) simply as normal modes. In this case the vertical structure equation may be solved first, which will give a separation constant as an eigenvalue, which can be expressed by a length *h*. For each *k* the Hough functions and mode frequencies *ω* are derived as the eigenfunctions and eigenvalues of the LTE, which expresses the dynamics of a shallow ocean of constant depth *h*. The separation constant *h* is commonly referred to as the “equivalent depth” of the atmosphere (e.g., Longuet-Higgins 1968; Chapman and Lindzen 1970).

The unbounded nature of the atmosphere in the vertical direction leads to a particular normal mode behavior. In contrast to bounded mechanical systems (simple examples are elastic cylinders or spheres), which have a set of normal modes that are complete in the sense that they can describe any inviscid motions, the atmosphere has a much smaller set of normal modes, since for a mode to ring freely requires that there be only horizontal energy propagation. This effectively limits such modes relevant for the lower atmosphere to be those trapped against the atmosphere’s lower boundary (what are referred to as Lamb waves). In the classical theory this issue is manifested in the VSE as the formulation of a radiation condition for the nominal “upper boundary” of the atmosphere. For the special case of an isothermal mean temperature structure the only normal mode solution (eigenfunction) to the VSE is well known to be the familiar Lamb wave, which has $Z\u2061(z)\u221de\kappa z/H$, where *H* is the atmospheric density scale height and *κ* is the ratio of the atmospheric gas constant to the specific heat at constant pressure. The equivalent depth for the Lamb wave solution is *Hγ*, where *γ* is the ratio of specific heat at constant pressure to the specific heat at constant volume. For *T* = 283 K, *H* = 8.3 km and the equivalent depth for the Lamb wave solution is *h* = 11 km. For more general *T*(*z*) the mathematical problem of solving for the eigenvalues and eigenfunctions of the VSE can be attempted with numerical methods, although there is some subtlety involved in appropriately imposing the radiation condition at the “top” of the atmosphere. Salby (1979) discusses this problem in detail and it is clear that for realistic *T*(*z*) profiles the only eigensolution relevant for the lower atmosphere is a solution trapped against the surface, which is a slightly modified version of the idealized Lamb wave. The computed equivalent depth obtained by Salby is about 9.6 km. Salby (1979) also discusses the possibility of one or two other solutions that correspond to vertical structures ducted in the mesosphere, but thus far no one has found evidence for such ducted global normal modes in the atmosphere and Salby (1980) concluded that such ducted mode solutions are of “dubious significance.”

Earlier Taylor (1929) noted that the atmospheric pressure pulse from the Krakatau eruption of 1883 could be tracked in microbarograph data as it circled Earth at least three times. He noted the observed phase speed of this wave pulse corresponded to that of a shallow ocean of depth roughly 10 km, suggesting that the dynamics of purely horizontally propagating waves in the atmosphere is analogous to that of a shallow water ocean of 10-km depth. More precisely, Chapman and Lindzen (1970) note that Taylor’s estimated phase speed for the pulse along the equator would correspond to *h* = 10.4 km.

The complete solutions to the LTE were obtained by Longuet-Higgins (1968) and are now straightforward to compute for any equivalent depth. Figure 2 shows our calculation of the normal mode frequencies for an atmosphere with an equivalent depth of 10 km, while Fig. 3 shows the corresponding Hough functions. The modes are classified into lower-frequency, westward-propagating Rossby modes, the eastward-propagating Kelvin modes, higher-frequency gravity modes, and the mixed Rossby–gravity mode. Since the classical theory ignores the mean winds, the theoretical predictions are expected to be most accurate for waves with zonal phase speeds much larger than the typical zonal mean wind speeds in the atmosphere. The zonal phase speeds of the Kelvin and gravity modes exceed 300 m-s^{−1}, while the Rossby modes have much slower phase speeds and so should be more distorted from the predictions of classical theory.

This beautiful theory produces the array of predicted frequencies and wavenumbers shown in Fig. 2. How will these predicted normal modes be manifest in actual atmospheric data? The space–time variability of the atmosphere can be characterized by variance spectra of observed meteorological fields and observational studies have shown that such spectra are generally fairly smooth and have an overall red behavior with variance density strongest at low frequencies and small wavenumbers. Even without detailed knowledge of all the relevant excitation mechanisms we can still hypothesize the normal mode oscillations to result in peaks of variance enhanced over this background. The differences of the observed central frequencies of these peaks from predictions should reflect the deviation of the real atmosphere from the simplified system described in the classical theory.

Evidence for these theoretically predicted modes in atmospheric observations has been sought in various studies over the last half century. Following pioneering observational studies of traveling large-scale waves by Deland (1965) and Eliasen and Machenhauer (1969), Madden and Julian (1972) found evidence for a westward-propagating zonal wavenumber one with period of about 5 days through Fourier analysis of gridded global maps of the surface pressure available once per day. They suggested that this oscillation corresponded to the gravest wavenumber-1 Rossby mode in Longuet-Higgins (1968) LTE solutions. Most later work related to atmospheric normal modes continued in this vein and dealt with low-frequency modes, including Madden and Julian’s “5-day wave” and other prominent modes now referred to as the “4-day wave,” “10-day wave” Rossby modes, etc. (Hirota and Hirooka 1984; Salby 1984; Madden 2007, 2019), which could be detected by global gridded surface analyses and satellite analyses with daily or twice-daily resolution.

Matsuno (1980) performed what he called a “trial search” for the higher-frequency Kelvin and gravity normal modes by computing spectra of long time series of hourly surface pressure at a few individual stations. He found an isolated peak in the power spectra at each station around a period of about 33 h, which he noted corresponded to the predicted frequency of the wavenumber-1 Kelvin mode for an atmosphere with 10-km equivalent depth. Hamilton (1984) and Hamilton and Garcia (1986) further investigated surface pressure spectra from individual stations and tentatively identified this “33 hour wave” as well as a handful of gravity modes with periods as short as 11 h. They also identified a peak in the surface pressure spectra that possibly corresponded to the “2 day wave” Rossby–gravity mode, which has itself been potentially linked to observed fluctuations in the winds at mesopause heights (e.g., Muller 1972; Gu et al. 2013). Matthews and Madden (2000) examined the 33-h wave in global reanalyses available four times per day and helped confirm the identification of this variance peak with the wavenumber-1 Kelvin wave normal mode. They also studied the vertical structure of the oscillation and showed that it was consistent with the expectation of a somewhat modified surface-trapped Lamb wave. Shved et al. (2015) have recently tentatively identified a limited number of gravity modes with peaks in spectra of single-station pressure records.

It should be noted that the results of the classical theory have also been loosely applied in recent years in a rather different context (i.e., for low-frequency motions appearing in a “convectively coupled” system, not the free modes of the dry atmosphere that are the main focus of this paper). Specifically, it has been found that at lower frequencies the variance of some atmospheric fields in the tropics have broad maxima that roughly fit the LTE solutions for much smaller equivalent depths; specifically peaks corresponding to *h* ~ 200 and *h* ~ 25 m have been found in several observational studies (e.g., Takayabu 1994; Wheeler and Kiladis 1999; Kiladis et al. 2009). The *h* ~ 200-m peaks exist, not due to the global resonance phenomenon, but due to a tendency of the vertical scale of deep convection in the tropics to establish a dominant vertical wavelength of forced equatorial waves (e.g., Salby and Garcia 1987; Horinouchi and Yoden 1996; Kiladis et al. 2009). The very low-frequency *h* ~ 25-m peaks are thought to result from systematic coupling between the tropical convection and the motion field, resulting in prominent moist convectively coupled waves (Kiladis et al. 2009). Note that the gravest modes for these small *h* values lie in a frequency range where they can be detected in standard global gridded analysis data based on conventional and satellite observations; by contrast, the global free modes, which are the main focus of the present study, lie in the much higher-frequency regions (Fig. 2) and thus have so far been challenging to observe.

The recent release of a long record of hourly global gridded fields as the ERA5 dataset (Hersbach et al. 2019) opens the possibility of a much deeper examination of the high-frequency global normal modes of the atmosphere. In our study we employed 38 years of the ERA5 data to compute the zonal wavenumber–frequency power spectrum of surface air pressure as well as of the 3D geopotential. We particularly focus on the tropical region, where the Hough functions for all modes (except for higher-order Rossby modes) retain significant amplitude (see Fig. 3) and where the signal to noise to detect short period oscillations is highest [see Zwiers and Hamilton (1986) for a detailed discussion of this signal-to-noise issue in the context of atmospheric tides]. Our results will demonstrate the observability of a very large array of the predicted global normal modes of classical theory. While these modes have only a small energy content, we will argue that the delineation of their properties is valuable as a diagnostic of the atmospheric circulation. Notably the amplitudes and widths of the normal mode spectral peaks contain information on the forcing mechanisms and energy dissipation for the modes, and the simulation of these properties for each of the many modes we have identified can serve as useful tests for global climate and weather prediction models.

The remainder of the manuscript is organized as follows. Section 2 describes the data and methods employed, while section 3 presents the basic results for the zonal wavenumber–frequency spectra. Section 4 discusses the shape of the spectral peaks corresponding to normal modes. Conclusions are summarized in section 5.

## 2. Data and analysis methods

### a. Data

We used ERA5, the latest atmospheric reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA5 was produced using 4D-Var data assimilation with the ECMWF’s Integrated Forecast System (IFS), which has 137 hybrid sigma–pressure model levels in the vertical with top level at 0.01 hPa. An artificial sponge layer is used in IFS to damp wave motions that propagate upward into the region near the model top (Shepherd et al. 2018). In this study, we mainly analyzed the hourly surface pressure dataset from 1979 to 2016 with horizontal resolution of 1°, while geopotential data on 37 pressure levels were also used to examine vertical variations. For more details about ERA5, see https://confluence.ecmwf.int/display/CKB/ERA5+data+documentation.

The unprecedented 1-h interval between reanalysis grids in ERA5 allows global space–time spectral analysis to investigate the high frequencies anticipated for the Kelvin and gravity normal modes. Of course, reanalyses are not raw observations, and it is reasonable to ask how well ERA5 can be expected to represent these normal modes. The traditional view is that reanalyses can adequately represent motions with large horizontal and vertical scales (both satisfied by the normal modes considered here). The somewhat novel aspect introduced in our investigation are the relatively high-frequency global-scale motions considered—how well are such motions represented in reanalyses? We have one set of special test cases for large-scale high-frequency (*T* ~ 12 h) motions, namely, the solar and lunar tides. The solar tide in surface pressure has been shown to be well represented in modern reanalyses (e.g., Díaz-Argandoña et al. 2016; Hamilton and Sakazaki 2017; Sakazaki et al. 2018). Also very relevant is the case of the lunar semidiurnal tide (L2), particularly as observed in the surface pressure. Since hourly pressure data at any station can be composited over long periods to isolate L2, we have a reasonable picture of the real observed L2 distribution from long records at ~100 individual barometric stations. One can compare the L2 then determined from the reanalyses data and find that there is agreement between analyses and the pattern obtained from compositing raw barometric data (e.g., Kohyama and Wallace 2014; Schindelegger and Dobslaw 2016). The L2 tide is exceptional in that the dynamical models used for reanalyses have not included the L2 gravitational forcing, so the tide in the reanalysis fields must come only from the actual data assimilated. So this shows that the reanalysis system is capable of taking the observed data and producing a good representation of the actual large-scale, high-frequency L2 tide. So it seems reasonable to suppose that the reanalysis will work well in representing other very large-scale, high-frequency (periods of a few hours) wave modes that we identified in our study.

While most radiosonde data are available only twice per day, ERA5 will assimilate hourly observations of surface fields from ~5000 aviation meteorology support stations (“METAR” stations) throughout the world as well as surface observations from another ~10 000 synoptic meteorological stations (SYNOPS) that report four times per day or more frequently (Haiden et al. 2018). Even with gaps over many parts of the oceans, these data are adequate to resolve large-scale waves such as zonal wavenumbers 1–4 in the tropics, which we analyze in our study. The reanalyses also assimilate very frequent asynoptic satellite radiometer data, and the Lamb wave modes we examine have very long vertical scales and so may be well resolved by such observations. It is also very notable, as explained in the introduction, that some of the Kelvin and gravity modes that we will study in our analysis of ERA5 data had been tentatively identified in earlier studies involving spectral analysis of time series of raw hourly single station barometric data (Matsuno 1980; Hamilton 1984; Hamilton and Garcia 1986). Our analysis enables us to clearly confirm that the earlier peaks are indeed the theoretically expected normal modes by revealing their spatial structure. So we have several examples of the same modes identified in the ERA5 and in raw data.

Note also that we do not expect the dynamical model in the IFS to generate spurious large-scale high-frequency motions that could somehow be interfering with the modes that appear due to the assimilation of real data. This might be a concern for very old-fashioned GCMs with low tops where spurious modes might be introduced by reflection of waves from the upper boundary (Lindzen et al. 1968; Kasahara 1976; Zwiers and Hamilton 1986; Žagar et al. 2015). However, the dynamical model used in IFS has its top full level at 0.01 hPa with a substantial damping region below the top; this should eliminate any concerns about reflections and it also ensures that the upper boundary will have negligible influence on the Lamb wave modes we discuss. We will show below that the signals of the normal modes for *h* = 10 km are clearly observed, but we did not find evidence for other internal global modes.

### b. Zonal wavenumber–frequency spectrum

A 2D-Fourier transform is applied for tropical hourly surface pressure data, which is a function of longitude and time. Data along single latitudes were first analyzed. The calculation was performed first for individual years (using the first 365 days each year, i.e., 365 × 24 hourly values even for leap years), after removing the linear trend over the year at each grid point. The spectrum, *q*(*k*, *f*), of the grid-point pressure values, *y*(*λ*_{i}, *t*_{j}), has been calculated as a function of zonal wavenumber *k* and frequency *f* (cpd):

Here, *k* = −*N*/2, …,−1, 0, 1, …, *N*/2 and *f* = 0, Δ*f*, …, *M*Δ*f*/2, where Δ*f* = 1/365 cpd. The results are shown for a power spectral density *P*(*k*, *f*) (Pa^{2} cpd^{−1} wavenumber^{−1}) which is defined as

The computed spectra (with individual year data) were averaged over the 38 years to reduce the noise. The single latitude results showed many isolated peaks that seemed likely to be identifiable as global normal modes, and our further analysis was guided by the expectations of theory for these modes.

The theoretical normal modes are meridionally coherent and each mode is either symmetric or antisymmetric about the equator (Figs. 2 and 3). To obtain the hemispherically “symmetric” results, the pressure at each longitude was averaged between 20°S and 20°N and then the spectral calculation was applied to the 38 years of data. To obtain the “antisymmetric” results the pressure average for 0°–20°N minus that for 20°S–0° was computed at each longitude and the 2D spectral calculation was repeated for this quantity. The latitude averaging procedure described here acts to emphasize meridionally coherent motions such as the global normal modes, and indeed we will see below that the signal to noise for the peaks of our meridionally averaged pressure variable is improved compared to the spectrum of raw pressure values for a single latitude.

### c. Estimation of observed spectral parameters

We examine the characteristics of the spectral peaks that correspond to the global normal modes by fitting the Lorentzian function to the frequency spectrum peaks for each zonal wavenumber. Detailed procedures are described below.

First, the background spectrum is estimated by fitting the AR(1) spectrum (“red noise” spectrum) to the frequency spectrum at each zonal wavenumber (Kikuchi, 2014). The spectrum of red noise is (Gilman et al. 1963)

where *σ* and *ρ* are both fitting parameters and *f* is the frequency (cpd); Δ*t* is the temporal resolution (=1/24 days); and *ρ* corresponds to the lag-1 autocorrelation coefficient for the time series of each zonal wavenumber component (again, this is a fitting parameter in our analysis). To keep weighting parameters relevant to high frequencies that have much less power, the fitting is done after taking the logarithm of the spectrum. Figure 4 shows, as an example, the frequency spectrum at eastward-propagating zonal wavenumber 1 for the meridionally averaged symmetric component (black curve), and the estimated background spectrum for this (red curve). This analysis has been done for zonal wavenumbers between −10 and +10 and those background spectra are shown in Fig. 5.

Next, around each identified peak the following Lorentzian function is fitted to the residual spectrum after the fitted background has been removed:

where *f* is frequency (cpd), and (*A*, *d*, *f*_{o}, *p*, *q*) are the fitting parameters. Here, *f*_{o} is the central frequency, *A* is a measure of spectral peak value, while *d* (cpd) is a measure of spectral width. Especially *d* is the measure of dissipation in a simple linear, weakly damped harmonic oscillator excited by a white-noise forcing. The actual fitting is made for the frequency region (*f*_{o} − *w*, *f*_{o} + *w*), where *f*_{o} is the central frequency and *w* is set to 0.2 cpd for the gravity modes and 0.15 cpd for the other modes. [More precisely, *f*_{o} is an unknown parameter so that the first iteration is done in the frequency region (*f*_{t} − *w*, *f*_{t} + *w*), where *f*_{t} is the theoretical frequency, to determine a tentative value of *f*_{o}$(fo\u02dc)$; then the fitting process is repeated in the region $fo\u02dc\u2212w,fo\u02dc+w$.] We consider the linear term (*pf* + *q*) to represent the “local” background spectrum near the spectral peaks, which is not completely captured by the red noise background spectrum determined by Eq. (4) for the 0–12-cpd range, where 12 cpd is the Nyquist limit of the hourly reanalysis output (e.g., Fig. 4) (again, *p* and *q* are also fitting parameters). Also, the frequency region around the solar tidal harmonics (*n* ± 0.05 cpd, where *n* is an integer) and lunar semidiurnal harmonic (24/12.42 ± 0.02 cpd) are not considered for the fitting.

### d. Regression analysis for deriving the meridional and vertical structure

As described above we use zonal wavenumber–frequency spectral analyses of tropical meridional surface pressure averages to identify the signature of the anticipated normal modes. To examine the meridional and vertical structures of these modal oscillations, we first establish a time series of an index for each mode based on narrow-band filtering of the hourly tropical pressure values. Then the latitudinal and vertical structures are determined based on regression with the index. For this analysis, we analyzed data for only 7 years between 2010 and 2016 to reduce the computation cost.

Specifically, we begin with the surface pressure averages described above—the average of 20°S–20°N is used for symmetric modes and the difference between 0°–20°N and 20°S–0° is for antisymmetric modes. The index time series of each free oscillation mode is then defined as the filtered time series at 0° longitude. The filtering is designed such that for each mode only spectral coefficients with the corresponding zonal wavenumber *k* and within the frequency range (*f*_{o} − 3*d*, *f*_{o} + 3*d*) are retained while the others are simply set to zero [*f*_{o} and *d* are estimated by fitting Eq. (5) to the raw spectrum; see Table 1 for these values] (cf. Kiladis et al. 2009). Note that for the Lorentzian function, the spectral energy at *f*_{o} ± 3*d* drops to 0.1 times the peak value.

The analysis has been applied for the modes that are well defined and are not contaminated by tidal harmonics. That is, putative normal modes are not considered for further analysis if any of the following occur: 1) one or more of the solar tidal harmonics and/or lunar semidiurnal are close to *f*_{o} (e.g., symmetric gravity mode for *k* = 5); 2) the spectral width appears to be too small for our procedure (<0.01 cpd; e.g., *k* = 1 of second symmetric gravity); 3) that the spectral fitting fails (e.g., the *k* = −5 Rossby mode) (see Fig. 11).

For each normal mode, the surface pressure data at all latitudes but filtered to include only the corresponding zonal wavenumber |*k*| (no attempt is made to filter out the eastward- and westward-propagating waves separately) are then regressed onto the index. Then, the regression coefficient *r*_{h} can be expressed as a function of longitude *λ* and latitude *θ* as

where *δ* is the phase relative to the index. The amplitude is finally rescaled by multiplying *R* by a factor $2\sigma t$, where *σ*_{t} is the standard deviation of index time series; that is, the amplitude shown later denotes the typical magnitude of each free oscillation (units: Pa). Note that we do not assume any a priori latitudinal structure for this analysis. The observed structure is then fitted by the theoretical Hough function of each mode as

where $R\xaf$ and $\delta \xaf$ are constant amplitude and phase, respectively. In our results, $2\sigma tR\xaf\Theta \u2061(\theta )$ (e.g., fitted amplitude) is shown to demonstrate how well the theoretical Hough function reproduces the observed meridional structure.

For the vertical structure, meridionally averaged geopotential data filtered to include the appropriate zonal wavenumber |*k*| (again, the average of 20°S–20°N is used for symmetric modes and the difference between 0°–20°N and 20°S–0° is used for antisymmetric modes) are regressed onto the same index time series. The regression is performed with data at 37 pressure levels spanning 1000 and 1 hPa and the regression coefficient *r*_{υ} can be expressed as a function of *λ* and pressure level *z* as

where *δ*_{g} is the phase relative to the index based on surface pressure. As for the case of the meridional structure, we finally multiply *R*_{g}(*z*) by a factor $2\sigma t$ to represent the typical amplitude of each mode (units: m^{2} s^{−2}). For phase, positive (negative) slope indicates the eastward (westward) tilting structure [see Eq. (8)], indicating an upward energy propagation for eastward (westward) modes.

## 3. Results

Before showing results for a meridionally averaged field, a typical result for the zonal wavenumber–frequency spectrum at a single latitude (10°N) is shown in Figs. 6a and 6b. Note that Figs. 6a and 6b display a red background spectrum but with a striking number of clearly identifiable isolated peaks superimposed, beyond the expected set of prominent lines corresponding to the astronomically forced solar (24 h and its higher-order harmonics) and lunar (12.4 h) tides. The spectral peaks for each zonal wavenumber occur in some clearly systematic distribution of frequencies, reminiscent of the tones in a spectrum of a vibrating musical instrument, and indeed we will show that these peaks correspond to the ringing of the free oscillations of the global atmosphere.

The theoretical normal modes are meridionally coherent and each mode is either symmetric or antisymmetric about the equator. So to emphasize the symmetric modes we computed spectra of the ERA5 pressure data averaged over 20°S–20°N, and then to emphasize the antisymmetric modes we computed spectra of the difference of the pressure averaged 0°–20°N and that averaged 0°–20°S. The results are shown in Figs. 6c–f. We again clearly see isolated peaks standing out from the background at many wavenumbers and frequencies. There is an impressive correspondence between the observed spectral peaks and the theoretical normal mode frequencies for *h* = 10 km (Figs. 6g–j). As reviewed in section 1 above, several individual modes have been tentatively identified in earlier observational investigations (mainly Rossby modes), but our study is the first to demonstrate the existence of the large array of normal mode oscillations with periods as short as 2 h. Both symmetric and antisymmetric modes are clearly responsible for the numerous peaks appearing in the results for a single latitude (Figs. 6a,b).

The regular-looking spacing of the peaks in Fig. 6 is quite striking. Of course, part of this effect comes from making a contour plot of a quantity that is only defined at discrete values of *k*. However, in the frequency domain the peaks are very well resolved in our long data record and we can even make inferences about the details of the shape of the individual peaks (Section 4 below). Also noteworthy is that the array of spectral peaks is quite apparent even in a simple unprocessed plot of the variance spectrum (Fig. 6), as long as a logarithmic contour interval is used so that the wide range of spectral power can be displayed together. In Fig. 7 we show the Fig. 6 results but now normalizing the local spectral variance by the fitted background spectrum (Fig. 5). This leads to an even clearer depiction of the normal mode peaks, at least in some ranges of spectral space. We regard the very beautiful Fig. 7 as a Fourier-space “portrait” of our atmosphere ringing like a bell! It may be worth noting that antisymmetric modes appear even more clearly in this figure than symmetric modes especially for gravity modes. This might be explained by the somewhat smaller background spectra for antisymmetric pressure variations (Fig. 5).

Frequency spectra for selected eastward and westward wavenumbers are shown in Fig. 8. We see that the spectral peaks for free oscillations are much broader than those of solar/lunar tides. Signatures of the antisymmetric Rossby modes and higher-order Rossby modes (e.g., 16-day wave) are not clearly apparent (left two spectra; see also Figs. 6 and 7), likely because their very low frequencies lead to a large distortion by extratropical mean winds (Salby 1981b); they are not considered in the following analyses.

With this simple space–time Fourier analysis we have demonstrated that Earth’s atmosphere is ringing with a broad array of discrete global-scale, free oscillations identifiable in long observational records. Over the wavenumber (0–5) and period [as small as 2 h (Nyquist frequency in the present study)] range we investigated, we found evidence for essentially the full expected array of normal modes for an atmosphere with *h* = 10 km and no evidence for modes with any other equivalent depth—except for the forced modes corresponding to small *h*, which we discuss next.

We will briefly show how the well-known equatorial waves believed to be directly related to tropical convection appear in the present analysis. Figures 7g–j additionally show the dispersion relationship for *h* = 200 m and *h* = 25 m. It is found that observed spectra also have broad, modest peaks around these two sets of dispersion relations (especially Kelvin, Rossby–gravity, and the gravest gravity mode); they are likely the signatures of equatorially trapped “dry” modes (*h* ~ 200 m) and “convectively coupled” modes (*h* ~ 25 m), respectively, associated with tropical convection and identified in earlier publications (e.g., Kiladis et al. 2009). These low-frequency convectively “forced” modes (not “free” modes) will not be examined further in this study.

To examine the meridional structure of the free oscillation modes (*h* = 10 km), we regressed the global surface pressure data (during the period 2010–16) onto an index time series of each free oscillation defined with the filtered tropical surface pressure data at 0° longitude (see section 2d). Figure 9 shows the meridional distribution of the amplitude and phase of the regression analysis, while the best-fit theoretical Hough functions derived from the LTE for each mode (Fig. 3) are superposed (red curves). The observed structures (black curves) show an impressive agreement with the theoretical Hough modes even for the fine meridional structure of second gravity modes. Symmetric gravity and Kelvin modes have a maximum at the equator, with an equatorially trapped structure that is narrower for the larger zonal wavenumber modes. Notably, only modes with zonal wavenumber 0 can retain a finite amplitude at poles (Chapman and Lindzen 1970). The agreement is less good for low-frequency Rossby and Rossby–gravity modes especially in the high latitude region, likely due to deviations from the idealized “classical” solutions attributable to the effects of mean winds.

We also investigate the vertical structure of each mode using geopotential data on 37 pressure levels during the period 2010–16. Figure 10 shows the vertical profiles of amplitude and phase determined by regression analysis (see section 2d). All modes except for some westward Rossby and Rossby–gravity modes show an exponential growth of amplitude with height at pressures smaller than 100 hPa. The phase is almost constant in vertical. These features are broadly consistent with the theoretical Lamb mode structure. By contrast, westward Rossby modes with *k* ≤ −2 and westward Rossby–gravity modes with *k* ≤ −1, which have relatively low frequencies, show an amplitude decay and a downward phase progression with altitude above the stratosphere, indicating upward energy propagation. Such a westward-tilting structure is consistent with those reported by previous studies of these Rossby modes (Hirota and Hirooka 1984; Salby 1984). These low-frequency Rossby modes deviate substantially from the “classical” solutions and their energy is likely not completely trapped near the ground but somewhat “leaky” in vertical (Salby 1984).

## 4. Analysis of individual spectral peaks

Thus far our results have demonstrated the expectations of classical theory for the free oscillations of the global atmosphere are broadly confirmed by the appearance of isolated spectral variance peaks at the predicted frequencies and zonal wavenumbers. Additionally, we have shown that the oscillatory motions near most of the spectral peaks have the theoretically predicted meridional and vertical dependence. Beyond this, however, there must be further information about the dynamics of the atmosphere contained in our results, and, in particular, in the details of the spectral peaks for each of our large array of identified modes. Below we discuss our empirical results for (i) the deviation of the spectral peak frequency from the predictions of classical theory and (ii) the width of each spectral peak, and we discuss the interpretation of our results within the expectations of linear wave theory. Note that in the following we consider only the first gravest modes for gravity wave modes.

Figure 11 shows the frequency spectra in the vicinity of each normal mode peak for each individual wavenumber. We found that each spectral peak can be reasonably fit with a Lorentzian function (blue curves; see section 2c for details).

In each case the empirically determined central frequency *f*_{o} is close to the theoretical frequency (red vertical lines) *f*_{t}; however, a small, but systematic, difference between the two is apparent especially for larger zonal wavenumber modes. Generally, *f*_{o} < *f*_{t} for westward waves (*k* < 0) and vice versa for eastward waves (*k* > 0), while *f*_{o} ~ *f*_{t} for the standing wave (*k* = 0). Figures 12a and 12b show the difference, Δ*f* = *f*_{o} − *f*_{t}, as a function of westward/eastward zonal wavenumber. We find for most of the modes an approximate linear relationship, with Δ*f* positive and negative for eastward and westward modes, respectively. Figures 12c and 12d show the difference for *f*_{t}, which is calculated with *h* = 9.5 and 10.5 km. Almost all observed frequencies are greater than theoretical ones when *h* = 9.5 km, and vice versa for *h* = 10.5 km. These findings again, consistent with aforementioned references and other observational studies, show that the choice of *h* = 10 km is a good one. The Δ*f* dependence on zonal wavenumber for *h* = 10 km (Figs. 12a,b) would be consistent with the Doppler shift in an idealized case of a mean state with constant relative mean angular momentum corresponding to eastward wind of 5–10 m s^{−1} at 30°. Of course, the real mean state is much more complicated with westward winds in the lower atmosphere in the tropics and strong eastward jets in the midlatitudes, and so it appears that the effect of the midlatitude jets dominates for these global modes. One oddity is that our determined Δ*f* for the Kelvin mode appears to be roughly independent of wavenumber, very close to zero for all *k*. The unique characteristics of the structure of the Kelvin mode (a pure equatorially trapped structure with no zero crossings or peaks off the equator, and almost no perturbation meridional wind) might be responsible, but this issue needs further investigation (cf. Kasahara 1980).

Figure 13 shows the spectral peaks after removing the “local” baseline of the background spectrum and being divided by the factor *A* [Eq. (5)]. The fitted Lorentzian functions are also shown in the right-hand panels. We see that, even though each peak is represented by a Lorentzian function quite well, the spectral width *d* differs among the modes.

The inviscid, adiabatic classical theory for the normal modes simply predicts the resonant frequency of each of the modes and has nothing directly to say about the width of each peak as realized in the real atmospheric circulation. We interpret our observed normal mode response in terms analogous to a simple linear harmonic oscillator subject to dissipation and excited by a broad spectrum of forcing. In this case we expect the spectrum of the response to be a Lorentzian function. As discussed earlier by Hamilton and Garcia (1986), the width of the spectral peak should be determined by (i) the strength of the dissipation, and possibly (ii) detuning of the peak frequency if the restoring force of the oscillator changes with time. In the atmospheric case we can expect there to be some detuning effect as the mean temperature and wind structure (that determine the resonant frequencies) will change with time. The most profound of such mean state changes are likely to be those of the seasonal cycle. We have repeated our spectral analysis for the surface pressure for four individual 3-month seasons (not shown) and found little seasonal change in the spectral width. So we conclude that the broadening of the resonant peaks seen in our spectral analysis is largely due to effective dissipation of the modes. Such dissipation could come from many mechanisms, including the surface friction and radiative cooling. For this real atmosphere case, however, the broadening might depend on other aspects of the flow that drain energy from the modal structure. This could include any upward propagation or “leakage” of energy to the upper atmosphere (where the energy would have to be ultimately dissipated) especially for low-frequency modes (see Fig. 10) or even the mixing effects of the large-scale flow. Thus, it seems reasonable to regard the spectral width for the modes to be indicators of the total energy dissipation acting on the modes. Specifically, our parameter *d* in Eq. (5) can be regarded as proportional to the reciprocal of the time scale for dissipation of disturbance energy (Hamilton and Garcia 1986). Figure 14 shows *τ* = 1/(4*πd*) days as function of zonal wavenumber. It seems that the characteristic damping time scale is the longest for small zonal wavenumbers (4–5 days for |*k*| = 0–1) and becomes shorter with increasing zonal wavenumbers (down to 1 day for |*k*| = 5). The damping time scale for the Kelvin wave with *k* = 1 is consistent with that reported by Hamilton and Garcia (1986).

The dependence of damping time scale on zonal wavenumber might indicate that horizontal mixing may be an important dissipation mechanism. This seems consistent with our finding that most modes have a nearly barotropic structure and so thermal damping may not be a leading dissipation process (Fig. 10). If the horizontal mixing is effective (and reasonably isotropic) instead, the damping scale should presumably depend on the “total horizontal” wavenumber including both zonal and meridional variations. Actually Figs. 14a and 14b show that the damping scale has a systematic variation between different types of modes even for the same zonal wavenumber (e.g., the scale for Kelvin mode is longer than that for gravity modes). To see this in detail, we define an “effective” total horizontal wavenumber for each of the Hough modes based on the global mean of the Laplacian acting on the mode ( appendix). As shown in Fig. 14c the damping scale is found to depend on that horizontal wavenumber reasonably well. This finding again suggests that horizontal mixing may be contributing to the dissipation, although other factors such as variability of background winds may contribute to such dependency (e.g., Salby 1981a).

By providing *globally integrated* measures of the dissipation of disturbance energy, our determined *d* values are an empirical diagnostic of dissipative processes that is unmatched by any other technique. While there is much to be understood about the processes that are actually dissipating the normal modes, the results here have important information and potentially represent a useful constraint on global model atmospheric simulations.

## 5. Conclusions

We analyzed the 38-yr surface pressure data from the latest reanalysis, ERA5, to identify the high-frequency free normal modes in the atmosphere. The wavenumber–frequency spectrum in the tropical region exhibits numerous spectral peaks, which we showed are the manifestation of a broad array of free normal modes predicted by classical theory with the equivalent depth of 10 km. The observed meridional and vertical structure of individual modes also agreed quite well with the corresponding Hough function and the Lamb mode, respectively. Therefore, our analysis has provided a striking basic confirmation of the classical theory of global atmospheric normal modes.

Our results provide another approach to estimate the equivalent depth of the atmosphere that supplements the earlier approaches of Taylor (1929) and Salby (1980) relevant for the horizontally propagating free modes. We found that *h* must be close to 10 km to provide the best fit overall for the full array of modes that we have identified. This value lies roughly midway between the 10.4 km based on observed propagation of the pressure pulse from the Krakatau eruption (Taylor 1929; Chapman and Lindzen 1970), and the 9.6 km computed by Salby (1979) as the eigensolution of a vertical structure equation including a realistic global mean temperature structure. The many clearly identifiable modes in Figs. 6 and 7 show just how complete is the mapping of atmospheric normal modes to the dynamics of Laplace’s idealized constant depth ocean.

While the classical solutions for *h* = 10 km provide the best fit with our observations, we showed there was a small difference between the observed frequencies and the classical predicted frequencies, and that this difference had a systematic dependence on the zonal wavenumber of the mode. This result is suggestive of a small modification of the expectations of classical theory due to a Doppler shift associated with a predominantly eastward mean flow.

While high-frequency modes we have identified have relatively little energy, the delineation of their properties can be valuable as diagnostics of the atmospheric circulation, just as observations of the free oscillations of Earth and the sun have been useful in advancing knowledge of the dynamics of those systems (Hill et al. 1996; Nawa et al. 1998; Nishida et al. 2000). Notably the amplitudes and widths of the normal mode spectral peaks may contain information on the forcing mechanisms and dissipation for the modes, and the simulation of these properties for each of the many modes we have identified can serve as tests for global climate and weather prediction models.

We have shown that the ERA5 dataset allows a much more complete characterization of the high-frequency normal mode variability of the atmospheric circulation than ever before. This, in turn, highlights the need to understand more about the modes, notably their mechanisms of excitation and dissipation. With more knowledge of these mechanisms we may be able to more precisely interpret the observations of normal mode spectra as diagnostics of the atmospheric circulation. The present authors are now working on these issues and hope to press forward with studies using general circulation models as well as the ever-improving global reanalysis datasets.

## Acknowledgments

T. S. thanks Dr. Keiichi Ishioka for providing a numerical code for solving Laplace’s tidal equation. We are grateful to three anonymous reviewers for their insightful comments on the original manuscript. T. S. was supported in part by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) through a Grants-in-Aid for Scientific Research (18K1361232). This is IPRC Publication Number 1446.

Data availability statement: ERA5 data are available through the Climate Data Store (CDS) (https://cds.climate.copernicus.eu/#!/home).

### APPENDIX

#### Effective Total Wavenumber

The results presented in Fig. 13 show that the spectral widths (and hence inferred characteristic damping times) vary substantially among the modes studied here. Figures 14a and 14b suggest that there is some tendency for the damping time scale to be shorter for higher zonal wavenumbers. If this indicates that possibly some of the energy dissipation occurs by large-scale horizontal mixing, then it suggests examining the dependence of the scale on the total horizontal wavenumber appropriate for each mode. The total horizontal wavenumber on the sphere is uniquely defined for a spherical harmonic on Earth, $Pnk$ [where *k* is order (zonal wavenumber) and *n* is degree], as $K=n\u2061(n+1)/a$, which satisfies Laplace’s equation $\u22072Pnk=\u2212K2Pnk$ (*a* is the radius of Earth). Unfortunately, individual Hough modes do not satisfy Laplace’s equation and some reasonable definition of an appropriate total horizontal wavenumber needs to be made. In our case we define for each of the Hough modes a horizontal wavenumber based on the global mean of the Laplacian acting on the mode:

where *ϵ* = 4*a*^{2}Ω^{2}/(*gh*) is the so-called Lamb parameter (e.g., Longuet-Higgins 1968), Φ is Hough function, and *χ* and *ψ* are velocity potential and streamfunction associated with the Hough mode, respectively, such that horizontal wind vector **v** of the mode is expressed as

Note that here these variables (and Laplacian) are nondimensionalized following Kasahara (1976). Figure A1 shows the total wavenumber calculated this way for the modes in our study.

## REFERENCES

*Mon. Wea. Rev.*

*J. Geophys. Res. Atmos.*

*Tellus*

*J. Atmos. Sci.*

*J. Geophys. Res. Atmos.*

*J. Meteor. Soc. Japan*

*J. Geophys. Res.*

*Quart. J. Roy. Meteor. Soc.*

*ECMWF Newsletter*

*Science*

*J. Atmos. Sci.*

*J. Meteor. Soc. Japan*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Climate Dyn.*

*Rev. Geophys.*

*Geophys. Res. Lett.*

*Mon. Wea. Rev.*

*Philos. Trans. Roy. Soc. London*

*Tellus*

*Bull. Amer. Meteor. Soc.*

*J. Atmos. Sci.*

*J. Meteor. Soc. Japan*

*J. Atmos. Sci.*

*Philos. Trans. Roy. Soc. London*

*Earth Planets Space*

*Science*

*Bull. Amer. Meteor. Soc.*

*Atmos. Chem. Phys.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Rev. Geophys. Space Phys.*

*J. Atmos. Sci.*

*J. Geophys. Res. Atmos.*

*Atmos. Oceanic Phys.*

*J. Meteor. Soc. Japan*

*Proc. Roy. Soc. London*

*J. Atmos. Sci.*

*Geophys. Res. Lett.*

*Geosci. Model Dev.*

*J. Geophys. Res.*

Atmospheric Tides: Thermal and Gravitational. D. Reidel, 200 pp