Wind-Driven Oscillations in the Meridional Overturning Circulation near the equator. Part I: Numerical Models

: The great ocean conveyor presents a time-mean perspective on the interconnected network of major ocean currents. Zonally integrating the meridional velocities, either globally or across basin-scale domains, reduces the conveyor to a 2D projection widely known as the meridional overturning circulation (MOC). Recent model studies have shown the MOC to exhibitvariabilityon near-inertial timescales, and also indicate a regionof enhancedvariabilityon the equator. We present an analysis of three integrations of a global conﬁguration of a numerical ocean model, which show very large amplitude oscillations in the MOCs in the Atlantic, Indian, and Paciﬁc Oceans conﬁned to the equatorial region. Theamplitudeoftheseoscillationsisproportionaltothewidthoftheoceanbasin,typicallyabout100(200)Sv(1Sv [ m s 2 1 ) in the Atlantic (Paciﬁc). We show that these oscillations are driven by surface winds within 10 8 N/S of the equator, and their periods (typically 4–10 days) correspond to a small number of low-mode equatorially trapped planetary waves. Furthermore, the oscillations can be well reproduced by idealized wind-driven simulations linearized about a state of rest.


Introduction
The large-scale global meridional ocean overturning circulation (MOC) gained widespread recognition through its depiction as the ''Great Ocean Conveyor Belt'' (Broeker 1987). Its portrayal as such gives the impression of a slow and steady transport of water, along with heat, salt, carbon, and other tracers around the global oceans. However, observations and numerical models show that this view is highly simplified and that the MOC exhibits variability on a wide range of time scales. The past 15 years have seen an observational focus on the Atlantic MOC (AMOC) (e.g., Cunningham et al. 2007;Mercier et al. 2015;Lozier et al. 2019;Frajka-Williams et al. 2019). Sustained observations by the RAPID-MOCHA-WBTS (Western Boundary Time Series) array at 26.58N (Rayner et al. 2011) revealed much greater MOC variability than was previously known to exist, with short-term values ranging from 25 to 35 Sv (1 Sv [ 10 6 m 3 s 21 ), around a mean state of ;17 Sv (McCarthy et al. 2015(McCarthy et al. , 2012Smeed et al. 2014Smeed et al. , 2018. Ocean models can simulate a large fraction of the observed AMOC variability (e.g., Duchez et al. 2014b;Xu et al. 2014;Blaker et al. 2015;Jackson et al. 2016). Models are therefore well suited to identifying the physical processes governing AMOC variability: short-time-scale Ekman (wind-driven) variability in the surface ocean; longer-time-scale wind-driven variability through mechanisms such as Sverdrup transport due to Ekman pumping at higher/lower latitudes (Duchez et al. 2014b); changes in the eastern and western boundary densities associated with local currents (Duchez et al. 2014a) or processes such as Rossby wave excitation and propagation ; variability of the Gulf Stream and Antilles Current (Hirschi et al. 2019;Meinen et al. 2019). Models also allow us to study how strongly ''chaotic'' intrinsic processes such as ocean mesoscale eddies project onto the MOC variability and to estimate how large this contribution is compared to the predictable sources of MOC variability Grégorio et al. 2015;Leroux et al. 2018).
Model-based work also suggests the existence of MOC variability linked to near-inertial waves (Blaker et al. 2012). The most striking feature of this type of MOC variability is its amplitude, which in the Atlantic at 268N can have peak-to-peak values of 50 Sv. The existence of MOC variability induced by near-inertial waves (NIWs) is yet to be confirmed in nature. However, theoretical considerations suggest that the MOC should exhibit a resonance at near-inertial frequencies (Sévellec et al. 2013) with resonance peaks in good agreement with the MOC variability caused by equatorward propagation of wind-driven lowmode NIWs presented in Blaker et al. (2012). NIWs are excited in all ocean basins with similar amplitude, and although not explicitly shown in Blaker et al. (2012), NIW-driven MOC variability is found in the Pacific and Indian Oceans in numerical model simulations.
The MOC variability linked to near-inertial waves only becomes obvious when model output is studied at a high (i.e., subdaily) temporal frequency. In addition to the NIW-driven MOC variability first described in Blaker et al. (2012), highfrequency model output has also revealed an even larger, and Denotes content that is immediately available upon publication as open access.
to-date unexplained, MOC variability occurring around the equator in all three ocean basins. Hirschi et al. (2013) show plots of the standard deviation of the MOC computed using 5-day mean output from a numerical model simulation (see their Fig. 6). In all ocean basins the MOC shows standard deviations of a few Sverdrups over most of the latitudinal range, with a marked increase confined to the region 6108N/S. Within 28-38of the equator the standard deviation increases considerably, with values for the global ocean MOC exceeding 60 Sv. This variability is not only much larger than the variability due to Ekman transports at other latitudes, which typically dominates the high frequency MOC variability, but it also has a very different vertical structure. The maximum standard deviation of the MOC occurs at a depth of about 1500-2000 m whereas the maximum MOC variability due to Ekman transports occurs in the top 100 m (Lee and Marotzke 1998). A more recent review of the Atlantic MOC in high-resolution models shows that the enhanced equatorial MOC variability occurs in a range of different numerical models and model configurations (Hirschi et al. 2020).
The equator has proven to be particularly challenging when observational data are assimilated into numerical models to improve forecast systems. The data assimilation can result in spurious vertical velocities and overturning circulations (Bell et al. 2004;Park et al. 2017). Consequently, some reanalyses have applied specific treatment to the equatorial region to address these problems (Balmaseda et al. 2013;Waters et al. 2014).
The studies introduced so far illustrate the rich spatial and spectral structure of the MOC and demonstrate how it is helpful to understand the origins of MOC variability whether or not it contributes significantly to the time-mean heat transport.
As the NIWs that drive the Atlantic MOC variability presented in Blaker et al. (2012) propagate equatorward, their period increases, remaining slightly superinertial. The analysis in Blaker et al. (2012) did not examine the equatorial region, and the role NIWs might play in the equatorial MOC variability is unknown. The existence of equatorially trapped motions near the equator, either in the form of vertically propagating waves or vertical normal modes, is well known (e.g., Gill 1982;Vallis 2017, chapter 8). The excitation and propagation of these waves by wind forcing was popularized by Lighthill (1969), reviewed by McCreary (1985), and discussed in depth by Clarke (2008). Much of the focus of this literature is on low-frequency planetary waves and Kelvin waves but the analysis applies also to the higher-frequency inertia-gravity and Yanai waves. Farrar and Durland (2012) present an analysis of TAO/TRITON array data that reveals periods characteristic of long-wavelength inertia-gravity and Yanai waves. One of the main questions we address in this study is whether the oscillations in the MOC found in our NEMO simulations can be attributed to those waves.
In this paper we explore the characteristics and driving mechanisms of the simulated MOC variability in the equatorial region and determine whether there is a link between it and the extreme near-inertial variability known to occur at higher latitudes (Blaker et al. 2012). We perform three numerical model simulations, storing output at 4-hourly mean intervals in order to resolve the near-inertial time scale at all latitudes in order to investigate the origin of the extremely large variability in the equatorial MOC. In an accompanying paper (Bell et al. 2021, hereafter Part II), we explore in much greater depth the role equatorial planetary waves play in generating zonally integrated meridional velocity anomalies and derive a second-order ordinary differential equation that describes these anomalies. We use this to produce idealized simulations of the equatorial MOC transports seen in the numerical ocean model simulations using the surface wind forcing applied to the numerical model. The remainder of the paper is presented as follows: section 2 describes the numerical model setup and experiment, together with a brief summary of the method developed in Part II; section 3a summarizes the spatial and spectral structure of MOC variations near the equator using 4-hourly time-mean data and how they depend on local and remote wind forcing; section 3b examines the spectra of the wind forcing; section 3c analyses the spatial and spectral structure of the dynamic height and vertical velocity fields; and section 3d discusses the structure of the vertical and meridional modes in the Pacific, Atlantic, and Indian Oceans and reconstructions of the MOC calculated using the idealized simulations. We conclude with a discussion in section 4.

a. Numerical simulations
We use NEMO v3.2 (Nucleus for European Modeling of the Ocean; Madec 2008) in the global ORCA025 configuration, a tripolar grid based on the semianalytical method of Madec and Imbard (1996). Horizontal resolution is 1/48 (1442 3 1021 grid points). South of 208N the model grid is isotropic Mercator, and north of 208N the grid becomes quasi-isotropic bipolar, with poles located in Canada and Siberia to avoid the numerical instability associated with the convergence of the meridians at the geographic north pole. At the equator the resolution is approximately 27.75 km, becoming finer at higher latitudes such that at 608N/S it becomes 13.8 km. Aside from an increase in the vertical grid resolution to 75 levels, and an update to the code version, the configuration is identical to the one presented in Blaker et al. (2012). A control integration of NEMO ORCA025 starting from rest in January 1958 forms the basis of the experiment. Climatological initial temperature and salinity conditions are taken from PHC2.1 (Steele et al. 2001) at high latitudes, MEDATLAS (Jourdan et al. 1998) in the Mediterranean, and Levitus et al. (1998) elsewhere. The simulations are forced with the CORE2 IAF dataset (Large and Yeager 2009). For the period 1 April 2006-5 April 2007 the frequency of model output was increased to 4-hourly means.
Using the configuration of NEMO described above we perform two additional simulations of 6-month duration with 4-hourly mean output frequency (Table 1). In these experiments the wind forcing is held constant over a range of latitude bands. For simulation EQ_CONSTANT the wind forcing was held constant between 108S and 108N. Poleward of 158N/S the wind forcing was fully varying and identical to CONTROL. A linear ramp between 158 and 108 was applied. Similarly, for simulation EQ_VARIABLE the wind forcing was held constant poleward of 158, fully varying between 108S and 108N, with linear ramping between 158 and 108. The effect of these adjustments on the zonal wind stress variability within the tropics is shown in Fig. 1. Note that all of these simulations used relative (not absolute) wind speeds to compute the wind stress, hence the small differences between CONTROL and the experiments where the curves overlay.

b. Analysis of the meridional overturning
We compute the meridional overturning circulation C on both depth z and density s 2 coordinates. The circulation C(y, z) is computed from the zonal-and depth-integrated model velocities y using (1) The zonal and meridional coordinates are referred to by x and y, respectively. The variable H max is the maximum depth of the basin cross section, and varies as a function of latitude.
The terms x W and x E are the western and eastern boundaries of the basin, and again vary as a function of latitude. The overturning on density coordinates C(s 2 ) is computed similarly, with the vertical integral being over density instead of depth.
To better understand the geographical nature of the variability in meridional velocities that gives rise to the large transport variability we also compute the cumulative transport across the ocean basins. This is done using which is similar to Eq. (1), but retains the zonal (x) integral at each grid point to give C c (x, y, z). Interpretation of this quantity is simpler when the vertical integral is done from the surface to the ocean floor.

c. Computation of equatorial modes
Earlier work has shown that the equatorial region where the large MOC oscillations are simulated is characterized by a very low ratio between intrinsic-chaotic and total MOC variability ). This provides a strong indication that one should be able to explain the simulated equatorial MOC variability from the surface forcing. How to project wind forcing onto equatorial MOC oscillations is described in detail in Part II.
The following provides the key steps which involve solving a second-order ODE for the amplitudes of the vertical m and meridional n normal modes,ỹ m,n (t), driven by the surface winds.
Part II assumes that the fluctuations can be represented as linear motions on a stably stratified state of rest driven by surface wind stresses and the pressure forces on the eastern and western boundaries and weakly retarded by Rayleigh damping (with coefficient r). It integrates the equations of motion zonally across the ocean basin using boundary conditions of no normal flow and then projects the forcing onto the vertical and meridional normal modes. The zonal integral of the meridional velocity anomalies y 0 is expressed as where p m (z) are the nondimensional baroclinic vertical modes, u m,n (ỹ m ) are the nondimensional meridional modes, and z andỹ m are nondimensional forms of z and y: with theỹ m coordinate depending on the phase speed c m of the baroclinic mode with which it is associated. The term b is the meridional derivative of the Coriolis parameter. The amplitudesỹ m,n of the zonally integrated meridional velocity in the mth vertical and nth meridional mode evolve according to the second-order ordinary differential equation 2 (X f m,n 1 gDh f m,n ), (5) whereỸ m,n is the projection of the meridional wind stress on the mode,X f m,n and gDh f m,n are the forcings due to the zonal wind and pressure forces on the mode, respectively, and the natural, unforced and undamped, frequency v m,n ofỹ m,n is given by v 2 m,n 5 bc m (2n 1 1). Equation (5) describes a lightly damped harmonic oscillator with natural frequency v 2 m,n driven by terms derived from the wind and pressure forcing.
As c m cannot be determined perfectly, Part II calculates solutions forỹ m,n for sections of the full time series, typically of 16-day duration, that are long enough to capture significant fluctuations but not so long that the solutions ''lose track'' of the phase of the oscillation. They also explore the magnitude of the uncertainties in c m . The results shown here were derived using the ''standard'' set of parameters selected by Part II and the full wind stresses and full pressures on the boundaries to drive the solutions. The simulations are referred to as D H (of type D and kind H) in Part II, which means that simulations include a constant offset and the frequency is allowed to vary. The sensitivity of solutions to the available options is explored in some detail in Part II. The code and input files to enable the reader to reproduce the results we present and explore the choice of modes are available via Zenodo/GitHub (Bell and Blaker 2021).

a. Volume transport
We start by analyzing the 1-yr ''control'' simulation for which output was stored as 4-hourly means. In agreement with Hirschi et al. (2013), we find extremely large variability of the MOC near the equator in all ocean basins. The time mean and standard deviations of the MOC computed in depth and density space for the Pacific, Atlantic, and Indian Ocean basins are presented in Figs. 2 and 3, respectively. The maximum standard deviation on the equator is 137, 47.5, and 55 Sv in the Pacific, Atlantic, and Indian Ocean basins, respectively. Maximum variability is around 1500 m in the Pacific, and 200-300 m deeper (2-3 model levels) in the Atlantic and Indian basins. The amplitude and structure of the standard deviation are not strongly affected by the choice of vertical coordinate on which the MOC is computed.
We select a point (08N, 1583 m) and examine the MOC transport in each basin more closely. The time series of the MOC shows transport ranging between 6400 Sv in the Pacific (Fig. 4a), between 6150 Sv for the Atlantic (Fig. 4d) and from 2200 to 150 Sv for the Indian basin (Fig. 4g). The time series all show that the MOC oscillations vary in both amplitude and period. Phases with large MOC oscillations alternate with phases where the amplitude of the oscillations is greatly reduced. The period of the oscillations also changes temporally, with the largest oscillations typically being associated with longer periods. A wavelet analysis highlights the presence of dominant periods in the time series, with significant power found at 10 and 4 days in all basins (Figs. 4b,e,h). The 4-day signal is more persistent in time, while the 10-day mode is transient. Figures 4c, 4f, and 4i show the global wavelet spectrum and Fourier power spectrum for the same range of periods. Two distinct peaks appear in the global wavelet spectrum, at approximately 4-and 10-day periods. The Fourier power spectrum shows greater detail, indicating that the peaks in the wavelet are composed of more than one signal. These modes are consistent throughout the water column, as suggested by the structure of the standard deviation in Figs. 2 and 3. Performing wavelet analysis at other depths shows similar results to those presented at 1583 m.
A disadvantage of examining the basin-integrated transport C(z) is that it provides no geographical information on where the main contribution(s) to this variability come from. Computing the cumulative MOC [C c , Eq. (2)], integrating from west to east across the basin, and its zonal derivative provides insight into the zonal structure of the source of the variability (Fig. 5). The amplitude in the easternmost cell of Fig. 5a (near 808W) is total basinwide MOC, and longitudes where the gradient is steepest (indicated by Fig. 5b) show where the largest contributions to the total transport arise. Interestingly, this reveals no robust pattern for the zonal structure of the meridional transport (velocity) anomalies that manifest in the basinwide integral transport. Contributions to the transport variability are not confined to the western boundary, or to topographic features. Instead we find the large MOC fluctuations accumulate from more or less random longitudes, with some accumulating from contributions in the eastern side of the basin (east of 1208W), while others begin closer to 1808W. There are also examples (e.g., 100-120 days into the simulation, corresponding to early to mid-July) where the cumulative MOC reaches 6 100 Sv by 1408W, and then reduces again so that by 808W the amplitude is greatly reduced. The zonal derivative of the cumulative transport ( Fig. 5b) shows that contributions may be modest and consistent from several tens of degrees longitude, but there are also smaller, localized contributions which appear to be constructive interference associated with the interaction of slower westward propagation of equatorially trapped Rossby waves. On average, we find the amplitude of the full basin-integrated MOC variability to be proportional to the width of each ocean basin, with all three ocean basins yielding a scaling of approximately 0.03 Sv km 21 when using the maximum amplitude of MOC, or 0.007 Sv km 21 when using the root-mean-square (rms) amplitude computed from the MOC time series. This scaling suggests that the meridional velocity anomalies could be very small. A uniform velocity of 1 cm s 21 across the width of the equatorial Pacific (15 500 km) and to the depth of maximum overturning (;1500 m) would be sufficient to yield a transport of 230 Sv.
In the absence of clear contributions from topographically constrained currents, the scaling of the transport variability by basin width implies a role either by the local surface forcing or from high-frequency variability propagating from higher latitudes. On these short time scales the former would most likely arise from the wind fields. Previous work has shown large highfrequency MOC variability caused by NIWs excited by wind variability at higher latitudes that propagate equatorward (Blaker et al. 2012). To isolate the source of the equatorial

650
MOC variability we next examine the simulations in which the wind forcing was held constant over given latitude bands.
In the control simulation with fully varying winds the spatial and temporal behavior of the equatorial MOC variability dominates the NIW forced variability at all other latitudes (Fig. 6a). The marked increase in the amplitude of the variability aligns with the standard deviation presented in Fig. 3. The signals are not purely symmetric about the equator. Some appear to propagate across the equator, aligning with the incoming NIWs from higher latitudes, while others appear to manifest locally. The dominance of ;10-day-period signals is evident during May (days 25-60) and September (days 150-180), while at other times in the simulation (e.g., during June-August, days 60-150) the shorter-period signals dominate. These also appear to have shorter latitudinal length scales than the ;10-day-period signals.
Holding the wind forcing constant near the equator (Fig. 6b) the variability quickly weakens, fading to much smaller amplitudes after 30 days. The ;10-day-period signal evident during May (days 25-60) is still discernible, but the same is not true for the events during September (days 150-180). From early July (day 90) onward the variability in the equatorial region no longer stands out compared with higher latitudes, and the prevailing signal is that of the NIWs propagating from higher latitudes, toward and across the equator into the opposite hemisphere. The time it takes for the decay in variability to occur may be a result of a natural resonance, as identified by Sévellec and Fedorov (2013), with weak damping. It may also be that some of the signal on the equator is driven by winds 108-158 remote from the equator. In this experiment wind variability is linearly ramped down, but still is present in the region 108-158N/S. To complete the picture, the complementary simulation (EQ_VARIABLE), in which only the winds within the equatorial region are permitted to vary with time shows that an equatorial response nearly identical to the one with fully varying wind forcing can be recovered. Under locally varying wind (Fig. 6c) the signals have a wider range of frequencies of variability, more asymmetry about the equator and more evidence of meridional propagation of the large amplitude variability across the equator than in the EQ_CONSTANT simulation. Detailed comparison of Figs. 6a and 6c shows that the main events near the equator in the two experiments agree well. Together, these simulations reveal that the large equatorial variability is forced by locally varying winds (Fig. 6), while a substantially smaller fraction of equatorial variability arises due to signals propagating from higher latitudes.

b. Analysis of surface wind stress components
Our simulations in which we limit the geographical region in which the surface wind forcing is allowed to vary show that the large equatorial MOC variability is predominantly driven by local wind variability. We therefore examine the wind stress fields to try to identify and characterize the source of the large equatorial MOC variability. The zonal and meridional components and the absolute magnitude of the wind stresses on the equator are presented in Fig. 7. Overall, the fields show wind stresses up to 0.1 N m 22 , and the modulus wind stress (Fig. 7e) is dominated by the zonal wind stress (Fig. 7a). The zonal wind stress is relatively constant throughout the year, although there is a seasonal signal due to the Walker cell over the Pacific [east-west movement of the strongest wind stress located in the central Pacific that can be seen in Figs. 7a and 7e, with the most westward extent in April (days 0-30) and most eastward extent in October (days 180-210)]. The meridional wind stress (Fig. 7c) peaks during the height of boreal summer (days 90-180), and reaches a minimum during boreal winter (days 240-330). Also noticeable is a ''banding'' that is near-constant in time and stationary in longitude, which is most clearly seen as vertical stripes in Fig. 7c. There are approximately three ''bands'' per 508 of longitude. These are symptoms of the T21 spectral atmosphere that was used to produce the CORE2 reanalysis. At higher frequency it is possible to make out predominantly westward propagating signals which cross the Pacific in around 1 month. At a fixed longitude these waves would appear to create signals with a time scale of a few days. To show these high-frequency signals more clearly, Figs. 7b and 7d display fast Fourier transforms (FFT) of the wind stress components. The power spectrum is computed for each longitude and presented showing the power (shaded) for periods between 3 and 18 days. This does not reveal any spectral peaks for periods between 4 and 10 days. The same computation was performed for wind stress curl (not shown) and this again exhibits no enhanced power for 4-10-day periods. The correlation coefficient between the time series of the Pacific zonal mean meridional wind stress and the MOC in the Pacific is 20.0244. Experimenting with smaller zonal means (east, central, and west Pacific) and with lags up to 67 days we found the correlations never exceed 0.2, so we conclude that there is essentially no correlation between any of the wind stress components on the equator and the MOC.

c. Analysis of ocean vertical velocity and SSH fields
Finding no clear link between the surface wind stress and the large equatorial MOC variability we now look at vertical velocities W and sea surface height (SSH). As stated earlier, there is a high degree of coherence in the variability throughout the water column, with the maximum variability centered around 1500-2000 m, reducing toward the bottom and surface of the ocean. Analysis in Blaker et al. (2012) found unfiltered vertical velocity W revealed some of the cleanest signals of NIWs. Here we compute power spectra (FFTs) for vertical velocity W at 500-m, 1-km and 2-km depth (Fig. 8). Distinct peaks in the power spectra are visible at all depths, with W 2km exhibiting the most energy, which is consistent with the vertical structure of the equatorial MOC variability (Fig. 2). There are peaks at 7-(;0.14 cpd), 5.5-(;0.18 cpd), and 4-day (;0.25 cpd) periods, as well as some higher-frequency modes of variability (around 3.3 and 3 days), which are most distinct in W 500m and W 1km . The 5.5-and 4-day modes dominate the W 500m , while the lower and higher frequencies appear deeper in the water column. The 7-day period is relatively weak in W 500m , but in W 2km it is the dominant period. To identify the surface signature of the equatorial MOC variability we compute power spectra for sea surface height anomalies (SSHA) and 500-m dynamic height (Fig. 9). SSH anomalies were calculated by removing a time-mean from the SSH for each longitude. Both the SSH and 500-m dynamic height exhibit distinct peaks that equate to 7-(;0.14 cpd), 5.5-(;0.18 cpd), and 4-day (;0.25 cpd) periods. As with W 2km , the 7-day period is also dominant in the 500-m dynamic height power spectrum. Signals in the Pacific are largely confined to the eastern half of the basin, and moving eastward from the center of the basin toward South America the peak period lengthens. The peak power in the Atlantic is located in the center of the basin, perhaps suggesting that they are associated with topography (the Mid-Atlantic Ridge). In the Atlantic and Indian basins there is no evidence of a change in the period of the peak power as a function of longitude.
While SSHA and the 500-m dynamic height exhibit more power at periods longer than 15 days, the power spectra of W shows no significant amplitude beyond a 10-day period. This reflects the role of different dynamics that influence the dynamic height fields, processes such as a long-term trend, seasonal cycle, and an imprint of processes such as atmospheric Yanai waves. Vertical velocities are instead dominated by the faster wave or wavelike processes in the ocean.
To obtain a better picture of the spatial structure of the dominant periods of variability in the region the dominant period is identified for each grid box for W 2km , SSHA, and 500-m dynamic height (Fig. 10). A zonal mean of the equatorial Pacific region for each variable is presented in Fig. 10d. The three variables yield different spatial patterns. The vertical velocity W 2km exhibits a period of around 2 days at 108N/S, and this increases toward the equator, peaking at around 5.25 days.
The dominant period for SSHA has an interesting structure. Outside of 6108N/S this generally defaults to (or is close to) the longest period in the range selected. This characteristic is most likely an artifact of the processing or the choice of window, and is not of concern for the equatorial region. Close to the equator some interesting meridional structure is revealed, which again becomes more apparent when presented as a zonal mean (Fig. 10). A well-defined peak in the period at around 5.4 days exists within 28 of the equator. Then, symmetric about the equator we find a pair of peaks, one in each hemisphere between 38 and 68N/S. There is longitudinal variability to this pattern, and similar structure in the Atlantic and Indian Ocean basins. The structure identified in the Pacific bears some similarity to the meridional structure of some of the equatorially trapped NIWs which Farrar and Durland (2012) identify as dominant in their analysis of the TAO/TRITON array. The dominant periods identified in the 500-m dynamic height (Fig. 10c) appear to show behavior similar to W 2km outside of the equatorial region (i.e., following the inertial period for the most part), and behavior similar to the SSHA within the equatorial region. The three peak structure is even more clearly evident in the 500-m dynamic height for the Atlantic and Indian Ocean basins than it is for the Pacific. Taking Fig. 10 as a whole, it would appear that on the equator there is a mode of variability that is strongly present at most/all depths, with a period close to that expected by theory. Away from the equator there appears to be discrete banding, and this fits with the meridional structure of equatorially trapped NIWs identified by Farrar and Durland (2012), although the absence of this banding from the W 2km analysis would suggest that it is confined more to the upper layers of the ocean. The nature of this analysis (selecting only the largest peak) masks the structure, or possible existence thereof, of secondary or tertiary peaks in the same region.

d. Excitation of overturning transports by wind-driven equatorially trapped waves
The existence of the structure seen in the analysis showing dominant periods leads us to investigate in more detail the family of equatorially trapped waves investigated by Farrar and Durland (2012). Theoretical dispersion curves in frequency-wavenumber space are computed for equatorially trapped waves (Fig. 11).
Focusing on values for long zonal wavelengths, with near-zero zonal wavenumber, the first baroclinic mode Yanai wave (solid blue curve) has a period of 9 days, while the first and second baroclinic mode, meridional mode-1 Poincaré waves (solid and dashed black curves) have periods of 5.25 and 7 days, respectively.
Using April-mean basin average density profiles at the equator for each ocean basin we compute the structure of the vertical modes. The first three vertical modes are presented in Fig. 12a. Mode 1 (the first baroclinic mode) has a single zero crossing at around 1200-m depth in the Pacific Ocean. The crossing is slightly deeper in the Atlantic and Indian Ocean basins. Mode 2 has two zero crossings, one around 200-m depth and a second around 1800 m. The deeper of these two crossings is ;100 m deeper in the Atlantic and Indian Ocean basins. Whereas mode 1 is largest at the surface in the Pacific Ocean, it is mode 2 that has the largest surface amplitude in the Atlantic and Indian Ocean basins. Mode 3 has three zero crossings at 150, 700, and 2200 m in the Pacific Ocean, again with the deepest of these ;100 m deeper in the Atlantic and Indian Ocean basins. It is possible to calculate higher modes, with each subsequent mode adding an additional zero crossing, but for clarity these are not plotted. It is worth noting that in Part II we tested the impact of using either ''classic'' baroclinic modes or ''surface modes'' (LaCasce 2017). We found that using surface modes resulted in a poorer fit to the NEMO simulations. The option to use either kind of vertical mode structure is included in the code repository. Individual wave modes comprise of both a vertical and a meridional structure. The first six meridional mode structures are also computed for each basin and each vertical mode. Figures 12b-d show the meridional mode structures in the Pacific for vertical modes 1, 3, and 5. The modes shown in Fig. 12b use vertical mode 1 and correspond to the solid curves shown in Fig. 11. They can be separated into two classifications: those which are symmetric about the equator (even number modes) and those which are antisymmetric about the equator (odd number modes). Mode 0, corresponding to the Yanai wave, is a single Gaussian shaped bump centered on the equator that decays to zero over 108. Away from the equator higher meridional modes provide larger contributions than low number modes. As the vertical mode number increases, the structure of the meridional modes contracts toward the equator, as illustrated by Figs. 12b-d.
To understand the role of these waves in generating the large-amplitude oscillations in overturning we simulate the zonally integrated meridional velocity field that arises from excitation of these wind-driven waves. This is done by solving the second-order ODE forỹ m,n (t) [Eq. (5)], driven by the surface winds for a range of vertical and meridional modes. Wind stresses and boundary pressures from our control model simulation are given as inputs to the calculation. The complete time series is split into sections of 400-h duration, and the amplitude ofỹ m,n (t) and its first time derivative at the start of each period is chosen to give the best fit toỹ m,n (t) in that period. The method is summarized in section 2, although the reader is referred to Part II in which a thorough description and analysis of the method is presented. Sensitivity of the solutions to changes in the density profile due to season and longitudinal location within the basin is explored in Part II. For the leading modes this sensitivity is of order 5%, and less for higher modes. The choice of the number of vertical and meridional modes presented here is arbitrary. Some modes are more readily excited and contribute more strongly to the solution. However, the importance of a given mode's contribution to the solution is dependent on both latitude and depth, as indicated in Fig. 12. Figure 13 shows idealized simulations of the MOC at 1583 m in the Pacific as a function of latitude and time. Figures 13a-d are idealized simulations using three vertical modes and one, two, three, and six meridional modes, respectively. For comparison, the MOC simulated by NEMO is presented in Fig. 13e. Using only a single meridional mode (Fig. 13a) it is possible to see the role of the first baroclinic mode Yanai wave, which sets the 10-day period. The larger-amplitude peaks and troughs in (Fig. 13a) correlate well with the regions of significant power in the 8-12-day scale-averaged power and the wavelet analysis presented earlier (Fig. 4). Note also that the signals shown in Fig. 13a are symmetric about the equator. Adding the second meridional mode introduces marked asymmetry to the solution, and there are now points in time (e.g., 240-250 days) which show much stronger MOC signals than with a single meridional mode. Progressively adding more meridional modes (Figs. 13c,d), not only does the detail close to the equator increase, but the strength of MOC anomalies away from the equator increases also. Comparison with Fig. 13e shows that it is possible to recover the variability of the MOC simulated by the numerical model to an extremely high degree with a small number of vertical and meridional modes.
It is possible to see how the relative contribution of a given mode to the solution will depend on the location of interest. For example, the higher meridional, lower vertical modes contribute more to the solution away from the equator (e.g., poleward of 58N/S), while the higher vertical, lower meridional modes contribute more at greater depth. To illustrate this, Fig. 14 shows the correlation between the true MOC in the equatorial Pacific and the wind-driven solution for three vertical and three meridional modes (top) and for six vertical and six meridional modes (bottom). Using three vertical and three meridional modes is sufficient to yield a convincing idealized simulation at the equator, with correlations at the equator around 0.9 in all basins (see also Table 2). Increasing the number of vertical and meridional modes to six, the largest increases in correlation occur between 58 and 108N/S and at depths below 2 km. Improvement in the idealized simulation of the MOC due to the additional vertical and meridional modes is marginal directly on the equator. We note there is an asymmetry about the equator to the pattern correlation, which is perhaps due to the relative role of different modes and the shape of the basin. Further, we find there to be some sensitivity to the choice of maximum depth used to fit the modes for the simulations. Increasing the depth from H 5 4191 m to H 5 5596 m improves the representation of the lower vertical modes and the corresponding correlations below 3500 m, but the fit for some of the higher modes is slightly worse, reducing the correlations in the upper and middepths. Further detail on how the choice of maximum depth affects the vertical modes is presented in section 4b of Part II.
Finally, we sum the contributions to the meridional velocity given by the first three vertical and meridional modes and construct the overturning transport in the simulation for each ocean basin and plot on the same axes as the MOC time series from NEMO in each ocean basin (Fig. 15). Increasing the number of modes to six vertical and meridional, or plotting at different depths, yields a very similar set of curves. To quantify the quality of the idealized MOC simulations, the rms amplitude of the equatorial MOC at depths of 1 and 2 km, together with the rms amplitude for idealized simulations using three vertical and meridional modes and six vertical and meridional modes are presented in Table 2. The fraction of variance explained (r 2 ) between the idealized simulations and true MOC time series are also given in brackets. Including three vertical and three meridional modes the idealized simulation yields 85% of the true MOC rms amplitude at 1-km depth, reducing to 80% at 2 km, and this is consistent across the three basins. The fraction of variance explained is around 0.85 at both depths and across all three basins. Increasing the number of vertical and meridional modes to six raises the rms of the idealized simulation MOCs to around 92% (1 km) and 87% (2 km), and the variance explained increases to around 0.9. Again these results are consistent across all three basins.

Conclusions
Enhanced variability of the MOC at the equator has been noted in previous model studies (e.g., Hirschi et al. 2013Hirschi et al. , 2020. To investigate the mechanisms behind this localized increase in variability we perform an experiment using NEMO, a widely adopted numerical ocean model used in ocean-only applications, including data assimilation, and a variety of coupled configurations such as the U.K. Met Office and ECMWF forecasting systems, as well as numerous CMIP5 and CMIP6 contributions. Examining the CONTROL simulation, we find that the large equatorial MOC variability observed in all basins exhibits peaks in the power spectra that correspond to periods of 10 and 4 days. Applying a time-invariant wind forcing over the equatorial region in numerical simulations results in a decay of this variability over a period of 3 months. Performing the complementary experiment, in which wind forcing away from the equator is held constant, results in equatorial MOC variability that is nearidentical to the fully forced simulation. These experiments demonstrate that the large equatorial MOC signals are forced by local (or near-local) wind variability. However, FFT analysis of the surface wind fields yielded no pronounced peaks in the power spectra to match those seen in the MOC. This suggests that the mechanism by which the wind variability drives the strong MOC oscillations involves resonant responses. Analysis of vertical velocities W, SSHA and 500-m dynamic height on the equator shows that these fields display prominent peaks at periods of 7, 5.5 and 4 days. These periods correspond to excitation of the equatorially trapped near-inertial (Poincaré) waves with vertical mode 2, meridional mode 1 (7 days), vertical mode 1, meridional mode 1 and vertical mode 2, meridional mode 2 (both around 5.25 days), and vertical mode 1, meridional mode 2 (4 days), and agree well with both theory and the analysis of observational 500-m dynamic height data from the TAO/TRITON array by Farrar and Durland (2012). The analysis of these fields also indicates that the contribution to the variability comes mainly from the eastern half of the Pacific basin, the central (Mid-Atlantic Ridge) region of the Atlantic, and contributions in the Indian Ocean are a combination of central and eastern sectors of the basin. The correspondence between the dominant frequencies seen in the vertical velocities, SSH, and dynamic height simulated in our model and the observational findings by Farrar and Durland (2012) is encouraging. It indicates that the model is simulating a real physical phenomenon. What is not evident in these fields is the enhanced power at the 10-day period. This is because this is the period of the first baroclinic Yanai mode. This mode has a meridional velocity profile which is symmetric about the equator. Its vertical velocity and dynamic height profiles are antisymmetric about the equator and hence zero there.
Analysis of the numerical model experiments, simple scaling arguments and the meridional velocities obtained through projection of low-mode waves all suggest that the strong overturning circulation identified in NEMO requires only small meridional velocities (order 1-2 cm s 21 ) distributed over a large fraction of the basin width, values which would be difficult to measure. To better understand the linkage between the time-varying winds and the large MOC oscillations in the model we simulate the meridional velocity field driven by low baroclinic and vertical mode equatorial planetary waves, and from this construct the wave induced overturning. We show that by exciting a small number of these modes with the surface wind stresses supplied to the numerical model we are able to simulate MOC transports of the correct amplitude and explain 90%-95% of the variance. Only a small number of vertical and meridional mode waves are required to obtain extremely good idealized simulations of the MOC variability on the equator. Including higher vertical and meridional modes significantly improves the idealized simulations away from the equator and at greater depths.
Part II defines slowly varying solutions toỹ m,n , using the period of the resonant frequency to create a running mean of the time series. These produce substantial volume transports on monthly to seasonal time scales, and consequently may result in non-negligible transports of heat. These slowly varying solutions therefore warrant further study.
Recent studies also suggest that the effects of the largeamplitude, high-frequency variability may be more than simply displacement of isopycnals. For instance, contributions to abyssal diapycnal mixing by equatorially trapped waves (Delorme and Thomas 2019) and indications of enhanced diathermal heat fluxes at the base of the mixed layer on the equator (Holmes et al. 2019) both point to changes in the background stratification of the ocean in this region. Results such as these would therefore imply a connection between the near-linear, highfrequency MOC variability presented in this study and the larger-scale, low-frequency circulation that arises from the background stratification.
Since the surface wind stresses imposed on the numerical model originate from observed wind fields, and we are able to achieve extremely good agreement between the MOC transport derived from equatorial planetary wave theory and that from our forced numerical ocean model simulation, we believe the modeled high-frequency MOC transports are likely to be representative of the true wind-driven overturning in the equatorial Pacific, Indian, and Atlantic Oceans on time scales shorter than 20 days.  2. The rms amplitude of the model MOC transport at the equator for depths of 1 and 2 km and the equivalent rms amplitude from idealized MOC time series using three vertical and three meridional modes (3v3m) and six vertical and six meridional (6v6m) modes (units are Sv). The fraction of variance explained (r 2 ) by the idealized simulations is given in parentheses.