TAO Data Support the Existence of Large High Frequency Variations in Cross-Equatorial Overturning Circulation

Large amplitude oscillations in the meridional overturning circulation (MOC) have been found near the equator in all major ocean basins in the NEMO ocean general circulation model. With periods of 3–15 days and amplitudes of ∼ ±100 Sv in the Pacific, these oscillations have been shown to correspond to zonally integrated equatorially trapped waves forced by winds within 10°N/S of the equator. Observations of dynamic height from the Tropical Atmosphere Ocean (TAO) mooring array in the equatorial Pacific also exhibit spectral peaks consistent with the dispersion relation for equatorially trapped waves. Here, we revisit the TAO observations to confirm that the amplitude of the oscillations is consistent with the simulations, supporting the modeled large amplitude MOC oscillations. We also show that the zonal structure of the frequency spectrum in both observations and simulations is predicted by changes in the baroclinic wave speed with variation in stratification across the ocean basin.

A recent two-part study by Blaker et al. (2021) (hereafter B1) and Bell et al. (2021) (hereafter B2) showed that large amplitude oscillations of the MOC found in an integration of the NEMO (Nucleus for European Modeling of the Ocean) global ocean model could be attributed to equatorially trapped Yanai and Poincaré waves driven by winds within ±10° of the equator.Using idealized simulations based on a zonally averaged formulation of the linear wave theory with zonal boundary conditions of no normal flow, they accurately reproduced the MOC oscillations in the realistic model with a small number of low-mode waves, thus explaining the modeled MOC variability on 3-15 timescales.The 10 day Yanai and 4 day Poincaré modes were shown to give the largest contributions to the equatorial MOC oscillations.However, the question remains as to whether such large oscillations occur in reality.
In this study, we compare output from the global ocean simulations of B1, 2 with the TAO mooring record to verify that zonally integrated fluctuations in dynamic height are consistent with those inferred from the TAO mooring observation data.An early indication of this is the good correspondence between high frequency, zonally averaged dynamic heights from NEMO and TAO data at the equator, plotted for the same time period in Figure 1a.Although the moorings are relatively sparse in latitude and longitude, FD12 showed that it is possible to estimate the temporal and zonal spectrum of dynamic height using least squares methods, finding that most of the power was concentrated at low zonal wavenumbers.Motivated by this and the findings of B1, 2 that the zonally averaged effect of these waves could have important consequences for variability of the MOC, we focus on zero zonal wavenumber modes and use similar methods to FD12 to fit the high frequency dynamic height data to a Fourier basis in time and a basis of theoretically derived meridional modes in latitude.We first compare the temporal/ meridional spectra of zonally averaged dynamic height in the model and observations, and show that they are largely consistent with each other and with linear theory.We then explain the zonal structure of the dynamic height power spectrum in terms of the east-west variation in stratification.

Equatorial Modes
Equatorially trapped waves can be modeled as linear perturbations upon a base state with zero background current and a background stratification N 2 (z).The solution for unforced waves on a beta plane can be found from the hydrostatic, incompressible, linear, unforced and inviscid momentum and buoyancy equations (Bell et al., 2021;Durland & Farrar, 2012;Wunsch & Gill, 1976).The meridional velocity is given by: where, p m (z) are a set of vertical modes and ϕ m,n (y) is the nth meridional mode for the mth baroclinic mode.Full expressions are given in Equations A1-A3.The frequency ω m,n and zonal wavenumber k satisfy the dispersion relation: where, c m is the mth baroclinic wave speed (see Equation A1), β is the Rossby parameter, and V m,n (k) determines the zonal wavenumber content of each vertical/meridional mode combination.In an unforced and undamped scenario, V m,n (k) would be set by the initial conditions.However, in reality and in the model the waves are both forced by surface winds and damped by turbulent processes (parameterized in the model).Durland and Farrar (2012) derived the relationship between the wind forcing and the resulting wave field, and B1, 2 showed using a global ocean model and idealized simulations that wind forcing between ±10°N is essential to the generation of the resulting MOC fluctuations.The effect of the wind forcing is to resonantly excite particular components of the free wave solution, dependent on the spatial structure and frequency content of the winds.We therefore expect to find waves corresponding to the unforced solutions (1) in the forced ocean.
Figure 1f shows the dispersion relation (2) for n = 0, 1, 2, 3, 4, 5, m = 1, 2. The notation used is the same as that of FD12: 'BC m' refers to the mth baroclinic mode, and 'M n' refers to the nth meridional mode.The curves in the lower left are Rossby waves, and the dispersion relations ω = c m k for Kelvin waves are not shown, since we are interested in long zonal wavelengths and high frequencies.The n ≥ 1 waves are Poincaré waves, and the n = 0 waves are Yanai waves.FD12 found spectral peaks in the TAO dynamic height record of long zonal wavelength (k ≃ 0) waves agreeing well with these dispersion relations.
We estimate the meridional spectrum of the waves by fitting to theoretical meridional modes, rather than a standard Fourier basis.Since dynamic height is proportional to pressure, instead of using the meridional velocity modes ϕ m,n (y), we find the corresponding modes for pressure as in Durland and Farrar (2012), using k = 0.
The first baroclinic mode is expected to be dominant, and we therefore fit to meridional modes with m = 1, although we will demonstrate that our choice of meridional modes still allows detection of the baroclinic mode 2. The baroclinic wave speeds c 1 and c 2 in the Pacific are shown in Figures 1c and 1d.These were calculated from the vertical eigenvalue Equation A1 using Galerkin methods and gridded neutral density data from the World Ocean Circulation Experiment atlas (Gouretski & Koltermann, 2010), and c 1 was verified against a similar calculation by Chelton et al. (1998).
Figure 1e shows the first six normalised meridional modes for pressure for the baroclinic mode 1 when k = 0.The black circles show the meridional locations of the TAO moorings, although in practise the meridional TAO data is much more sparse due to missing data and no moorings at some longitudes.This limits our ability to properly resolve the meridional structure, and it was concluded by FD12 that this resolution was insufficient to project the TAO data onto meridional modes -they instead focused on symmetry/antisymmetry about the equator.However, we demonstrate some success in using the TAO moorings to uncover the meridional structure of the waves.

TAO Data
We use 26 years of data (1 January 1993-31 December 2018) from the TAO mooring array, with locations spanning the equatorial Pacific (Figures 1c and 1d).Daily averages of surface dynamic height relative to 500 dbar are used throughout; the data set is described in detail in FD12.We exclude some TAO and nearby TRITON moorings on the western side of the array (not shown in Figures 1c and 1d), since there isn't coverage both sides of the equator at these longitudes.For "like-for-like" comparison with the year-long NEMO output, we take the corresponding year (1 April 2006-5 April 2007) in the TAO timeseries, labeled in plots as "TAO: sampled".When the full time series is used, it is labeled "TAO: full data".

NEMO Simulation
We use output from the NEMO v3.2 global ocean model at 1/4° horizontal resolution (Madec, 2008).The configuration and the time period of one year from 1 April 2006 to 5 April 2007 are identical to the control simulation described in B1.We use daily average fields of surface dynamic height relative to 500 dbar for consistency with the TAO data set, and the spatial domain is bounded zonally to contain the most western and eastern TAO moorings, and meridionally by ±12°N, shown as the red dashed box in Figures 1c and 1d.The model is forced using the CORE2 IAF data set (Large & Yeager, 2009), with wind stress based on the NCEP reanalysis (Kalnay et al., 1996), thus winds forcing the model are intended to represent the true wind field for the same time period.
For "like-for-like" comparison with the TAO data set, we sample the NEMO output at the locations of the TAO moorings, and remove data that are missing in the corresponding TAO time series.These data are labeled in plots as "NEMO: sampled", and the full spatial data are labeled as "NEMO: full data".

Methods
We estimate the power spectra of the dynamic height data in NEMO and TAO, first for the zonally averaged data, and then at individual longitudes.We follow FD12 and use a least squares fitting method.FD12 fitted in time and longitude; we instead fit in time and latitude.If there were no missing data, the temporal fit could be accomplished using a discrete Fourier transform, but given the nature of the mooring data a least squares fit to Fourier modes in time is more appropriate and allows us to use all of the data.This also has the advantage of allowing us to fit in latitude to the theoretical meridional basis of pressure modes explained in §2.

Least Squares Fit
We use a tapered least squares method, described in detail in FD12.The general method can be explained for 1D data; 2D (time and space) data are first flattened to one dimension.First, L modes are chosen as a basis for the data h of size N, and stored in a matrix E of size N × L. The model can then be written as where, x (size M) contains the coefficients of the modes, and n (size N) is the noise in the data.The tapered least squares solution for x minimizes the sum of squares and is given by: where, the parameter  =  2  ∕Δ 2  is the "noise to signal ratio", with   2  the expected variance in the noise, and  Δ 2  the expected variance in the coefficients of the model.

Temporal Modes
The temporal extent of the NEMO and sampled TAO datasets is 370 days.We define a full set of 185 frequencies equally spaced between ω min = 1/370 cpd and ω max = 1/2 cpd with a resolution of ω min .We then create a full set of 370 temporal modes consisting of a sine and cosine at each of these frequencies.
The full 26 year TAO time series is split into non-overlapping periods of 370 days, and the resulting spectrum for each is averaged to give the "TAO: full" spectrum.The temporal modes are thus the same as for the other datasets, allowing comparison to spectra of the shorter time series whilst providing a more reliable estimate.

Meridional Modes
The meridional modes are found from the theoretical solution for the baroclinic mode 1 (Equation A2), converted for pressure at k = 0. Whilst the meridional velocity modes are orthogonal, the pressure modes are not (Zang et al., 2002).We use the first 6 modes, as we find that this gives the best balance between representing the higher mode wave signals present in the data and minimizing overfitting associated with too many model parameters.

Fitting Procedure
Before performing the least squares fit, we apply a high pass filter at 20 days (0.05 cpd) to remove any low frequency Rossby waves (see Figure 1f) and other non-wavelike structures.The method is described in detail in Appendix B. This step is particularly important when zonally averaging the TAO data, as the average can easily be contaminated by missing data coming into/out of the record.
We then find the temporal/meridional power spectrum both for zonally averaged fields (e.g., Figure 2) and at each longitude separately (e.g., Figure 3).When taking a zonal average, we remove any times/latitudes with data at only a single longitude.The set of temporal frequencies is truncated to those higher than 0.05 cpd, since those below this have already been removed in the high pass filter, leaving 167 of the 185 original frequencies.The 167 × 2 = 334 temporal and six meridional modes are cross-multiplied and flattened to give a set of 334 × 6 = 2004 modes, and these form the matrix E. The noise to signal ratio r is taken to be 35; see Appendix C for justification.
The coefficients x of the modes are then estimated using Equation 5.There are two modes at each frequency and meridional mode number, corresponding to the temporal sine and cosine.The power is then defined as the square root of the sum of these two coefficients squared.

Results
We first directly compare the model data to the observations.Figure 1a shows the high pass filtered (see Appendix B) and zonally averaged 500 dbar dynamic height in NEMO and TAO at the equator, both sampled for a "likefor-like" comparison.Most of the time, the similarity between the fields is remarkable -the peaks and troughs often agree well, showing that the model captures much of the varying phase and amplitude of the variations.
Although the ocean model is forced by realistic winds, this level of agreement is somewhat surprising and lends credence to the idea of wave excitation being deterministic rather than stochastic.
The Pearson correlation coefficient of the zonally averaged NEMO and TAO time series at the equator (shown in Figure 1a) is 0.39, which is significant at the 1% level.This correlation decreases to ∼0.25 at a latitude of ±8° N. Agreement at longer timescales is generally better, with a correlation coefficient of 0.66 at the equator when the high pass filter cut-off is at 50 days, and 0.75 at 100 days.Figure 1b confirms that there is no obvious relationship between the amplitudes of the meridional and zonal wind stress and the dynamic height; see Durland and Farrar (2012) for a derivation of the complex dependence of wave amplitudes on wind stresses.
Note that although the dynamic height time series does not appear to agree with that of the MOC stream function at the equator from NEMO shown in B1 (their Figure 4a), this is expected since meridional velocity modes have opposite symmetry about the equator to those of dynamic height.The time series of dynamic height at the equator does not contain a contribution from the 10 day (n = 0) mode since the corresponding pressure mode is zero at the equator (see Figures 1e and 1f), whereas the 10 day mode is contained (and at times dominant) in the time series of the MOC stream function (B1, 2) since the n = 0 mode for meridional velocity has a maximum at the equator.Furthermore, the MOC time series (B1, 2) is at 1,583 m depth, whereas the 500 dbar dynamic height is calculated by integrating the specific volume anomaly over the top 500 m of the ocean, and is thus more representative of upper ocean dynamics and wave structures.
Figure 2 shows the temporal/meridional spectra of the zonally averaged dynamic height data.The black dotdashed line showing the spectrum of the full NEMO simulation exhibits well defined peaks at (or just above) the expected frequencies, particularly for the first 3 meridional modes.There are also peaks in each spectrum corresponding to other meridional modes of the same symmetry -this indicates an imperfect fit of the basis functions to the data, and is expected (for example, we are using a baroclinic mode 1, k = 0 meridional basis to represent all baroclinic modes).Due to the completeness and high spatial resolution of this data set, there is less power at the higher frequencies than in the other data sets, which predictably pick up more "noise", due to both the limited meridional resolution and the missing data in the zonal average.
The blue line shows the sampled NEMO data, and this should be compared to the red sampled TAO data for a "like-for-like" comparison between the model and observations.Both exhibit peaks at the predicted frequencies (using k = 0 in Equation 2) for at least the first 4 modes, and they show excellent similarity at most frequencies and meridional modes.In particular, they have identical peaks at both the baroclinic modes in the meridional Magenta and gray dashed lines show the predicted k = 0 frequency of the baroclinic mode 1 and mode 2 respectively, calculated using the dispersion relation (2) and the baroclinic wave speeds shown in Figures 1c and 1d, averaged over the red-dashed NEMO subdomain.modes 0, 2, and 3.There is a discrepancy of a factor of ∼2 in the amplitude of the BC1M1 and BC2M1 peaks between the TAO and NEMO data, which is consistent between both their respective 'full' and 'sampled' data sets, suggesting that this is a true difference between the model and observations.
Overall, the comparison of the zonally averaged fields in Figure 2 shows that the modeled large amplitude oscillations of the Pacific MOC found by Bell et al. (2021); Blaker et al. (2021) are largely consistent with observed amplitudes, although NEMO appears to overpredict the amplitude of the BC1M1 and BC2M1 modes by a factor of two.The comparison between the full spatial and spatially sampled datasets allows us to evaluate the usefulness of the TAO mooring array in representing these waves over the full equatorial ocean.Although the spectrum calculated using the full data (black dashed) deviates from the sampled data (blue lines) somewhat at higher where, He n is the nth Hermite polynomial, and β is the Rossby parameter.The frequency ω m,n and zonal wavenumber k satisfy the dispersion relation (Equation 2 in the main text).

Figure 1 .
Figure 1.(a) High-pass filtered and zonally averaged 500 dbar dynamic height at the equator for sampled Nucleus for European Modeling of the Ocean (NEMO) data (blue) and Tropical Atmosphere Ocean (TAO) mooring data (black).(b) High-pass filtered and zonally averaged zonal (τ x ) and meridional (τ y ) wind stress at the equator used to force the NEMO simulation.(c) Baroclinic wave speeds c 1 and (d) c 2 , calculated from World Ocean Circulation Experiment neutral density data.The red dashed box shows the subdomain of the global NEMO simulation used, and black circles show locations of TAO moorings.(e) Meridional structure of pressure modes for baroclinic mode 1, c 1 = 2.7 m s −1 , k = 0. TAO mooring latitudes shown as black circles.(f) Dispersion relation for unforced equatorial modes.The first (second) baroclinic modes are shown in black (blue dashed), for baroclinic wave speeds c 1 = 2.7 m s −1 and c 2 = 1.6 m s −1 (their values averaged over the red dashed subdomain in [c, d]).The red line shows k = 0, and the labeling convention is such that "BC2M1" corresponds to baroclinic mode 2, meridional mode 1.

Figure 2 .
Figure 2. Temporal and meridional spectrum of zonally averaged 500 dbar dynamic height.Each panel shows a meridional mode, and a rolling mean over five frequencies is shown.Red shows Tropical Atmosphere Ocean (TAO) data for the duration of the NEMO simulation, blue shows NEMO data sampled at locations and times of available TAO data.Black dot-dashed shows NEMO data with full spatial resolution, and black solid shows TAO data from the full 26 year timeseries.Magenta and gray dashed lines show the predicted k = 0 frequency of the baroclinic mode 1 and mode 2 respectively, calculated using the dispersion relation (2) and the baroclinic wave speeds shown in Figures1c and 1d, averaged over the red-dashed NEMO subdomain.

Figure 3 .
Figure 3. Zonal structure of dynamic height power spectrum, split in columns by meridional mode (top) Nucleus for European Modeling of the Ocean (NEMO) data set with full spatial resolution (upper middle) NEMO data set sampled at Tropical Atmosphere Ocean (TAO) locations and times (lower middle) TAO data set for the duration of the NEMO simulation, and (bottom) full TAO data set.The top three rows show a rolling mean over five frequencies, and the bottom row is calculated as an average over spectra of 1 year moving windows.Red and black dashed lines show the predicted baroclinic mode 1 and 2 frequencies respectively, calculated using the baroclinic wave speeds averaged in latitude between ±12°N.