In the 1970s and 1980s, there was considerable interest in near-equatorial variability at periods of days to weeks associated with oceanic equatorial inertia–gravity waves and mixed Rossby–gravity waves. At that time, the measurements available for studying these waves were much more limited than today: most of the available observations were from scattered island tide gauges and a handful of short mooring records. More than a decade of the extensive modern data record from the Tropical Atmosphere Ocean (TAO)/Triangle Trans-Ocean Buoy Network (TRITON) mooring array in the Pacific Ocean is used to reexamine the internal-wave climate in the equatorial Pacific, with a focus on interpretation of the zonal-wavenumber/frequency spectrum of surface dynamic height relative to 500 decibars at periods of 3–15 days and zonal wavelengths exceeding 30° of longitude. To facilitate interpretation of the dynamic height spectrum and identification of equatorial wave modes, the spectrum is decomposed into separate spectra associated with dynamic height fluctuations that are symmetric or antisymmetric about the equator. Many equatorial-wave meridional modes can be identified, for both the first and second baroclinic mode. Zonal-wavenumber/frequency spectra of the zonal and meridional wind stress components are also examined. The observed wind stress spectra are used with linear theory of forced equatorial waves to provide a tentative explanation for the zonal-wavenumber extent of the spectral peaks seen in dynamic height. Examination of the cross-equatorial symmetry properties of the wind stress suggests that virtually all of the large-scale equatorial inertia–gravity and mixed Rossby–gravity waves examined may be sensitive to both zonal and meridional wind stress.
Equatorially trapped inertia–gravity waves were first identified in the ocean by Wunsch and Gill (1976) through analysis of Pacific Ocean island tide gauge records. Peaks in the frequency spectra at periods of about 3, 4, and 5.5 days were found to be common to islands over a large range of latitudes and longitudes. The oscillations at these periods were coherent with large-scale equatorial winds, but the available frequency spectra of winds showed either no peaks at these special periods or peaks that were much broader than the oceanic peaks. Furthermore, the coherence phase between wind and sea level was found to have 180° shifts near the oceanic spectral peaks, as would be expected from a resonant response. The peaks in the frequency spectrum of sea level were thus interpreted as oceanic resonances forced by the wind, and Wunsch and Gill showed that the distinct frequencies coincided with the modal-minimum frequencies predicted by linear equatorial wave theory for the lowest few meridional modes of baroclinic mode-1 inertia–gravity waves. They also showed that the latitudinal dependence of the power in the peaks at each period matched reasonably well with the latitudinal dependence of the pressure variance predicted for linear, baroclinic mode-1 inertia–gravity waves with basin-scale zonal wavelengths.
Equatorially trapped inertia–gravity waves and mixed Rossby–gravity waves (MRG) have since been observed in all of the major equatorial basins (e.g., Weisberg et al. 1979; Luther 1980; Eriksen 1980; Ripa and Hayes 1981; Eriksen 1982; Garzoli and Katz 1981; Chiswell et al. 1987; Weisberg and Hayes 1995; Gilbert and Mitchum 2001). Improvements on the analysis of Wunsch and Gill (1976) have been attempted. Luther (1980) used an expanded tide gauge dataset together with a more extensive analysis of the Pacific Ocean wind field to refine the results and arguments of Wunsch and Gill. Eriksen (1982) used cotemporal tide gauge records to estimate wavenumbers as well as frequencies. Using simultaneous records from two instruments having a zonal separation of about 300 km, Chiswell and Lukas (1989) were able to identify several zonal-wavenumber/frequency points along the MRG dispersion curve and two near the meridional mode-2 inertia–gravity wave dispersion curve of the first baroclinic mode.
Improvements over the analysis of Wunsch and Gill (1976) have been only incremental because of inadequacies in the spatial and temporal distribution of the available data. The error bars on the wavenumber estimates of Eriksen (1982) and Chiswell and Lukas (1989), for instance, were such that wavenumbers could have been either positive or negative. Consequently, some of the issues raised by Wunsch and Gill remain unresolved.
One such issue is the question of why the elevated energy levels were found near the frequency minima of the various inertia–gravity wave modes. Wunsch and Gill (1976) offered two hypotheses: first, that the ocean preferentially resonates at the wavenumber-frequency locus where the group velocity vanishes, or second, that the atmospheric forcing is strongest on basin scales, causing the oceanic response to be concentrated at small zonal wavenumbers. Given the available data, the frequencies predicted by the two hypotheses could not be distinguished from each other, leading Eriksen (1982) to estimate wavenumber–frequency spectra. However, the limited dataset still did not allow any distinction between an enhanced response at the wavenumber of vanishing group velocity and an oceanic response to basin-scale winds.
The mathematical analysis of Wunsch and Gill (1976) favored the low-wavenumber-forcing hypothesis as an explanation for the peaks they observed. Nevertheless, they suggested that one might expect to find elevated sea level variance at the wavenumber associated with vanishing group velocity, and this concept continues to be popular. In a companion paper (Durland and Farrar 2012), we show that a more detailed examination of the theoretical solutions suggests that we should not necessarily expect to see an enhanced sea level response at the wavenumbers associated with vanishing group velocity. One goal of the present study is to investigate Wunsch and Gill’s second hypothesis, that the character of the equatorial inertia–gravity wave field is a consequence of the large-scale nature of the wind forcing, by developing a more detailed characterization of the zonal-wavenumber/frequency spectra of dynamic height and wind stress in the equatorial Pacific Ocean.
Another issue raised by the early investigations concerns the extent to which baroclinic modes higher than the first might contribute to near-equatorial variability at periods of days to weeks. On theoretical grounds, Wunsch and Gill (1976) suggested that contributions to sea level variability from the higher modes should be minimal, and indeed, all of the nontidal spectral peaks that they observed could be explained by theoretical solutions for the first baroclinic mode. Luther (1980) detected a 7-day peak in sea level that coincided with the minimum frequency of a second baroclinic mode, first meridional mode inertia–gravity wave, but he ultimately opted against this interpretation based on sparse information about the meridional variation of the signal’s amplitude. Virtually all of the identifications of Pacific Ocean inertia–gravity waves mentioned above were for baroclinic mode 1. Using a gridded sea level product made from satellite altimetry data, Shinoda (2010) found relatively elevated power on the oceanic MRG dispersion curves for both baroclinic mode 1 and baroclinic mode 2. This suggests that the more extensive datasets now available may yield evidence of higher baroclinic mode inertia–gravity waves that were not detected previously.
The Tropical Atmosphere Ocean (TAO)/Triangle Trans Ocean Buoy Network (TRITON) mooring array, with its long, overlapping records from sites spanning the Pacific Ocean, allows a markedly improved estimate of the zonal-wavenumber/frequency spectrum of equatorial inertia–gravity waves and MRG in the Pacific. In this paper, we use daily-average estimates of surface dynamic height relative to 500 decibars (hereafter “dynamic height”) and surface winds from the TAO/TRITON array to estimate wavenumber–frequency spectra of dynamic height and wind stress in the 3–15-day period band. When compared to theoretical dispersion curves, the dynamic height spectra show evidence for resonant excitation of MRG and inertia–gravity waves of meridional modes 1 through 3, for both baroclinic mode 1 and baroclinic mode 2. Comparisons of the meridional dependence of spectral energy density with theoretical predictions support these interpretations. Symmetry filtering even provides evidence that the 5-day peak in the dynamic height spectrum contains roughly equal contributions from baroclinic mode 1, meridional mode 1, and baroclinic mode 2, meridional mode 2; modes of opposite cross-equatorial symmetry overlapping in the same portion of wavenumber–frequency space. Some of the inertia–gravity wave spectral peaks extend to zonal wavenumbers distinctly removed from the dispersion curve minimum, calling into question the necessity of a vanishing group velocity for producing the elevated energy levels.
Wavenumber–frequency spectra calculated for zonal and meridional wind stresses allow us to confirm and refine two important inferences made by Wunsch and Gill (1976) and Luther (1980). Over most of the period range of interest, both the zonal and meridional winds of the equatorial Pacific (8°S–8°N) are dominantly symmetric about the equator, and the wavenumber content is indeed concentrated in relatively narrow, low-wavenumber bands. The symmetry of the winds about the equator suggests that symmetric dynamic height signals should be forced primarily by the zonal winds, and antisymmetric dynamic height primarily by the meridional winds (e.g., Wunsch and Gill 1976; Durland and Farrar 2012), although we show that even small amounts of asymmetry in the winds can complicate the situation. Nonetheless, a simplistic comparison of the oceanic and atmospheric spectra based on the assumption that the wind field is perfectly symmetric about the equator does show that the most energetic peaks in dynamic height generally occur at wavenumbers and frequencies where the relatively energetic parts of the wind stress spectrum overlap with the free oceanic modes that can be resonantly excited.
While the meridional spacing of the TAO/TRITON mooring array cannot adequately resolve the meridional structure of oceanic equatorial waves and is thus not suitable for a detailed study of the forcing-response relationship, the zonal-wavenumber resolution is sufficient to demonstrate that vanishing zonal group velocity is not a requirement of oceanic resonance and that the wavenumber content of the resonance appears to be largely determined by the wavenumber content of the wind stress. The spectra also show strong evidence for second baroclinic mode inertia–gravity wave variability that has apparently not yet been observed and identified as such. This is perhaps not surprising, because the second baroclinic mode has oppositely phased pressure extrema at the surface and 500-m depth, making the surface dynamic height relative to 500 decibars a quantity that is almost perfectly “tuned” for detection of the second baroclinic mode (section 2).
2. Theoretical background
The strong equatorial currents are known to have a significant effect on the dispersion relations and meridional structures of equatorial Rossby waves (e.g., McPhaden and Ripa 1990; Chelton et al. 2003; Durland et al. 2011), but less is known about the effects of the currents on the higher-frequency MRG and inertia–gravity waves. The work that has been done on these higher-frequency (and higher-phase-speed) waves suggests that they are less sensitive to the equatorial currents. The Doppler shift becomes relatively less important with increasing phase speed, and the modification of the background potential vorticity gradient by the second meridional derivative of the mean currents, which has a large effect on Rossby waves (Durland et al. 2011), presumably has less of an effect on the inertia–gravity wave dynamics.
In a 1.5-layer model, McPhaden and Knox (1979) showed that the effect of the equatorial Pacific currents on the dispersion relations of low-wavenumber MRG and inertia–gravity waves was relatively small and could be represented largely by a constant frequency shift. The relative frequency change was a few percent at most for the lower-frequency MRG, and even this discrepancy might have been removed by a different choice for the mean layer depth (which is a somewhat arbitrary quantity in the presence of a complex current system and equatorially trapped waves: see Durland et al. 2011). McPhaden and Knox (1979) also found that, although the meridional structures of the waves’ zonal velocity were distorted by the mean zonal currents, the meridional structures of the waves’ pressure and meridional velocity were only minimally affected (also see McPhaden 1990). We will therefore interpret our results under the assumption that the dispersion relations and meridional structures of the pressure signals associated with low-wavenumber MRG and inertia–gravity waves are reasonably well described by linear theory. The viability of this assumption will be tested by comparisons of the theoretical predictions with the observations.
The dispersion curves and meridional structures associated with a particular baroclinic mode are dependent on the nonrotating gravity wave speed associated with that vertical mode through the vertical eigenvalue problem (e.g., Moore and Philander 1977). We estimated the gravity wave speeds for baroclinic modes 1 and 2 (Fig. 1) following Chelton et al. (1998), but using the time-averaged, mapped hydrographic data of the World Ocean Atlas (WOA) 2005 (Locarnini et al. 2006; Antonov et al. 2006). There are spatial and temporal variations in the hydrographic structure of the equatorial Pacific and the associated gravity wave speeds, but we will use 2.8 m s−1 (baroclinic mode 1) and 1.7 m s−1 (baroclinic mode 2) as nominal values representative of the time-mean state of the central equatorial Pacific.
The left panel of Fig. 2 shows the linear equatorial-wave dispersion curves over a slightly broader zonal wavenumber–frequency (k–ω) range than the one we will investigate, for baroclinic mode 1 and baroclinic mode 2. To simplify subsequent discussions, we will use an abbreviated notation to refer to the meridional and baroclinic mode numbers that uniquely identify a particular equatorial wave mode; for instance, we will identify baroclinic mode 2, meridional mode 0 as bc2m0. The curves in the bottom left corner, at frequencies less than about 0.03 cycles per day (i.e., periods longer than 30 days), are for the Rossby modes, with which we will not be concerned. The curves that have frequency minima (at the intersections with the red curve) are for the various modes of inertia–gravity waves (m = 1, 2, 3, …), and the two curves with no frequency extrema are for the MRG (m = 0) of baroclinic modes 1 and 2.
The right panel of Fig. 2 shows the latitudinal structures of the meridional velocity V(m)(y), zonal velocity U(m)(y), and pressure P(m)(y) for unit-amplitude free waves of baroclinic mode 1, meridional modes m = 0 through 4, all calculated at k = 0.1 The corresponding structures for baroclinic mode 2 are compressed toward the equator by about 20%. We retain the notation of Durland and Farrar (2012) in which the parentheses around the mode number subscript indicate a function of y (latitude) as opposed to a Hermite expansion coefficient. In keeping with the common convention (e.g., Moore and Philander 1977), V(m) is real while U(m) and P(m) are imaginary. This just means that the U(m) and P(m) structures seen in Fig. 2 are in phase with each other and they lag the corresponding V(m) structure by a quarter of a period. The vertical black lines in the right panel are at the nominal latitudes of the TAO/TRITON moorings (8°S, 5°S, 2°S, 0°, 2°N, 5°N, and 8°N)– the difficulty of using the TAO/TRITON measurements to resolve the meridional structure of the equatorial waves is apparent.
An important point to note in Fig. 2 is that for odd numbered modes U(m) and P(m) are symmetric about the equator while V(m) is antisymmetric. The opposite is true for the even numbered modes: U(m) and P(m) are antisymmetric while V(m) is symmetric. Because we will analyze dynamic height (proportional to pressure), we will use the symmetry of P(m) as the designation of the wave symmetry: the odd-numbered meridional modes are referred to as being symmetric and the even-numbered modes as antisymmetric.
A second important point illustrated in Fig. 2 is that for P(m) at small wavenumbers, the local extrema closest to the equator are larger than the more poleward extrema. Consequently, for meridional modes m ≤ 4 the power in dynamic height variability (proportional to |P(m)|2) is largely confined to within 7°–8° of the equator. The power in the dynamic height signals associated with these waves at the 8°N and 8°S TAO/TRITON moorings will be small compared to the power closer to the equator and away from zeros of P(m) (see also profiles of |P(m)|2 in Figs. 9–10). The structure of |P(m)(y)|2 is wavenumber dependent, but the concentration of energy in the equatorward extrema applies throughout the band of wavenumbers examined here.
The proportionality between the amplitude of a resonantly forced wave and the wind stress can be expressed as (e.g., Durland and Farrar 2012)
V(m)(y) and U(m)(y) are the meridional structures of meridional and zonal velocity associated with a unit-amplitude free wave (see Fig. 2 and footnote 1), while p(m)(y) is the pressure signal of the forced wave. The equatorially symmetric and antisymmetric parts of the meridional wind stress are τys and τya, respectively, and the corresponding parts of the zonal wind stress are τxs and τxa. The proportionality is wavenumber and frequency dependent, but (1)–(2) express all of the forcing dependence. Basically, the efficiency with which the wind excites a particular oceanic mode depends on the projection of the meridional structure of the wind stress onto the meridional structure of the free wave’s velocity. The symmetry properties of the free waves allow the separation of the forcing expressions for symmetric and antisymmetric oceanic modes seen in (1)–(2): symmetric modes (m odd) are forced by symmetric zonal winds and antisymmetric meridional winds; antisymmetric modes (m even) are forced by symmetric meridional winds and antisymmetric zonal winds.
When the equatorial winds are symmetric about the equator (τya = τxa = 0), the symmetric oceanic modes (m odd) are forced only by the zonal wind stress, and the antisymmetric oceanic modes (m even) are forced only by the meridional wind stress. This greatly simplifies the matter: we can make predictions about a particular oceanic mode with knowledge of only one component of the wind stress. When the wind stress field has sufficient asymmetry about the equator, the situation becomes more complicated. Consider an antisymmetric oceanic mode forced by a symmetric meridional wind stress. Suppose we now add an antisymmetric zonal wind stress that is coherent with the meridional wind stress and which projects equally efficiently onto the wave’s meridional structure, but which has only th the power of the meridional wind stress. One might be tempted to ignore the small amount of power in the zonal wind stress and consider the forcing to be essentially due to a symmetric meridional wind stress. The amplitude of the zonal wind stress, however, is roughly that of the meridional wind stress, and the addition of what appears to be a relatively small amount of antisymmetric zonal wind stress could lower the power in the resonant response by over 50% or increase it by over 70%, depending on the phase relations between τys and τxa in (2).
When (or ) is nonnegligible, even rough quantitative predictions about the oceanic resonances require that we have sufficient knowledge of the meridional structure of the winds to carry out the projections in (2) [or (1)]. We must also know the phase relations between τx and τy. The limited meridional spacing and extent of the TAO/TRITON data used here prevents us from being able to make meaningful estimates of the projections in (1) and (2). For now, we note based on inspection of (1) and (2) that the ratios and can provide a rough measure of the importance of zonal and meridional stress in forcing the waves, assuming that the meridional structures of τx and τy project onto the structures of V(m)(y) and U(m)(y) with equal efficiency.2 With the above caveats in mind, some useful connections can still be made between the observed variability in the winds and the oceanic dynamic height (section 4e).
A final point to make regarding Fig. 2 is that, in contrast to the situation with P(m), the maximum amplitudes of both V(m) and U(m) are in the most poleward extrema. Winds at the most poleward TAO/TRITON latitudes (8°S and 8°N) are likely to be important in forcing meridional modes 2 and higher for baroclinic mode 1, and in forcing meridional modes 4 and higher for baroclinic mode 2, even though the pressure signals of these modes are much more tightly trapped to the equator.
We will draw comparisons between our results and the sea level, or sea surface height (SSH), analyses of previous authors, and accordingly we note here that the relation between SSH and dynamic height relative to 500 decibars depends on the baroclinic mode under consideration. Using the eigenfunction calculations described above (for estimating the gravity wave speed of the baroclinic modes), we computed the ratio of surface perturbation pressure (a proxy for SSH) to the difference in perturbation pressure between the surface and 500 m (a proxy for dynamic height relative to 500 decibars) for baroclinic modes 1 and 2 (Fig. 3, top two panels). This ratio gives a rough conversion factor for estimating the SSH amplitude that would be expected for a given amplitude of dynamic height relative to 500 decibars, if the baroclinic mode producing the dynamic height signal is known. The ratio for mode 2 is less than one and the ratio for mode 1 is greater than one, meaning that, in comparison to SSH, measurements of dynamic height relative to 500 decibars will enhance the second mode and suppress the first mode. This can be easily understood by noting that the vertical structure of pressure for both modes has a maximum amplitude at the surface, but the second baroclinic mode has antiphased extrema at the surface and 500 m, while the first baroclinic mode has pressure signals that are in phase at the two depths. Thus, while the surface-pressure (or sea level) signal of baroclinic mode 1 is larger than the dynamic pressure difference between the surface and 500 m (or dynamic height), the reverse is true for baroclinic mode 2. The bottom panel of Fig. 3 shows the ratio of the amplitude of the baroclinic mode-1 SSH oscillation to the amplitude of the mode-2 SSH oscillation, when the dynamic height amplitudes of the two modes are equal. Throughout the equatorial Pacific, equivalent dynamic height amplitudes for the two baroclinic modes would imply roughly twice the amplitude, or four times as much power, in the baroclinic mode-1 SSH signal as in baroclinic mode 2.
3. Data and methods
The spectral analysis presented here is based on 12 years (1997–2008 inclusive) of dynamic height estimates from the TAO/TRITON mooring array, which spans the equatorial Pacific waveguide with a nominal longitudinal spacing of 15° (e.g., McPhaden et al. 1998). All of the data used are from nine longitudes between 156°E and 95°W (indicated in Fig. 4). We use the daily-average estimates of the surface dynamic height relative to 500 decibars provided by the TAO Project Office of National Oceanic and Atmospheric Administration (NOAA)’s Pacific Marine Environmental Laboratory. The daily interval provides sufficient resolution for the period range that we consider, and the space-time coverage of the daily data retrieved by satellite telemetry is better than that of the higher-frequency, internally recorded data, which are not always recovered. Salinity is not typically measured on the TAO/TRITON moorings, so the TAO Project Office used a climatological temperature–salinity regression for each location and depth to estimate salinity. In cases where temperatures at 500-m depth were unavailable (e.g., instrument failure) and temperatures at 300-m depth were available, temperatures at 500 m were estimated by a regression analysis between available measurements from 300 and 500 m at that site. Further details of the dynamic height estimation procedure are available on the TAO Project Office web site (http://www.pmel.noaa.gov/tao/data_deliv/dyn.html); Busalacchi et al. (1994) also provide a discussion of the errors in dynamic height estimated by this approach.
We also use wind measurements from the TAO/TRITON surface buoys (also from 1997–2008) for interpreting the observed variability in dynamic height. The wind stress was estimated by application of the Large and Pond (1981) bulk formula to the measured winds, assuming neutral stability of the atmospheric boundary layer. Use of a stability-dependent drag coefficient requires coincident measurements of sea surface temperature, air temperature, and humidity. Estimates of wind stress made using the stability-dependent Coupled Ocean–Atmosphere Response Experiment (COARE) bulk flux algorithm (Fairall et al. 2003) are readily available for the TAO/TRITON array (http://www.pmel.noaa.gov/tao/oceansites), but the combined effect of the data gaps in the various required parameters leads to a marked reduction in the space-time coverage. We performed our analysis with both versions of the wind stress; the results were acceptably similar, but we judged it preferable to use all of the available wind data and an assumption of neutral stability in order to obtain better statistical stability in the spectral analysis.
b. General discussion of methods
The particular types of spectral analysis used here are modeled after an analysis of lower-frequency equatorial waves (periods >20 days) by Farrar (2008) and previous work (e.g., Wunsch and Gill 1976; Wheeler and Kiladis 1999; Roundy and Frank 2004). In particular, we will estimate the zonal wavenumber-frequency spectrum at each TAO/TRITON latitude and compute the latitudinally averaged spectrum by averaging spectral estimates from several latitudes (typically 5°S, 2°S, 0°, 2°N, and 5°N). The main reason for doing this is that different vertical and meridional modes have different meridional structures, so the zonal wavenumber–frequency spectrum at some particular latitude will be influenced in a complicated way by both the overall energy in the various modes and the meridional structure of those modes. For example, if a particular meridional mode had a zero crossing at the latitude under consideration, that mode would make no contribution to the spectrum at that latitude, even if the mode is more energetic than all the others. We present the latitudinally averaged spectrum in an attempt to provide a summary description of several modes at once. Given that the 2°–3° meridional spacing of the TAO/TRITON moorings cannot adequately resolve the meridional modal structures of baroclinic equatorial waves (especially higher meridional and vertical modes), our latitudinally averaged spectra are an imperfect representation of the true latitudinal average. Averaging over a number of TAO/TRITON latitudes is nonetheless beneficial for the present purposes, because our interest is not in the true latitudinal average but instead in identifying the modes that are present and understanding how the amplitude of a given mode is distributed in wavenumber and frequency. The principal risk associated with the lack of meridional resolution in this analysis is that the relative contribution of some modes to the latitudinally averaged dynamic height variance will be over or underestimated to the extent the TAO/TRITON latitudes tend to coincide with extrema or nodes of those modes.
We will also use a simple filtering procedure to decompose the variability into motions that are symmetric or antisymmetric about the equator (e.g., Yanai and Murakami 1970; Zangvil and Yanai 1980; Tsai et al. 1992; Wheeler and Kiladis 1999; Farrar 2008). Though our implementation of this procedure is slightly different than previous ones we have seen (next subsection), the basic idea of the filtering procedure is simply to add (or subtract) signals from opposite latitudes about the equator to cancel meridionally antisymmetric (or symmetric) contributions. Examination of the spectra of these symmetry-filtered fields, together with theoretical expectations about the symmetry properties of the various equatorial waves, can be useful for understanding which equatorial wave mode is responsible for an observed spectral peak, even if the observed modes do not have the perfectly symmetric or antisymmetric meridional structures predicted by classical equatorial wave theory (Farrar 2008).
Spectral analysis is clearly not the most direct way of gleaning information on some other properties of equatorial waves, such as their variability in time and space. There is surely spatial and temporal variability in the amplitude, propagation, and meridional structure of particular wave modes. For example, Luther (1980) showed that the amplitude of 3–5-day sea level oscillations was significantly lower at the Galapagos station (east Pacific) and west of the Marshall Islands than it was in the central Pacific. Seasonal modulation of wave amplitude seems likely, too. The results for a given spectral band represent a composite wave, averaged over time and longitude. As such, the spectral results here may be viewed as a description of the “climatology of equatorial waves” in the Pacific (using the terminology of Roundy and Frank 2004). The variability of equatorial waves in time and space is of interest, but it will need to be addressed in a different study.
The in situ data from the TAO/TRITON array, going back to the 1980s, make up a remarkable and unprecedented dataset, but there are numerous data gaps owing to the challenges of long-term deployments of surface moorings at sea (e.g., instrument failure, mooring failure, and fishing vandalism). The various data gaps in time and depth lead to a large number of prolonged (>10 days) gaps in estimates of dynamic height (which involves a vertical integral over depth and thus requires data from sensors at each depth). The time–longitude dynamic height data coverage is shown in Fig. 4. Our choice of a particular method for estimating the zonal-wavenumber/frequency spectrum was driven by a desire to use as much of the data as possible, despite the data gaps. The time period for this analysis (1997–2008, inclusive) was chosen because it is a time of relatively few data gaps. We will first discuss how we estimated the Fourier coefficients from longitude-time sections of data on each latitude, before discussing symmetry filtering and estimation of the zonal wavenumber–frequency spectra.
If data were available with uniform spacing in time and longitude, the most straightforward way of estimating the zonal-wavenumber/frequency spectrum would be to use a two-dimensional fast Fourier transform (FFT) and standard spectral methods (e.g., Farrar 2008, 2011). The dynamic height sampling is not uniform in time or longitude, primarily because of missing data, but also because the longitudinal spacing of the TAO/TRITON moorings is not regular in the Western Pacific. If only “small” data gaps (i.e., gaps shorter than the periods of interest) were present, it would be sensible to interpolate across the gaps and proceed with an FFT-based spectral estimate. We attempted this approach, and obtained results consistent with the ones presented below, but even with very aggressive interpolation (i.e., interpolating gaps longer than one month and interpolating across 30° gaps in longitude), we were only able to obtain about three, 180-day longitude–time sections for analysis at each latitude (on average) from more than 20 years of data.
Rather than discard so much hard-won data, we opted to estimate the Fourier coefficients by fitting sinusoids to the data in time and longitude using a variant of least squares. In other words, we suppose that each data point, hn, reflects contributions from M zonally and temporally periodic waves (and noise),
where ω and k are the angular frequency and zonal wavenumber, x and t are longitude and time, and am and bm are Fourier coefficients. Each ω–k pair is assigned an index m, and each longitude–time pair (i.e., each data point) is assigned an index n. The fit is performed independently for each latitude.
We could equivalently write,
where h is an N × 1 vector of the dynamic height observations from each longitude and time, n is an a posteriori estimate of the noise in the observations (i.e., n is the model–data misfit), x is a 2M × 1 vector of the Fourier coefficients, and is an N × 2M matrix containing columns of sin(kmxn + ωmtn) and cos(kmxn + ωmtn) evaluated at each m (i.e., k–ω pair) and n (i.e., the longitude-time of each sample). The form of and x is given explicitly by Wunsch (1989). The number of Fourier harmonics we have chosen to fit (described below) is less than the number of data points, making the inversion for x formally overdetermined. There are many viable means by which to estimate the Fourier coefficients– one could choose from methods like least squares, maximum entropy, Gauss-Markov, singular value decomposition, and Backus–Gilbert estimates, among many others.
After experimentation with several methods for estimating the Fourier coefficients in Eq. (4), we decided to use the tapered least squares approach (e.g., Wunsch 1989, 1996). We made this choice because the results were essentially the same regardless of the method used, and the tapered least squares approach is a relatively transparent and computationally efficient method that allows for noise in the data. In its general form, the tapered least squares approach is also known as ridge regression (e.g., Draper and Smith 1998, pp. 387–400) or as Tikhonov regularization (e.g., Aster et al. 2005, pp. 89–118), and with the particular choice of tapering parameter used here (based on an expected signal-to-noise ratio), the estimate could also be considered a Gauss-Markov estimate (and hence an optimal one, in the sense that it would have the smallest possible mean-square error; e.g., Wunsch 1996), if the a priori statistics of the Fourier coefficients and the noise met certain conditions.3 The tapered least squares solution to Eq. (4) is
which minimizes the sum of squares,
Minimization of the cost function J simultaneously minimizes the solution variance (proportional to xTx, the sum of squares of the Fourier coefficients) and the mean square model–data misfit (proportional to nTn), with a relative weighting, or tapering parameter, given by . In the general implementation of tapered least squares, the tapering parameter could be assigned any value, but, in viewing the tapered least squares method as a special case of the Gauss–Markov theorem (Wunsch 1989, 1996), it is clear that the tapering parameter can be considered to be an inverse signal-to-noise ratio (SNR), expressing the relative variance of the measurement noise and the squared Fourier coefficients. We thus chose to be the expected noise variance, and based on an expectation that Parseval’s relation should hold approximately, we chose to be the sample variance on each latitude divided by the number of Fourier coefficients estimated.4 We took the standard deviation of the noise (σn) in the daily-average dynamic height to be 0.5 cm. The results are almost completely insensitive to the specific choice of , because the SNR is high.
Prior to estimation of the Fourier coefficients for each latitude using Eq. (5), we performed some data preprocessing. We first removed the time-mean dynamic height at each location, which necessarily makes the longitude–time mean at each latitude zero and also removes the time-mean zonal structure (the most substantial trend in the data). It was at this point that we computed the data variance on each longitude-time section for specification of the SNR in the tapered least squares problem formulation. We then tapered the longitude-time sections to zero using Hann (i.e., squared cosine) windows to reduce spectral sidelobes (e.g., Thompson 1971; Scargle 1982). Because we use only nine TAO/TRITON longitudes, tapering the data records to zero in longitude would eliminate about 22% (i.e., ) of the data, which is clearly unacceptable, so we instead added fictitious, zero data at longitudes 15° east and west of the array (i.e., we zero padded with one longitude point on each side). After estimating the Fourier coefficients on each latitude, we renormalized them so that the spectrum would have the same total variance as the dynamic height data did before the taper was applied.
A key choice to be made in formulating the problem posed in Eq. (4) is the choice of frequencies and zonal wavenumbers at which estimates of the Fourier coefficients are sought. In principle, one could attempt to estimate the Fourier coefficients for any zonal wavenumber and frequency, but the record length and sample spacing (and data gaps) limit the spectral resolution that can be attained. We invested substantial effort in assessing the resolution of the time-longitude sampling and the sensitivity of our results to the frequencies and zonal wavenumbers used in and settled on the following, commonsense choice of uniformly spaced frequencies and zonal wavenumbers. With the 1-day average dynamic height data used here, the shortest resolvable period is two days, so the maximum frequency included in the fit is ωmax = 2π/2 days. With most of the TAO/TRITON array having a zonal spacing of 15° of longitude (Fig. 4), the shortest wavelength that we will attempt to fit is 30° (i.e., kmax = 2π/30°). The remaining choices concern the fundamental harmonics that determine the spacing in wavenumber and frequency. The longest zonal wavelength that we could reasonably expect to resolve is equal to the zonal extent of the array (109°), but we actually chose to use a slightly longer fundamental wavelength (135°, or Δk = 2π/135°) as a compromise that minimizes spectral sidelobes because it is a multiple of the most common zonal spacing (15°) and nearly fits the zero-padded zonal domain (139°) (serving to make the discrete sinusoids nearly orthogonal). The consequence of using a fundamental wavelength longer than the data array is a slight decrease of statistical independence between adjacent wavenumber bands, as also happens when the spectrum is computed for a uniformly spaced, zero-padded or tapered record using an FFT. The highest frequency resolution that we could reasonably expect is one cycle per record length (12 years), but, to limit the size of the matrix inversion in Eq. (5), we chose a fundamental frequency of -years (i.e., Δω = 2π/6-years). Using twice as large a fundamental frequency should be similar to computing a two-band average of the spectrum estimated with the smaller fundamental frequency. The frequencies and wavenumbers used in the fit thus span wavenumbers of ±kmax with a spacing of Δk and frequencies of 0 to ωmax with a spacing of Δω. With this choice of frequencies and wavenumbers, each latitude has about three times more data points than spectral bands, making the fit formally overdetermined. After computing the spectrum, we averaged over 13 frequency bands to further increase the degrees of freedom and statistical reliability of the estimate.
Given an estimate of the Fourier coefficients, the power spectral density can be expressed
where the angle brackets indicate averaging over latitude and/or frequency, Δk and Δω are the wavenumber and frequency bandwidth, and the factor of two in the denominator yields the mean square of the sine and cosine components ( and ). For the averaging indicated by angle brackets, most results shown are averaged over blocks of 13 frequency bands and over five latitudes (5°S, 2°S, 0°, 2°N, and 5°N); the main exception is when spectral density is compared at different latitudes, in which case no latitudinal averaging is done. The factors of 2π are included so that the power spectral density is given in units of dynamic-centimeters squared per cycle-per-day per cycle-per-degree (not in radians). We have not included data from 8°S and 8°N in the averages because the spectral power in dynamic height is expected to be relatively low at these latitudes for all modes of interest (see Section 2 and Fig. 2). Including data from these latitudes would only serve to lower our signal-to-noise ratio.
One approach for examining motions that are symmetric or antisymmetric about the equator (e.g., odd and even numbered equatorial modes) is to filter the data for cross-equatorial symmetry by adding or subtracting data from opposite latitudes about the equator. Doing this in the time–space domain with the TAO/TRITON data is not desirable, as it would nearly double the number of missing data, but, given estimates of the Fourier coefficients at opposite latitudes, the spectrum of symmetry filtered dynamic height can be expressed
where use of the minus sign, for example, yields the spectrum of equatorially antisymmetric variability formed from latitudes y° from the equator. Normalization by a factor of four, rather than by the factor of two used in Eq. (7), has no special significance, but it ensures that the average of the symmetric and antisymmetric spectra equals the average of the spectra computed from the two latitudes separately. Note that this convention has the effect of doubling the spectral density when the signal is either purely symmetric or purely antisymmetric about the equator. For example, when the variability is purely symmetric, . evaluated only on the equator is always twice the value of Ψm evaluated on the equator, because the antisymmetric part is zero there by definition.
A formal estimate of the error in the estimated Fourier coefficients is a natural part of the tapered least squares solution [e.g., Wunsch 1989, his Eq. (11b)]. Our primary interest here is not in the error of the Fourier coefficients themselves, but in the error of their mean square, i.e., in the power spectral density [Eqs. (7) or (8)]. We can estimate the random spectral error by noting, in analogy with the usual error estimates for Fourier spectra, that if the input data (e.g., dynamic height) were a normally distributed random variable, the linear estimate of the Fourier coefficients would also be random and normally distributed, so their mean-square value would follow a chi-square distribution, with percentage confidence intervals for the spectral density following directly (e.g., Bendat and Piersol 2000, p. 98). We take a reasonable, but somewhat optimistic, estimate of the number of degrees of freedom for each estimated Fourier coefficient to be the ratio of the number of data points to the number of Fourier coefficients estimated. Because of the different data gaps at the different latitudes, the estimated number of degrees of freedom per Fourier coefficient varies between latitudes; it was 1.5 ± 0.05 for all latitudes on 5°S to 5°N, while estimates at 8°S and 8°N were lower (1.0 and 1.4). Taking spectral errors from different latitudes and frequency bands to be independent, the total number of degrees of freedom for the estimate in a single band [i.e., Eq. (7)] will be increased in proportion to the number of latitudes times the number of bands averaged. For example, for the spectrum shown in Fig. 5 (section 4), 13 frequency bands and five latitudes were averaged, so the spectrum was estimated to have 195 degrees of freedom (≈1.5 × 2 × 13 × 5). This estimate of the degrees of freedom is clearly optimistic, but we have found the resulting confidence intervals to be a useful measure of the statistical uncertainty, and we note that the stated 95% confidence intervals would be about the same size as 90% confidence intervals if the estimated number of degrees of freedom were reduced by 30%. Note that adjacent frequency and wavenumber bands are not independent because of the tapering and zero padding, but spectral bands separated by a few points in wavenumber or frequency can be considered independent.
a. Dynamic height spectra
The power spectral density of dynamic height, averaged over the 5°S–5°N part of the TAO/TRITON array, has spectral peaks that coincide with the dispersion curves of several distinct equatorial modes (Fig. 5). Recall our abbreviated notation for meridional and baroclinic mode numbers (section 2): baroclinic mode 2, meridional mode 1, for example, is referred to as “bc2m1.” We see spectral peaks on the dispersion curves for bc1m0, bc2m0, bc2m1, bc1m2, bc1m3, and near 5-day periods, where the dispersion curves of bc1m1 and bc2m2 overlap. There is also a suggestion of elevated spectral power near higher modes (e.g., bc1m4).
Taking advantage of the symmetry properties of linear waves, we can gain further confidence in the mode identification by calculating separately the spectra of dynamic height fluctuations that are symmetric or antisymmetric about the equator (Fig. 6). The spectral peaks in these symmetry-filtered spectra agree with expectations from linear theory. In the antisymmetric spectrum (Fig. 6, right panel), the peaks for the antisymmetric modes bc1m0, bc2m0 and bc1m2 are enhanced relative to the background, while the peak for the symmetric mode bc2m1 has vanished. Spectral levels in the antisymmetric spectrum remain elevated near 5-day periods, suggesting that the antisymmetric mode bc2m2 contributes to the peak seen in the unfiltered spectrum near 5-day periods. In the symmetric spectrum (Fig. 6, left panel), the 7-day peak associated with the symmetric mode bc2m1 (near 0.15 cpd) is enhanced relative to the background, and we see weak peaks near 3.5 and 4.5 days, apparently associated with the symmetric modes bc1m3 and bc2m3. We also see elevated spectral levels in the symmetric spectrum near 5 days (0.2 cpd), suggesting that the dynamic height variance at 5-day periods may contain roughly equal contributions from the symmetric mode bc1m1 and the antisymmetric mode bc2m2.
The symmetric–antisymmetric identification is facilitated by examination of the logarithm of the ratio of the symmetric spectrum to the antisymmetric spectrum (Fig. 7). This quantity expresses the relative magnitude of symmetric and antisymmetric variability; for example, a value of +0.3 means that the symmetric power exceeds the antisymmetric power by a factor of 100.3 ≈ 2. One virtue of presenting the symmetry-filtered spectra this way is that, because the logarithm of the ratio is identical to the difference of the logarithms, displaying a confidence interval is relatively straightforward, so that one can easily assess whether there is a statistically significant difference between symmetric and antisymmetric energy.5 There is a statistically significant preference for meridionally symmetric dynamic height fluctuations near the theoretically symmetric Kelvin waves (bc1 and bc2), bc2m1, bc2m3, and bc1m3 modes (Fig. 7). There is a statistically significant preference for antisymmetric dynamic height fluctuations near the theoretically antisymmetric MRG (bc1m0 and bc2m0) and the bc1m2 mode. There is also a suggestion of an antisymmetric preference near the bc2m4 and bc1m4 dispersion curves at 3.8 and 3 days, although the preference is not significant at 95% confidence.
The strong 5-day peak, seen in both the symmetric and antisymmetric spectra, is not evident in Fig. 7 because of the nearly equal amounts of symmetric and antisymmetric spectral energy density. There is a preference for meridionally symmetric variability on the bc1m1 dispersion curve near the larger positive wavenumbers, where this dispersion curve separates from that of antisymmetric mode bc2m2. The interpretation that the 5-day peak contains two distinct modes is further supported by an examination of the spectral power as a function of latitude, which we carry out in the next subsection.
b. Meridional structures
The relatively sparse latitudinal spacing of the TAO/TRITON array does not lend itself to a clear identification of the meridional structure of high meridional modes, but it is worth examining the latitudinal dependence of the 5-day peak, which lies on the dispersion curves of two low meridional modes of opposite symmetry. The total, symmetric, and antisymmetric spectra were calculated for each of the TAO/TRITON latitudes (8°S, 5°S, 2°S, 0°, 2°N, 5°N and 8°N). The spectra were averaged over the wavenumber–frequency box shown in the right panel of Fig. 8 and plotted against latitude in the three left-hand panels. Superimposed on the top two panels is the theoretical meridional structure of the squared pressure fluctuation (“P2(y)”) for a symmetric bc1m1 wave at k = 0, the central wavenumber of the averaging box. The theoretical structure is scaled, somewhat arbitrarily, to have a maximum value 1.2 times that of the largest spectral density estimate in this wavenumber–frequency band. The squared pressure structure of the antisymmetric mode bc2m2 is superimposed on the bottom-left panel, similarly scaled. Although the mooring array cannot adequately resolve the meridional structure of the waves, the meridional variations in spectral levels of the symmetric and antisymmetric signals are consistent with the theoretical structures of the symmetric bc1m1 and antisymmetric bc2m2 modes, giving support to the interpretation that both modes are present near 5-day periods.
The 5-day peak was identified in tide gauge records by Wunsch and Gill (1976) and Luther (1980) as the symmetric mode bc1m1, based in part on examination of the sea level spectrum as a function of latitude, which showed good agreement with the theoretical structure of bc1m1. The meridional structure of the total dynamic height spectrum in this band (upper-left panel of Fig. 8) does not resemble the structure of bc1m1, but this reinforces our interpretation that the antisymmetric energy in dynamic height is likely due to bc2m2. As noted in section 2, when dynamic height relative to 500 decibars is converted to SSH, any baroclinic mode-2 energy should be attenuated by roughly a factor of 4 relative to baroclinic mode-1 energy. If the peak in the antisymmetric spectrum is due to baroclinic mode 2, then the meridional variation of the total SSH spectrum should look more like the middle left-hand panel, in agreement with the SSH analyses of previous authors.
Figure 9 shows the meridional structures of the dynamic height spectrum (without symmetry filtering) at the wavenumbers and frequencies where the dispersion curves of the lowest five meridional modes (MRG and four lowest inertia–gravity modes) cross the k = 0 axis, for both baroclinic modes 1 and 2. (Each spectral estimate shown is an average over the 13 frequency bands and two wavenumber bands nearest to the wavenumber and frequency where the dispersion curve crosses k = 0.) Figure 10 shows the structure of the symmetric or antisymmetric spectra at the same wavenumbers and frequencies, with the symmetric or antisymmetric spectra being chosen for display depending on the theoretical symmetry of each mode. The higher meridional modes cannot be resolved by the TAO/TRITON array, and measurements at individual latitudes become more sensitive to slight shifts in the meridional structure. Nevertheless, the agreement between the observed variability and the theoretical wave structures is reasonably good in all cases except bc1m3, near 3.5 days. The agreement with the theoretically symmetric meridional structure for this mode is no better for the symmetric spectrum (Fig. 10) than for the total spectrum (Fig. 9). We speculate that this band may contain contributions from another theoretically symmetric mode with a nearby dispersion curve (i.e., bc2m5), making the symmetry filtering ineffective for isolating the bc1m3 mode.
c. Wavenumber content of the oceanic variability
In addition to the identification of second baroclinic mode inertia–gravity wave variability and the consistency checks with linear theory via symmetry filtering, the other major novelty of Figs. 5–7 is their depiction of the wavenumber content of the oceanic spectral peaks. The wavenumber limits of the plots are the nominal Nyquist wavenumbers of the TAO/TRITON array (±1 cycle per 30° longitude), and the free-wave dispersion curves are near their minimum frequency throughout this wavenumber range. Nevertheless, it is possible to make some important observations, particularly with respect to meridional modes 1 and 2, for which the wavenumber at the frequency minimum (locus of vanishing group velocity) is the most distinctly removed from k = 0.
For meridional mode 1, the spectral power in both baroclinic modes (1 and 2) remains high across most of the resolvable wavenumber band (Fig. 6, left panel). For meridional mode 2, the spectral power in both baroclinic modes is more confined to small wavenumbers (Fig. 6, right panel). There is a weak spectral peak near the frequency minimum of bc1m1 (e.g., Fig. 6, left panel), but it is only about 10% higher than the power level across most of the rest of the wavenumber band, which is not a statistically significant difference. For bc1m2 and bc2m2 (Fig. 6, right panel), the spectral peaks are at the smallest resolvable negative wavenumber, which is less negative than the wavenumber associated with vanishing group velocity. However, again the difference between the spectral power at these peaks and that at the vanishing-group-velocity wavenumbers is not significant with 95% confidence. For bc2m1 (Fig. 6, left panel), the spectral peak is actually at a positive wavenumber, noticeably removed from the vanishing-group-velocity wavenumber. The power at this peak is higher than at the vanishing-group-velocity wavenumber with 95% confidence. The presence of significant spectral power and even spectral peaks at locations on the dispersion curves associated with nonnegligible group velocity suggests that any oceanic preference for enhancing resonant energy at the locus of vanishing group velocity is not the dominant mechanism for setting the wavenumber of the spectral peaks. We have left for last the obvious presence of enhanced spectral energy density on the MRG dispersion curves, which have no loci of vanishing group velocity. The MRG group velocity at vanishing wavenumber is several times that of the inertia–gravity waves, yet we see evidence of strong variability in both baroclinic modes near k = 0. The remaining question is whether the locations of the oceanic peaks can be explained by the wavenumber content of the wind forcing.
d. Wind stress spectra
Keeping in mind the expected correspondence between wind components and oceanic modes seen in (1) and (2), we show in Fig. 11 the wavenumber–frequency spectra for the symmetric and antisymmetric parts of the zonal wind stress (τxs and τxa) and the symmetric and antisymmetric parts of the meridional wind stress (τys and τya). In light of the expected importance of higher-latitude wind forcing (Section 2), these spectra are averaged over 8°S–8°N. Also shown are dispersion curves for atmospheric equatorial waves with a 30-m equivalent depth (magenta lines), a value near the center of the equivalent-depth range identified by Wheeler and Kiladis (1999) in a study of convectively coupled atmospheric equatorial waves using outgoing (top-of-atmosphere) long-wave radiation measurements.
Both τx and τy show a ridge of spectral power centered near k = −1 cycle per 75°–90° longitude (planetary wavenumber 4–5). The most energetic part of the meridional wind stress spectrum is found where the atmospheric MRG dispersion curve crosses this ridge. The zonal wind stress has an additional ridge of spectral power at positive wavenumbers that follows the dispersion curve of the atmospheric Kelvin wave, and we will refer to this ridge as such. Although there is some frequency dependence of the spectral power along all of the ridges, the prominent peaks found in the oceanic spectra are not evident in the atmospheric spectra (as noted by Wunsch and Gill 1976, and Luther 1980). The spectral power in the ridges falls off with increasing frequency at periods shorter than about 4 days.
The atmosphere, with no clearly defined upper boundary, cannot support vertically standing modes as the ocean can, so the comparatively diffuse spectral ridges are to be expected (e.g., Philander 1978). Still, the apparent utility of the 30-m equivalent depth for describing the distribution of energy in surface wind stress suggests that the surface wind signals are related to the heavily studied atmospheric convectively coupled equatorial waves (e.g., Wheeler and Kiladis 1999). The τxs ridge paralleling the Kelvin wave dispersion curve is even more pronounced when the atmospheric spectra are averaged over the more equatorially confined range of 5°S-5°N (not shown), and the ridge exists only in the symmetric spectrum of zonal wind stress; both facts support the interpretation of this variability as Kelvin waves.
Both τx and τy are dominantly symmetric over the wavenumber-frequency range of interest (cf. left and right panels of Fig. 11). We saw in (1)–(2), however, that for the purpose of assessing the importance of symmetry in the forcing, a better diagnostic than the symmetry of an individual component is the ratio or (Fig. 12). Of the places where the atmospheric spectral ridges cross dispersion curves of oceanic modes that could be excited, only at the intersection of the atmospheric Kelvin wave with the bc2m1 oceanic curve (near 0.15 cpd) does the ratio of spectral power in the symmetric part of the wind stress (|τxs|2) to the power in the corresponding antisymmetric part (|τya|2) exceed 10 (just barely). As discussed in section 2, the power in the response can be sensitive to even this small amount of antisymmetric forcing. At other locations for the inertia–gravity modes, the ratios are about 6 or less, meaning that the amplitudes of the symmetric components of the wind stress are not more than about 2.5 times the amplitudes of the relevant antisymmetric components. Throughout the MRG wavenumber-frequency range, the amplitudes of τys and τxa are nearly equal. With ratios this low, the antisymmetric stress component could make a significant contribution to, or even dominate, the net forcing, depending on how the symmetric and antisymmetric wind stress components project onto the equatorial wave structures [(1)–(2)].
A full understanding of the relative importance of antisymmetric and symmetric forcing requires evaluation of the projections in (1) and (2), which we cannot presently do given the marginal meridional resolution of the TAO/TRITON array. We can make a rough comparison between forcing and response based on the simplifying assumption of total symmetry of the forcing, which we do in the next section. In doing so, however, we must keep in mind that we do not really have all the information we need, and it is reasonable to expect discrepancies.
e. Forcing-response comparisons based on symmetric forcing
Figure 13 shows contours of normalized symmetric wind stress spectra superimposed on the oceanic spectra that they would be expected to excite: zonal stress with the symmetric dynamic height spectrum (left panel) and meridional stress with the antisymmetric dynamic height spectrum (right panel). Within any frequency band containing elevated oceanic spectral power, we are interested in the relation between the wavenumber content of the oceanic spectrum and that of the appropriate symmetric wind stress component. For display purposes then, we divide the spectral density of wind stress at each frequency by the average spectral density at that frequency and refer to the result as the normalized wind stress spectrum. This allows for fewer contours and less clutter in the plots, while capturing the basic pattern of wavenumber dependence in the forcing at any given frequency. The normalized wind stress spectra are contoured at levels of 0.8, 1.2, and 1.8, so these contours indicate the wavenumber where the spectral density is 80%, 120%, and 180% of the average spectral density for each frequency band.
The qualitative agreements seen in Fig. 13 are reasonable given our uncertainty about the meridional projections of wind stress onto wave velocities and the degree of asymmetry in the wind. The double ridge in τxs spans most of the resolvable wavenumber band at periods longer than about 5 days (0.2 cpd), and the spectral peaks in the two oceanic modes expected to be excited by these winds, bc2m1 and bc1m1, also span most of the wavenumber band (left panel, Fig. 13). In contrast, the peaks for the antisymmetric oceanic modes bc2m2 and bc1m2 are somewhat more narrowly confined, like the single τys ridge. These latter oceanic peaks appear displaced toward positive wavenumbers relative to the τys ridge, but in the case of bc1m2, the peak in |τys|2 is also displaced from the center of the ridge (Fig. 11, bottom left panel), and in both cases τxa is large enough to influence the response. With the exception of the peak on bc1m0 near the most negative wavenumbers (12 days period, 0.085 cpd), the bc1m0 and bc2m0 MRG peaks fall within the wavenumber range of the most energetic part of the τys ridge.
Where the atmospheric Kelvin wave dispersion curve crosses the bc2m1 dispersion curve near 0.15 cpd, the ratio exceeds 10 and the assumption of purely symmetric forcing is at its best. Our estimates of the peak locations for both the oceanic and atmospheric spectra are at the same wavenumber and frequency near this crossing, although the confidence intervals are large enough that either peak could actually lie in nearby wavenumber or frequency bands. The strong peak in the positive-wavenumber part of the oceanic spectrum spans wavenumbers for which the group velocity of bc2m1 ranges from 17%–35% of the Kelvin wave speed. This is perhaps the strongest evidence that the oceanic peaks are determined more by the atmospheric spectrum than by the condition of vanishing group velocity for the oceanic modes. The correspondence between the τys ridge and the peaks in the MRG waves, with group velocities approaching 50% of the corresponding Kelvin-wave speeds, further supports this view.
The weak peak in the symmetric oceanic spectrum near the bc1m3 dispersion curve (about 0.3 cpd) appears to be centered on a local minimum of the τxs spectrum (Fig. 13, left panel). This appearance in Fig. 13, however, is due in part to the atmospheric Kelvin wave ridge that has most of its power in the near-equatorial latitudes and hence has a minimal effect on the bc1m3 mode. An inspection of the U(3) structure in the right panel of Fig. 2 shows that, of the wind measurements available from the TAO/TRITON array, bc1m3 should respond most strongly to symmetric zonal winds at 8°S and 8°N, slightly less strongly to winds at 5°S and 5°N, and hardly at all to the more equatorially confined winds. A spectrum of τxs estimated using only the wind stress at 8°S and 8°N (not shown) does not exhibit a minimum at the location of the bc1m3 peak.
The largest discrepancy between the comparisons of Fig. 13 and our expectations for the correspondence of wind stress and oceanic resonances is the strong peak on the bc1m0 dispersion curve near the negative Nyquist wavenumber, where the power in τys is falling off. It is possible that other forcing mechanisms contribute to the oceanic spectrum, and this could be one example. Nevertheless, the qualitative comparison made under the simplifying (but only marginally valid) assumption of purely symmetric forcing provides a reasonably satisfying degree of success in matching the wavenumber content of the oceanic peaks to the wavenumber content of the wind stress. More work is needed to better define the oceanic variance driven by fluctuating winds, but a more detailed assessment will require wind measurements with meridional resolution sufficient to estimate the projections in (1)–(2).
A long data record from the TAO/TRITON mooring array was analyzed to allow a more detailed assessment of the spectra of dynamic height relative to 500 decibars and wind stress in the equatorial Pacific at periods of 3–15 days. In the dynamic height spectra, elevated spectral density is found on or near the dispersion curves predicted by linear theory for the lowest 3–4 meridional modes of baroclinic modes 1 and 2. The baroclinic mode 2 variability had not been identified in studies of SSH variability, but it shows up clearly in the TAO/TRITON dynamic height data. The reference depth (500 m) is below the first zero crossing and near the second extremum in the vertical structure of perturbation pressure for baroclinic mode 2, so the surface dynamic height signal of this mode is amplified relative to the surface-pressure and SSH signals. The meridional resolution of the TAO/TRITON array is limited, but the meridional dependence of dynamic height variability within the spectral peaks generally agrees well with the predictions of linear theory. Filtering the variability with respect to its equatorial symmetry also produces good agreement with linear theory. The spectral peaks associated with a particular symmetry (symmetric or antisymmetric) are enhanced or eliminated in accord with expectations from linear equatorial-wave theory. The symmetry filtering proved particularly helpful in examination of the spectral peak near 5-day periods, which had been identified in frequency spectra of SSH variability as a baroclinic mode 1, meridional mode 1 wave (Wunsch and Gill 1976; Luther 1980). The symmetry-filtered spectra suggest that the dynamic height signal at this period contains roughly equal contributions from the symmetric baroclinic mode 1, meridional mode 1, and the antisymmetric baroclinic mode 2, meridional mode 2. As suggested by previous authors (e.g., Wunsch and Gill 1976; Luther 1980), we find that both the zonal and meridional winds within 8° of the equator are dominantly symmetric about the equator and concentrated in low-wavenumber bands in the 3–15-day period range.
In general, we find spectral peaks in dynamic height where the spectral ridges in the wind stress cross the dispersion curves of the free oceanic modes. Some understanding of the connection between the patterns of forcing and response in wavenumber-frequency space is gained via the expedient of considering the winds to be purely symmetric about the equator. There are discrepancies, and we show that despite the dominantly symmetric winds, there is almost nowhere in the wavenumber-frequency range of interest where the possible contributions of the asymmetries in the winds can be ignored. The only exception is where the atmospheric Kelvin wave ridge crosses the dispersion curve of the baroclinic mode 2, meridional mode 1 oceanic wave, and in this case the estimates of peak energy in the oceanic variability and the atmospheric Kelvin wave ridge are indeed found at the same wavenumber.
Other mechanisms could of course contribute to the equatorial variability at periods of 3–15 days and zonal wavelengths exceeding 30°, but the agreement between the patterns in the atmospheric and oceanic spectra, when considered in light of the limits of the information available to us from the TAO/TRITON array, suggests that direct wind forcing be examined in more detail before resorting to alternate explanations for the oceanic variability.
Our results do not support the expectation of finding elevated dynamic height or sea level variance at the wavenumber and frequency of vanishing group velocity, unless the wind stress happens to fluctuate strongly at those wavenumbers and frequencies. The results do support the prediction of Wunsch and Gill (1976) that the frequency spectra they observed were likely produced by atmospheric forcing confined to low-wavenumber ridges, narrow enough to limit the oceanic response to frequencies not too far removed from the modal minima. Our results also confirm the prediction of Wunsch and Gill that the atmospheric wavenumber ridges must also be wide enough for the MRG response to be spread across a wide band of frequencies, in order to account for the fact that they did not see distinct peaks from these waves in frequency spectra. It is remarkable that the qualitative nature of the wavenumber–frequency spectrum of the wind stress, complete with upper and lower bounds on the widths of the wavenumber ridges, was inferred from scattered tide gauge records.
A thorough analysis that includes the effects of asymmetries in the forcing will require measurements with meridional spacing and extent sufficient for estimating the meridional projection of the wind stress components onto the meridional structure of velocity for the mode in question, as shown in (1)–(2). We expect that satellite scatterometer and altimeter data will allow a more detailed investigation of the relationship between the wind field and the oceanic response, and one purpose of the present study is to provide a baseline of information at relatively high frequencies with which to compare analyses of satellite data. In any case, the spectra shown here present an interesting qualitative picture: the equatorial Pacific Ocean appears to be somewhat like a tightly strung guitar– details of the forcing aside, it rings at the frequencies (and zonal wavenumbers) of its free modes of oscillation.
We are grateful to Dudley Chelton for numerous helpful and stimulating discussions during the course of this research. We thank Mike McPhaden, Billy Kessler, and an anonymous reviewer for helpful comments on the manuscript. We thank Michael Schlax for helpful critical comments on an early version of the manuscript. We also benefitted from discussions with Carl Wunsch, Jay McCreary, Ken Brink, Breck Owens, Alexey Kaplan, Deepak Cherian, and Jake Gebbie. We thank the TAO Project Office of NOAA’s Pacific Marine Environmental Laboratory for providing the dynamic height and wind data, and we thank the TAO Project and TRITON Project for their sustained efforts to collect the data. This research was funded by NASA Grant NNX10AO93G and is a contribution of the Ocean Vector Winds Science Team.
As the zonal wavenumber goes to infinity, the inertia–gravity wave modes increasingly resemble pure, zonally propagating gravity waves, and the amplitude of V(m) becomes small compared to that of U(m). In this limit, the ratios and have no special significance for the excitation of the waves– the waves will only respond to τx, regardless of the strength of τy (Durland and Farrar 2012). For MRG, this is also true as k → ∞, but as k → −∞, the waves only respond to τy. Within the wavenumber range considered here, the relative amplitudes of V(m)(y) and U(m)(y) are comparable for both wave types, so the antisymmetric part of the wind stress cannot be neglected a priori.
These conditions are that the measurement noise be white and the best-guess a priori estimate of the Fourier coefficients be one that assigns them all equal value, with a sum of squares equal to the expected dynamic height variance. The latter would correspond to a white a priori spectrum, but with a very particular phase arrangement. The dynamic height spectrum is not expected to be white, so our estimate is not optimal, but we have chosen the tapered least squares approach in part out of a spirit of agnosticism, to avoid giving the impression that our expectations have influenced the result.
Parseval’s theorem states that, for uniformly spaced samples, the integrated spectral power equals the variance of the data. This result relies on the orthogonality of sinusoids with the usual Fourier frequencies on a uniform sampling grid, but sinusoids on an irregular grid are not in general orthogonal, so Parseval’s theorem does not generally hold for a spectrum estimated from gappy data. In addition, a key premise of the tapered least squares method is that there is noise in the data and that a solution that represents all of the variance in the data is thus undesirable. Nonetheless, given our reasonably good SNR, we view an acceptable solution as one that has about the same amount of variance in the spectrum as in the raw data.
The confidence interval displayed in Fig. 7 is for judging whether there is a statistically significant difference in symmetric power and antisymmetric power at a particular wavenumber and frequency, which is done by using the confidence interval to determine whether the value plotted is bounded away from zero (i.e., whether the ratio of the symmetric spectrum to the antisymmetric spectrum is significantly different from 100 = 1). The upper and lower confidence intervals are each about 2.5 contour intervals, so any values in the figure that are three contour intervals away from zero represent a statistically significant symmetry preference.