Audrey Delpech
email: audrey.delpech@legos.obs-mip.fr

Claire Ménesguen

Yves Morel

Leif Thomas

Frédéric Marin

Sophie Cravatte

Sylvie Le Gentil

Intra-annual Rossby waves destabilization as a potential driver of Low Latitudes Zonal Jets -Barotropic dynamics

At low latitudes in the ocean, the deep currents are shaped into narrow jets flowing eastward and westward, reversing periodically with latitude between 15 • S and 15 • N. These jets are present from the thermocline to the bottom. The energy sources and the physical mechanisms responsible for their formation are still debated and poorly understood. This study explores the role of the destabilization of intra-annual equatorial waves in the jets formation process, as these waves are known to be an important energy source at low latitudes. The study focuses particularly on the role of barotropic Rossby waves as a first step towards understanding the relevant physical mechanisms.

It is shown from a set of idealized numerical simulations and analytical solutions that Non-Linear Triad Interactions (NLTI) play a crucial role in the transfer of energy towards jet-like structures (long waves with short meridional wavelengths) that induce a zonal residual mean circulation. The sensitivity of the instability emergence and the scale selection of the jet-like secondary wave to the forced primary wave is analyzed. For realistic amplitudes around 5-20 cm s -1 , the primary waves that produce the most realistic jet-like structures are zonally-propagating intra-annual waves with periods between 60 and 130 days and wavelengths between 200 and 300 km. The NLTI mechanism is a first step towards the generation of a permanent jet-structured circulation, and is discussed in the context of turbulent cascade theories.

Introduction

The deep equatorial and tropical circulation is organized into systems of alternating eastward and westward jets [START_REF] Firing | Deep zonal currents in the central equatorial pacific[END_REF][START_REF] Firing | Equatorial subthermocline currents across the pacific[END_REF][START_REF] Johnson | Temporal and spatial structure of the equatorial deep jets in the pacific ocean[END_REF][START_REF] Ollitrault | Zonal intermediate currents in the equatorial atlantic ocean[END_REF][START_REF] Ascani | Deep equatorial ocean circulation induced by a forced-dissipated yanai beam[END_REF][START_REF] Cravatte | Intermediate zonal jets in the tropical pacific ocean observed by argo floats[END_REF]Ollitrault and Colin de Verdière 2014;[START_REF] Qiu | Generation of the north equatorial undercurrent jets by triad baroclinic rossby wave interactions[END_REF][START_REF] Cravatte | Subthermocline and intermediate zonal currents in the tropical pacific ocean: Paths and vertical structure[END_REF]. We distinguish in particular: (i) meridionally alternating off-equatorial jets with a meridional scale of ∼ 3 • within the 15 • S-15 • N latitude range, which have a large vertical scale (quasi-barotropic), called Extra Equatorial Jets (EEJs); and (ii) vertically alternating equatorial jets trapped equatorward of ∼ 2 • , which have a small vertical scale of ∼ 350m, called Equatorial Deep Jets (EDJs). These jets are present from below the thermocline down to at least 2000 m, with some evidence of jets extending to the ocean bottom [START_REF] Firing | Equatorial subthermocline currents across the pacific[END_REF]Delpech et al. 2020a).

The underlying physical mechanisms generating these systems of jets are still poorly understood.

Different physical mechanisms have been proposed to explain their formation, relying on a cascade of mechanisms transferring energy from a deep energy source (generally generated through the propagation at depth of atmospheric variability or currents instabilities) to the mean jet-structured circulation (see Fig. 2 of [START_REF] Ménesguen | Observations and mechanisms for the formation of deep equatorial and tropical circulation[END_REF].

Earlier studies have shown that two-dimensional turbulence induces an inverse cascade, with energy transferred towards larger scales. On a rotating planet, this cascade results in the emergence of zonal structures due to the anisotropy of the Coriolis parameter (i.e. the beta effect) [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF][START_REF] Vallis | Generation of mean flows and jets on a beta plane and over topography[END_REF]. Most of these studies have however been conducted in a theoretical framework and have not addressed the source of the turbulence. The initiation of an inverse cascade at depth at low latitude in the oceans would indeed require the presence of strong turbulence which has not been documented yet.

One of the major energy sources present in the ocean at low latitude is associated with planetary waves. It has been shown that in the Atlantic and Pacific ocean, a large source of energy in the deep is associated with annual and semi-annual Rossby waves, as well as intra-annual waves: 30-day, 1000-km mixed Rossby-gravity waves, associated with surface Tropical Instability Waves; and short-scale variability, with periods around 70 days and wavelengths around 500 km [START_REF] Bunge | Variability at intermediate depths at the equator in the atlantic ocean in 2000-06: Annual cycle, equatorial deep jets, and intraseasonal meridional velocity fluctuations[END_REF][START_REF] Von Schuckmann | Generation of tropical instability waves in the atlantic ocean[END_REF][START_REF] Tuchen | Deep intraseasonal variability in the central equatorial atlantic[END_REF][START_REF] Lyman | Tropical instability waves as a resonance between equatorial rossby waves[END_REF][START_REF] Lyman | Distinct 17-and 33-day tropical instability waves in subsurface observations[END_REF][START_REF] Kessler | The annual wind-driven rossby wave in the subthermocline equatorial pacific[END_REF][START_REF] Eriksen | An estimate of equatorial wave energy flux at 9-to 90-day periods in the central pacific[END_REF][START_REF] Farrar | Intraseasonal variability near 10 n in the eastern tropical pacific ocean[END_REF][START_REF] Farrar | Barotropic rossby waves radiating from tropical instability waves in the pacific ocean[END_REF][START_REF] Farrar | Wavenumber-frequency spectra of inertia-gravity and mixed rossby-gravity waves in the equatorial pacific ocean[END_REF][START_REF] Lee | Monitoring and interpreting the Tropical Oceans by Satellite Altimetry[END_REF]Delpech et al. 2020b) (See Fig. 1 for an overview of the spectral characteristics of the observed variability). Equatorial waves, when generated at the surface, propagate to depth. In the intra-annual range in particular, waves have steep propagation ray paths [START_REF] Cox | Generation and propagation of 30-day waves in a numerical model of the pacific[END_REF] and are thus very efficient at transferring energy to the deep ocean.

These waves can develop instabilities and the destabilization of some particular planetary waves have been shown to be a potential mechanism for the formation of a jet-structured circulation [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF][START_REF] Hua | Destabilization of mixed rossby gravity waves and the formation of equatorial zonal jets[END_REF][START_REF] Connaughton | Modulational instability of rossby and drift waves and generation of zonal jets[END_REF][START_REF] Qiu | Generation of the north equatorial undercurrent jets by triad baroclinic rossby wave interactions[END_REF]d'Orgeville et al. 2007;[START_REF] Ménesguen | Dynamics of the combined extra-equatorial and equatorial deep jets in the atlantic[END_REF][START_REF] Ascani | Deep equatorial ocean circulation induced by a forced-dissipated yanai beam[END_REF]. In an idealized model configuration forced with different intraannual baroclinic Yanai waves (Fig. 1 grey dots), d [START_REF] Orgeville | Equatorial deep jets triggered by a large vertical scale variability within the western boundary layer[END_REF] and [START_REF] Hua | Destabilization of mixed rossby gravity waves and the formation of equatorial zonal jets[END_REF] were able to reproduce the EDJs-like structures, whose vertical scale was mainly dependent on the initial wave period. In a similar configuration forced with baroclinic Yanai waves and short barotropic Rossby waves at 50-day period (Fig. 1 magenta dots), [START_REF] Ménesguen | Dynamics of the combined extra-equatorial and equatorial deep jets in the atlantic[END_REF] were able to reproduce the EDJs and the first meridionally-alternating EEJs at 2 • S and 2 • N in equatorial dedicated experiments. [START_REF] Ascani | Deep equatorial ocean circulation induced by a forced-dissipated yanai beam[END_REF] were able to reproduce meridionally-alternating EEJs in the 5 • S -5 • N latitude band, with a vertically propagating beam of a Yanai wave with a period of 30 days and a wavelength of 1000 km (Fig. 1 red dot). Finally, in an idealized reduced gravity model forced with a baroclinic mode 1 annual Rossby wave (Fig. 1 blue dot), [START_REF] Qiu | Generation of the north equatorial undercurrent jets by triad baroclinic rossby wave interactions[END_REF] were able to reproduce meridionally-alternating EEJs in the 5 • N-40 • N latitude band. The configurations and characteristics of the wave forcing as well as the physical interpretation of the processes involved in these numerical experiments differ, however, from study to study. A general framework for interpreting the effects of planetary waves instability on the mean circulation is still missing. In the framework of quasi-geostrophic dynamics, [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF] studied the instability of barotropic Rossby waves on a beta plane, focusing on two particular asymptotic cases: the case of very weak nonlinearities (corresponding to resonant triad interactions) and the case of very strong non-linearities. He showed that in the limit of strong non-linearities, Rossby waves can transfer their energy towards a zero zonal wavenumber, a non-zero meridional wavenumber and zero frequency mode, the socalled zonal jet-like structures. [START_REF] Hua | Destabilization of mixed rossby gravity waves and the formation of equatorial zonal jets[END_REF] adapted the barotropic theory of [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF] to baroclinic equatorially-trapped waves. They showed that equatorially-trapped baroclinic waves can also be destabilized into vertically-alternating and meridionally-alternating jet-like structures, thus extending Gill's results to lower latitudes. [START_REF] Ménesguen | Dynamics of the combined extra-equatorial and equatorial deep jets in the atlantic[END_REF] interpreted the destabilization of the barotropic Rossby waves into EEJs-like structures in their simulation using the theory of [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF] in the strong non-linear limit, while [START_REF] Qiu | Generation of the north equatorial undercurrent jets by triad baroclinic rossby wave interactions[END_REF] interpreted the formation of EEJs in their simulations as the destabilization of annual Rossby waves through resonant triad interactions (corresponding to the weak non-linear limit of Gill's theory) and to further non-linear adjustments involving potential vorticity fluxes. In addition, all the numerical studies exploring the formation of EEJs from waves have focused on a limited set of wave periods. A general continuous framework to analyze the potential contribution of all planetary waves in generating the mean jet-structured circulation as a function of their characteristics has yet to be formulated.

In this study, our goal is to provide new insights into the mechanisms that drive the deep tropical circulation and its jet-structuring, focusing on the formation of EEJs-like structures by the destabilization of planetary waves. In particular, we aim to extend previous numerical experiments to other type of waves, with a focus on the intra-annual waves, and to create a general theoretical framework to compare and interpret the different experiments. This study focuses on the early stage of jet formation. We thus emphasize the mechanism of energy transfer from a primary energy source to EEJs-like structures with quasi-barotropic zonal velocities that reverse over a meridional scale of ∼350 km and that have a long temporal coherence. We leave for further studies the long-term equilibration of these jet-like structures. We will in particular address the following questions: (1) What is the potential for intra-annual waves to create off-equatorial, meridionally-alternating, jet-like structures ? (2) Is there a preferential frequency and wavenumber of the waves for reproducing realistic jets ? (3) What are the processes responsible for the transfer of energy from waves to jet-like structures ?

To answer these questions, we perform idealized numerical experiments of the equatorial ocean forced by an oscillating wind-stress that acts as a wave maker and generates waves with predetermined characteristics. The idealized configuration is designed to evaluate the sensitivity of the response of each forced wave to the period and wavenumber of the forcing. As a first step, we will study the destabilization of barotropic waves using two-dimensional numerical simulations, leaving the generalization to fully three-dimensional dynamics for subsequent work. The barotropic mode offers the advantage that it is not equatorially-trapped, which makes the problem invariant with latitude. It also radiate energy rapidly away from forcing regions, making it easier to simulate freely propagating waves. Evidence for purely barotropic waves in the equatorial oceans are limited to a few studies [START_REF] Farrar | Barotropic rossby waves radiating from tropical instability waves in the pacific ocean[END_REF][START_REF] Rohith | Basin-wide sea level coherency in the tropical indian ocean driven by madden-julian oscillation[END_REF][START_REF] Farrar | Long-distance radiation of rossby waves from the equatorial current system[END_REF]. The approach described in this study can however be generalized to baroclinic modes.

We show that the formation of jet-like structures is triggered by the non-linear terms in the momentum equations, which allow for wave-wave interactions. A full description of the non-linear physics is not possible, due to its intrinsic chaotic nature. We thus investigate how the waves can destabilize and produce jet-like structures within the framework of a truncated non-linear wave interaction theory. This theory, which can be viewed as a generalization of [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF], has been fully described and investigated with applications in geophysics and plasma physics [START_REF] Connaughton | Modulational instability of rossby and drift waves and generation of zonal jets[END_REF] and is here adapted to study the instability of planetary waves. It allows for solving analytically the non-linearity of the system, assuming a limited number of wave interactions, with a continuity from a weakly non-linear regime to a strongly non-linear regime, for any primary wave spectral characteristics.

The remainder of this paper is organized as follows. The numerical model configuration, forcing and wavefield used in the experiments is described in Section 2. The spectral method using wavelets transform that is used to analyze the waves in the simulations is described in Section 3. The results of the numerical experiments are described in Section 4. The instability that develops in the non-linear simulations is compared to the prediction from the 3-mode truncation Non-Linear Triad Interaction (NLTI) analytical model in Section 5. The sensitivity of the creation of jet-like structures to the wavenumber and frequency of the forced wave is investigated in Section 6. We end the article with a discussion of how the results support and complement previous studies, and summarize the main results and provide perspectives for future work.

Numerical Experiments

a. Model description

The idealized simulations performed in this study are run with the Coastal and Regional Community model (CROCO). The CROCO model solves the primitive equations (Shchepetkin andMcWilliams 2005, 2009) with a time splitting method between the fast barotropic mode and the slow baroclinic modes. The model is used in its barotropic configuration, which solves the vertically-integrated momentum equations (Shchepetkin and O'Brien 1996;[START_REF] Shchepetkin | Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations[END_REF].

(ℎ ) + (ℎ 2 ) + (ℎ ) -ℎ + ℎ = F + D (1a) (ℎ ) + (ℎ ) + (ℎ 2 ) + ℎ + ℎ = F + D (1b) ℎ + (ℎ ) + (ℎ ) = 0 (1c)
where ( , ) is the barotropic velocity vector, ℎ( , , ) = + ( , , ) is the total water depth, is the free surface elevation, is the water depth at rest, is the gravity constant and the Coriolis parameter. F and D represent forcing and frictional terms, respectively.

b. Model configuration

The basin configuration is meant to represent an idealized ocean at low latitudes. Its size is 140 • ×70 • ×5000 m (longitude × latitude × depth) and it is centered on the equator. The horizontal resolution is 0.25 • ×0.25 • . The Coriolis parameter follows the variations of an equatorial betaplane, = where = 2.3 × 10 -11 (ms) -1 is the planetary vorticity gradient and the meridional position from the equator.

The model is forced with a surface stress / 0 = ( , ), which is incorporated in the model equations (Eq. 1) as a surface momentum flux (SMF) F = / 0 , F = / 0 . As our objective is to generate freely propagating intra-annual waves, the SMF is localized inside a specific region and takes the form of a wavemaker (Eq. 2).

= -0 ( ) ( ) (k 0 • x -0 )( 0 / 0 ) (2a) = -0 ( ) ( ) (k 0 • x -0 ) (2b)
where ( ) and ( ) are envelope profiles. With this choice of wavemaker, will generate a plane wave of streamfunction = 0 (k 0 • x -0 ). Unless specified differently, we take 0 = 0 for all experiments, which leads to = 0. The generated waves thus have a strong meridional component of the velocity and propagate zonally, as are the characteristics of the observed intra-annual waves.

The SMF wavenumber 0 , frequency 0 and amplitude 0 differ for each experiment (Table 1).

The containment of the wavemaker inside a localized forcing region ( ) ( ) is essential to have an unforced region where the different waves can freely propagate without interference by the forcing. It is also more realistic as most intra-annual waves observed in the ocean are found to be locally generated and vary spatially (e.g., [START_REF] Tuchen | Deep intraseasonal variability in the central equatorial atlantic[END_REF]Delpech et al. 2020b). The envelope of this forcing region is defined following Eq. 3, it is similar to a Blackman window of size ( , )

centered on ( 0 , 0 ).

( ) = tanh -0 -0.5 -tanh -0 + 0.5 (3a) ( ) = tanh -0 -0.5 -tanh -0 + 0.5 (3b)
where ( 0 , 0 ) is the position of the center of the forcing region, and are the widths of the forcing region, and and are the tapering extent of the forcing region. Unless specified differently, the forcing region is centered in the middle of the basin ( 0 , 0 ) = (70 • , 0 • ). The other parameters are taken by default as = 20 • , = 15 • , = 5 • , = 2 • . Note that the exact location of the forcing region does not influence the outcome of the simulations. As the barotropic Rossby deformation radius is very large, the propagation of barotropic Rossby waves does not depend on latitude and are not subject to beta dispersion [START_REF] Schopf | Beta-dispersion of low-frequency rossby waves[END_REF]). Test simulations have been conducted with 0 = 15 • , without substantial modifications to the results (not shown). Similarly, the conclusions of this study are not sensitive to the extension of the forcing region. Test experiments have been conducted with = 5 • and = 25 • (not shown). The size of the forcing region controls the energy input into the ocean and changes the amplitude of the wave response. This sensitivity is thus similar to the sensitivity to the SMF amplitude ( 0 ), which is fully investigated in Section 6. The characterization of the oceanic response to the forcing is investigated in Section 2c and examples are given in Section 4 and illustrate that for the range of parameters considered, the envelope does not affect the wave response.

All lateral boundaries are closed with free slip boundary conditions. A sponge layer is implemented on the northern and southern boundaries to avoid the propagation of artificial coastal waves. In the sponge layers, the lateral viscosity is increased towards the northern and southern boundaries. The bottom is flat and a linear bottom drag law is applied to balance the energy input from the forcing. Therefore, the frictional terms in the equations are: D = -+ 0 (tanh(( -)/ ) + tanh(( + )/ ))/2, D = -+ 0 (tanh(( -)/ ) + tanh(( + )/ ))/2 , with = 1.5 × 10 -4 ms -1 and 0 = 5000 m 2 s -1 , = 34 • and = 1 • . No explicit horizontal viscosity is included in the basin interior. The horizontal momentum equations are discretized using a third-order upstream scheme (UP3), which has implicit diffusion.

The initial condition is a state at rest with zero velocity and no free surface elevation. The model is used in non-linear and linear configurations. In the case of linear simulations, the non-linear advection terms u • ∇u are discarded from Eq. 1.

c. Wave generation

To study the instability of barotropic Rossby waves and their role in the formation of jet-like structures, we aim to create freely propagating waves using the surface momentum flux described in Eq. 2. Because the direct response to the SMF is not in geostrophic equilibrium, an adjustment occurs by radiating many different types of waves. The objective is to ensure that the wave with the characteristics of the forcing (i.e. with wavevector k 0 and frequency 0 ) is the wave that contains the most of the energy, which we will call hereafter the primary wave. If the wave vector that contains most energy and the forcing frequency 0 satisfy the dispersion relation for barotropic Rossby waves

= - 2 + 2 + (4)
then this wave will be resonantly forced, where ( , ) is the wave vector, the wave frequency, and the inverse of the Rossby deformation radius squared. In this study, we focus on the barotropic mode, = (1/ 2 ,0 ), with ,0 the barotropic Rossby deformation radius. In the equatorial regions, ,0 is very large (i.e.

,0 ∼ 4500km at 20 • , the poleward limit of the domain we are using), thus the upper limit of in this region is ∼ 10 -14 m -2 , which is smaller than the wavenumbers considered (typically |k| 2 ∼ 10 -10 m -2 ). In the following, we thus consider ∼ 0.

Because the SMF wave maker is modulated by the envelope, its spatial spectral footprint is spread over a range of wavenumber (Fig. 2). Its energy maximum in wavenumber space is however found close to k 0 , and intersects the Rossby wave dispersion relation for the frequency 0 , thus satisfying the resonant condition. It is ensured that all the experiments conducted in this study (Table 1) satisfy this condition.

In the case of resonant SMF, the primary wave in the simulations is a barotropic wave with the characteristics of the forcing function (i.e. wavevector k 0 and frequency 0 ) that propagates along the ray paths defined by the integration of the group velocity vector ( , ).

= ( 2 -2 ) ( 2 + 2 ) 2 (5a) = 2 ( 2 + 2 ) 2 (5b)
Because the initial state of the simulations is at rest, there is a spin up time which corresponds to the time it takes for the forcing to transfer energy to the oceanic response and balance the different dissipation terms. The spin up time and the amplitude of the response can vary for the different experiments. It is important to evaluate both quantities to get the best possible description of the generated primary wave, since the characteristics of the primary wave determine the properties of the secondary waves (Sections 5, 6).

We define the spin up time for a wave forced at k 0 and 0 (Eq. 2) as k 0 , 0 the time at which the energy reaches 60% of its maximal value.

k 0 , 0 = min ∈ KE( ) -0.6 × KE max > 0 ( 6
)
where KE is the spatial average of the kinetic energy over the spatial domain of the simulation (Ω), is the total duration of the simulation. The spin up time decreases slightly towards higher forcing frequencies but the values remain around 100-150 days (Fig. 3a).

We define the amplitude for a wave forced at k 0 and 0 (Eq. 2) as k, , the average of the highest velocities values in the forcing region.

k 0 , 0 = 1 Ω k 0 , 0 ∫ k 0 , 0 0 ∬ Ω k 0 , 0 ( , , ) ( k 0 , 0 ( , , ) - * k 0 , 0 ) (7)
where k 0 , 0 is the meridional velocity in a simulation forced by k 0 and 0 (Eq.2), * k 0 , 0 is the 95 th percentile of the distribution of in the forcing region ( ) ( ) (Eq. 3), is the Heaviside function and Ω the simulation domain. The amplitude of the forced wave varies linearly with the magnitude of the SMF (Fig. 3b).

d. Experimental Design

Table 1 summarizes all the experiments used and described in this study. The simulations are chosen to span the intra-annual period range from 30 to 130 days with primary waves that follow the dispersion relation for barotropic Rossby waves with zero meridional wavenumber (experiments a to l in Fig. 1). These = 0 primary waves are found to be the optimal waves to produce jet-like structures, as will be shown in Section 6b. For this reason, our simulations focus on this particular wave type. All experiments (a to l) are run for a nominal SMF amplitude of 0 / 0 = 5 × 10 -5 m 2 s -2 (simulations 1). Experiments a, c, e, g and i are run with for 0 / 0 = 5 × 10 -6 m 2 s -2 (simulations 2), 0 / 0 = 2.5 × 10 -5 m 2 s -2 (simulations 3), 0 / 0 = 7.5 × 10 -5 m 2 s -2 (simulations 4), 0 / 0 = 1.0 × 10 -4 m 2 s -2 (simulations 5) and 0 / 0 = 1.5 × 10 -4 m 2 s -2 (simulations 6). All simulations are run with a non-linear physics (NL) and simulations 1a and 1l are repeated using a linear physics (L). The simulations are run for 1800 days (5 years) with output every 10 days.

Spectral Analysis Method

One objective of the analysis of the simulations is to quantify the primary wave destabilization and transfer of energy between wave modes by NLTI. Spectral analyses are the ideal tool for this task. As the properties of the waves and their spectral content are expected to change in space and time, we compute local spectra using wavelet transforms [START_REF] Graps | An introduction to wavelets[END_REF][START_REF] Gargour | A short introduction to wavelets and their applications[END_REF]).

a. Wavelet transform

A continuous wavelet transform using a complex Morlet wavelet [START_REF] Lee | Pywavelets: A python package for wavelet analysis[END_REF]) is applied to the velocity field ( , ). This method results in an estimation of the amplitude of the velocity signal contained in each frequency ( ( ; , , )) and wavenumber ( ( ; , , ), ( ; , , )) at each location ( , ) ∈ Ω of the domain Ω and each time ∈ of the simulation, over the time interval .

We define the spatialy-averaged one-dimensional spectral kinetic energy as:

KE x ( x ; , , ) = 1 Ω ∬ Ω | x ( x ; ; ; )| 2 + | x ( x ; , , )| 2 (8)
where x ∈ { , , } and is the spectral variable along x dimension (including frequency). Note that this definition of the one-dimensional spectral kinetic energy differs from the wavelet transform of the kinetic energy (( 2 + 2 )/2), and has the advantage of retaining the spectral scales of the velocity field. The spectral dimensions and their resolution are chosen as follows: we estimate the temporal wavelet transform for periods from 10 days to 350 days with a 3.3 days interval, and the spatial wavelet transform for wavelengths from 90 km to 3700 km with a 18.5 km interval. The wavelet transform amplitude cannot be accurately estimated for locations close to the boundaries of the simulation domain and for times close to the beginning and the end of the simulations (especially for long wavelengths and periods). This is known as the edge cone of influence [START_REF] Torrence | A practical guide to wavelet analysis[END_REF]. Estimates that lie within the cone of influence are thus discarded.

b. Two-dimensional spectra from wavelet transform

A wavelet analysis provides a one-dimensional spectrum as the spectral transform is computed on each physical dimension separately. In order to correlate this information, we construct crossproduct spectra x x of a quantity along the spectral dimension x x as defined by Eq. 9.

x x ( x , x ;

) = 1 Ω ∬ Ω | x ( x ; , , )| 2 | x ( x ; , , )| 2 (9)
where x ,x ∈ { , , } and is the spectral variable along the x dimension. Note that x x is equivalent to the two-dimensional Fourier transform, but also varies in time.

We define the two-dimensional spectral kinetic energy as:

KE x x ( x , x ; ) = x x ( x , x ; ) + x x ( x , x ; ) (10)

Numerical solutions a. Propagation of the forced waves in the linear simulations

The linear simulations illustrate the key features of each experiment and help to ensure that the dominant signal generated by the SMF is associated with the forced primary wave (Section c).

In simulations 1d and 1l (Fig. 4), waves with a strong signature in meridional velocity propagate eastward of the forcing region. These waves have the characteristics (i.e. wavenumber and frequency) of the forcing (not shown). Their energy propagation follows the theoretical ray paths predicted from the dispersion relation. The propagation is much faster for high-frequency, long-wavelength forcing (Fig. 4d-f) than for low-frequency, short-wavelength forcing (Fig. 4ac). Note that there is also a slight spreading of the energy in the meridional direction when the wave propagates eastward. Barotropic waves are plane waves and are thus invariant in the direction perpendicular to their propagation direction (the meridional direction in the case of the two examples). The relaxation of the forced solution (which is not an exact plane wave due to the non-zero associated with the envelope ( ) ( )) towards a free plane wave solution outside of the forcing region can thus explain the meridional spread of the wave signature. These linear simulations confirm that the SMF excites principally a barotropic Rossby wave with the forcing characteristics, as expected from the chosen idealized setup (Section b).

b. Propagation of the forced wave and secondary waves in non-linear simulations

When the non-linear terms are taken into account in the model (Eq. 1), the key features of some simulations is modified. For example, in simulation 1d, in addition to the eastward propagation of energy associated with the forced primary wave (Section a), there is also westward energy propagation (Fig. 5a-i). This westward propagation is associated with a wave with quasi-zonal phase lines and has a dominant signature in the zonal velocity and sea surface height. This suggests that an additional wave (aside from the forced wave) is also present in the simulation. Conversely, for other simulations forced by different wavenumbers and frequencies, such as simulation 1l, the non-linear response does not differ much from the linear response (Fig. 5j-r and Fig. 4d-f).

c. Resulting jet-like mean flows in the non-linear simulations

The secondary wave that develops in simulation 1d generates zonal velocities that are coherent over a large zonal extension (i.e. from 0 • to 60 • ) and that change signs over short meridional scales (i.e. every 2.5 • ) (Fig. 5f). When temporally-averaged, this wave results in a jet-like mean flow (Fig. 6). By jet-like mean flow, we specify a circulation that shares common characteristics with the observed jet-structured circulation in the equatorial and tropical oceans (e.g., [START_REF] Ollitrault | Zonal intermediate currents in the equatorial atlantic ocean[END_REF][START_REF] Cravatte | Intermediate zonal jets in the tropical pacific ocean observed by argo floats[END_REF][START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF], that has meridionally-alternating zonal velocities with amplitudes around 0.03 to 0.05 m s -1 , sea surface height anomalies of around 2 to 3 cm that alternate sign on short meridional scales (typically 1 to 2 • ) and that have zonal and temporal coherence over thousand kilometers and several years respectively (see also Section 1).

The average over the last three years of the simulations shows that simulation 1d develops jet-like mean flow (Fig. 6a), even if its characteristics are not fully realistic (Fig. 6c); namely the meridional scale (2 • ) is slightly larger than what is observed in the Pacific (1.5 • ) and the amplitude (0.01 m s -1 ) is weaker than in the observations (0.06 m s -1 ). In contrast, simulation 1l does not develop a jet-like mean flow (Fig. 6b). The generation of jet-like mean flows and their characteristics (meridional scale, amplitude) thus strongly depend on the characteristics of the forced waves.

The aim of the remainder of this paper is to understand the origin of the secondary wave that generates jet-like mean flows, and to understand in what conditions they develop and what are their characteristics. We investigate the physical mechanisms in simulation 1d in depth before examining the sensitivity of this mechanism to the primary wave characteristics using multiple experiments (Section 2d).

Evidence for non-linear triadic instability in the numerical simulations

In this section, we focus on one particular experiment, simulation 1d (Table 1), to describe the secondary wave and explain it in the framework of a 3-mode truncation NLTI model. An analytical solution can be derived for the evolution of the primary wave streamfunction, by limiting its non-linear interactions to a single triad, with three waves exchanging energy. The growth rate ( ) of this instability, which determines the two secondary waves with wavevectors q and pq that interact with the primary wave with wavevector p in a triad (p = q + pq), is given by Eq. 11

and is derived analytically in Appendix A, following [START_REF] Connaughton | Modulational instability of rossby and drift waves and generation of zonal jets[END_REF],

= Im 1 2 (Δ ± Δ 2 -4| p | 2 (q, p, q -p) (p -q, p, -q)) , (11)
where Δ = pqp-q is the resonant triad condition (the triad is resonant if

Δ = 0), (k, k i , k j ) = - (k i ×k j ) ( 2 -2 ) 2
and Im designates the imaginary part. The growth rate depends on the wavevector (p), frequency ( p ) and amplitude of the streamfunction (| p |) of the primary wave. Note that the triad instability described by this model is more general than the resonant triad instability (e.g., [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF]). The approach is very similar to that of [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF]. But while [START_REF] Gill | The stability of planetary waves on an infinite beta-plane[END_REF] explored the weakly non-linear limit | p | 3 / 1 (leading to triad resonance) and the strongly non-linear limit | p || 3 / 1 as two separate instabilities, the analytical solution used here is a continuous function of the primary wave amplitude | p |. Given the characteristics of the primary wave, Eq. 11 allows one to determine the secondary waves wavevectors, q and pq, that yield the maximum growth rate. These waves are the ones most likely to emerge from non-linear interactions.

a. Theoretical prediction for instability growth rate and secondary waves

We compute the growth rate given by Eq. 11 for the primary wave p of simulation 1d (Fig. 7).

The secondary waves q and pq that will develop as predicted by theory are determined by the maxima of this growth rate. The growth rate maximum is found numerically using a wavevector grid spanning -5 × 10 -5 to -5 × 10 -5 m -1 in both the and directions with a resolution of 10 -7 m -1 . The two secondary waves resulting in the fastest growth have different characteristics.

The first one has a long zonal wavelength (∼5000 km), a short meridional wavelength (∼430km) (Fig. 7a), and a long period (∼1.5 years) (Fig. 7b), called hereafter the long secondary wave. The second one has a short zonal wavelength (∼340 km, similar to the one of the primary wave), a short meridional wavelength (∼430 km) (Fig. 7a), and a short period (∼4 months) (Fig. 7b), called hereafter the short secondary wave. The exact wavenumbers and frequencies associated with each of these waves are listed in Table 2. The long secondary wave thus has characteristics similar to the jet-like structures.

The growth rate of the secondary waves is max = 0.032 day -1 (Fig. 7). The primary wave is very unstable to NLTI because -1 max is of the same order of magnitude as its period. We define the time of emergence of the instability as:

p = 5 max ; such that q ( p ) = 5 q (0) 10 2 q (0) ( 12
)
where max is the maximum growth rate of the instability (Eq. 11). p corresponds to the time it takes for the amplitude of the secondary wave q (the fastest growing perturbation of the primary wave p), which has initially a small amplitude q (0), to increase by a factor 5 10 2 since the secondary waves grow as max (Appendix A). The factor of 5 that has been chosen will be justified a posteriori. We thus expect to detect the secondary waves after 5 -1 max 150 days. The group velocity, which determines the direction of energy propagation, is also different for the three waves (Fig. 8). It is eastward for the primary wave, westward for the long secondary wave and poleward for the short secondary wave. We thus expect to detect the signature of the primary wave eastward of the forcing region, the signature of the long secondary wave westward of the destabilization region (which corresponds approximately to the forcing region in this case since the primary wave has a slow group velocity and the secondary wave has a fast growth rate), and the signature of the short secondary wave poleward of the forcing region.

b. Comparison with the spectral analysis of the non-linear simulation

We perform spatial and temporal wavelet analyses of the velocity in the simulations to infer the spectral content at the time and location where we expect the primary and secondary waves to be strongest.

The primary wave is observable after a few wave periods east of the forcing region (see ray tracing Fig. 4 and group velocity Fig. 8). The wavelet transform in this region (red square in Fig. 5c) shows a spectral peak at a period of 79 days (Fig. 9a; solid red curve), a zonal wavelength of 250 km (Fig. 9b; solid red curve) and a long meridional wavelength of 3000 km (Fig. 9b; dashed red curve). This matches the characteristics of the expected forced primary wave (Table 2). The long secondary wave is detected west of the forcing region, the spectral characteristics measured from the wavelet transforms in this region (green square in Fig. 5f) show a long period and a long zonal wavelength (not precisely quantifiable within the spanned spectral range (Fig. 9; solid green curve)). The meridional wavelength is around 500 km (Fig. 9b; dashed green curve). This wave has spectral characteristics consistent with the prediction from the NLTI theory (Table 2).

The short secondary wave is difficult to identify in the simulations because its period and zonal wavenumber are similar to the primary wave's (Table 2). It also has a slower group velocity than the long secondary wave (Fig. 8), and its energy propagation is directed poleward, in a region where the primary wave is also present. There is thus no region where the short secondary wave can be observed alone. Its presence can however be identified in the two-dimensional spatial spectrum (Fig. 11c).

The time at which the instability occurs in this simulation is inferred from the wavelet scalogram, which shows the evolution of the mean spectral characteristics averaged in the basin KE ( ), KE ( ), KE ( ) (Fig. 10). The instability occurs when there is a shift in the spectral characteristics from the forced frequency and wavelengths to the emerging secondary wave, which is around = 275 days in the simulation (dashed white line in Fig. 10). The shift is particularly visible for KE ( ) (Fig. 10c) where the energy goes from a large to small meridional wavelength. The energy in the forced, small period and small zonal wavelength is concurrently strongly attenuated (Fig. 10a,b). The theoretical growth rate corresponds to an instability emergence time of ∼ 150 days (Section a). Based on Fig. 10, the spin up period, which corresponds to the time before the detection of the spectral characteristics of the primary wave, is around 130 days. The instability emergence time in the simulation (275 days) corresponds thus approximately to the sum of the spin up period and theoretical instability emergence time. The spatially-averaged two-dimensional spectral kinetic energy shows that basin wide, the velocity field is dominated by the forced wave before the instability occurs (Fig. 11 a,b) and by the two secondary waves after the instability occurs (Fig. 11 c,d). Note that the two-dimensional spectrum KE allows one to separate the contribution of the short secondary wave (Fig. 11 c), and reveals that it is less energetic than the long secondary wave that is dominant in both KE and KE .

We have thus shown that the secondary waves that develop in the non-linear simulations are consistent with a non-linear triad instability. The growth rate and wavenumbers of the instability are well predicted using a non-linear truncated theory. In the example shown, the long secondary wave has jet-like characteristics.

Application to the generation of jet-like structures

The objective of this section is to explore the sensitivity of the formation of jet-like long secondary waves to the characteristics of the primary wave (i.e. the frequency, wavevector and amplitude).

From the simulations that have been performed, it appears that not all primary waves are able to trigger instability, nor to create waves that resemble jet-like structures.

a. Instability threshold and emergence time

Some non-linear simulations are not subject to instability, like for example simulation 1l, in which the forced primary wave propagates without spectral modification (Fig. 5j-r). Theoretically, all primary waves are subject to triad instability (Appendix A). But the growth rate can be close to zero, in which case it would take a very long time for the instability to develop, and thus it might not be realizable in numerical simulations. To determine whether a wave will be readily unstable or not, we define an instability emergence ratio (Eq. 13) following [START_REF] Qiu | Generation of the north equatorial undercurrent jets by triad baroclinic rossby wave interactions[END_REF], which is a combination of the instability growth rate and the primary wave energy propagation velocity. We define the instability emergence ratio as: p corresponds thus to the ratio of the propagation time scale to the instability emergence time.

The instability emergence time depends on the primary wavenumber, frequency and amplitude (Fig. 12a). As predicted by the theoretical instability emergence ratio (Eq. 13), simulations 1a, 1d and 1f, should undergo triad instability (Fig. 12a). These simulations have a primary wave amplitude around 0.07 m s -1 (see Fig. 3b for the relation between the forcing strength and the primary wave amplitude). The emergence of the instability in these simulations is visible from the switch of spectral properties in KE ( ) (Fig. 12b-d). In contrast, the primary wave in the simulation 1h remains stable and KE ( ) shows a dominant meridional scale that is constant in time (Fig. 12e). This shows that the truncated NLTI theory describes well what happens in the different simulations and can be used to predict the stability or instability of the different barotropic waves in the ocean. It is interesting to note that some waves will not destabilize unless they have an unrealistic amplitude (Fig. 13). This is the case for barotropic waves with a long wavelength.

For example, a wave with a TIW-like spectral signature, with a 30-day period and a 1000-km zonal wavelength destabilizes only for amplitude greater than 75 cm s -1 (Fig. 13). For other waves, and in particular short wavelength primary waves, the amplitude threshold needed to reach p 1 is very low.

The time at which the instability occurs in the simulation can be approximated by the time at which the energy in the forced spectral range reaches a maximum and starts to decrease. For all unstable simulations ( p 1) this time is found to be inversely proportional to the theoretical growth rate expected from the 3-mode truncation triad interaction theory (Fig. 14). Primary waves with a larger amplitude (big dots in Fig. 14) have a higher instability growth rate and thus a lower instability emergence time than primary waves with a smaller amplitude. Similarly, simulations with a higher primary wave number (blue dots in Fig. 14) have a higher growth rate and lower instability emergence time than simulations with lower primary wavenumber (orange dots in Fig. 14). The coefficient of proportionality between the inverse of the growth rate and the emergence time is approximately equal to 5, which justifies the definition of the emergence time (Eq. 12).

b. Selection for the secondary waves' spatial and temporal scales

The non-linear simulations show that the secondary waves resulting from the triad instability differ from one simulation to another, for example, in their meridional wavelengths (Fig. 12b-d).

The characteristics of the long secondary wave depend on the forced primary wave's wavenumbers and frequency, as inferred from the growth rate computation at a fixed amplitude (Fig. 15). The correspondence between the scale for the long secondary wave in the numerical simulations and the theoretical prediction in Fig. 15 shows that the dominant instabilities in the simulations are consistent with NLTI theory.

The long secondary wave meridional wavelength varies with the period of the primary wave (2 / p ). It is shorter (∼ 100-300 km) for long periods (2 / p > 75 days) and much longer (up to 1000-3000 km) for shorter periods (2 / p < 50 days). At fixed p , the secondary wave meridional wavelength also increases with decreasing (Fig. 15a). The long secondary wave zonal wavelength varies non-monotonically with along the primary wave's spectral characteristics.

At fixed p , it is large for short and long primary waves. It reaches a maximum for short intra-annual primary waves, with periods between 50 and 145 days, along the dispersion relation curves corresponding to low meridional wavenumbers (Fig. 15b). The long secondary wave's period is longer for short primary waves (large ) and shorter for long primary waves (small ). Similar to the zonal wavelength, the secondary wave's period reaches a maximum for short intra-annual primary waves with periods comprised between 50 and 145 days, along the dispersion relation curves for low meridional wavenumbers (Fig. 15c). The emergence time of the instability (defined by Eq. 12) increases from about a hundred days to several years from high to low respectively. Short primary waves are thus much more unstable. The short intra-annual waves described in the paragraphs above in particular have very short emergence times (Fig. 15d).

Realistic jet-like structures have a meridional wavelength of about 350 km, a long zonal wavelength (greater than 8000 km) and a long period (more than 10 years). Given the dependence of the secondary wave's properties on the primary wave's characteristics (Fig. 15a-c), the "optimal" primary waves that yield realistic jet-like secondary waves have periods between 70 and 90 days, zonal wavelengths from 200 to 300 km, and low meridional wavenumber (as indicated by the squared box in Fig. 15). From the simulations that have been performed, simulation 1c seems to have the most realistic jet-like secondary wave. This is consistent with Fig. 6.

The analysis of the dependence of the scales of the long secondary wave to the primary wave's properties described above assumed a constant primary wave amplitude of 0.07 m s -1 , corresponding to the amplitude of the primary wave used in the numerical simulations. The optimal primary wave spectral characteristics shown in Fig. 15 are fairly robust, regardless of the primary wave's amplitude. Indeed, Fig. 16a shows that the long secondary wave's meridional wavelength is not very sensitive to the amplitude of the primary wave, in particular for the range of simulations that have been performed here. The variation coefficient, defined as the ratio of the standard deviation of the secondary wave's properties across an amplitude range to the mean of these properties, is close to zero for most of primary waves, except for waves with very high meridional wavenumbers (Fig. 16d). The secondary wave zonal wavelength varies with the amplitude of the primary wave (Fig. 16 b,e). These variations are however not very large with respect to the values of the wavelengths: even if they vary by one or two thousands of kilometers, they are still long waves (i.e. with wavelengths larger than 3000 km). The variation coefficient remains less than 0.5 almost everywhere in the spectrum (Fig. 16e). The period of the long secondary wave is the most sensitive to the amplitude of the primary wave. Especially along the first barotropic dispersion relation ( = 0), where the variation coefficient is close to 0.8 (Fig. 16f) and the secondary waves' period can go from 200 days to a few years (Fig. 16c). In the simulations, the secondary wave meridional wavelength is also found to be not sensitive to the primary wave amplitude (Fig. 16a), illustrating again the similarity between the numerical simulations and the NLTI theory in complement of Fig. 15a.

Discussion, Summary and Perspectives

a. Beyond the 3-mode truncation theory

This study shows that NLTI theory can predict the instabilities that can form jet-like structures in the equatorial oceans. The theoretical framework was however developed under strong assumptions, in particular non-linear interactions were limited to three waves. In the model, the dynamics is not truncated and the non-linear interactions may involve additional triads of waves. However, it is interesting to note that the 3-mode truncation captures the early stages of the simulations well and provides a simple description of the instability process. However, on longer time scales, secondary instabilities may arise, which involve additional waves. The spectra show evidence of this (e.g., Fig. 11d) as they indicate that more than 3 waves are clearly present in the simulations. Spectral analysis of the simulation 1d on long time scales shows that energy is present at various scales and is exchanged through a broad continuum of waves, similar to a turbulent cascade, rather than being contained in a few discrete waves (Fig. 17). Energy tends however to accumulate in regions where = 0, which correspond to jet-like structures. This might be explained by the stability of jet-like waves to NLTI. Namely, the growth rate for any jet-like primary wave, p = (0, ), is identically equal to zero (Fig. 18), such that jet-like secondary waves are very stable and will likely not undergo further triad interactions. The jet-like waves thus act as an attractor of this "wave cascade", and every wave undergoing NLTI might eventually cascade towards jet-like structures.

This is analogous to turbulent cascades on a beta plane, where the energy cascades from being isotropically distributed in spectral space towards accumulation zones at low frequencies and small zonal wavenumbers, which correspond to zonal jets [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF][START_REF] Vallis | Generation of mean flows and jets on a beta plane and over topography[END_REF].

b. Beyond barotropic dynamics

Because we are only considering a two-dimensional dynamics, we neglect the possible interactions between different vertical modes. A full description of the non-linear dynamics of equatorial waves should consider baroclinic modes and their possible non-linear interactions. The two-dimensional problem investigated in this paper is well posed for the barotropic mode only.

Non-linear interactions of barotropic waves indeed remain barotropic (not shown). In this case the two-dimensional reduction of the problem is a good approximation and 2D simulations are a good tool to study the NLTI. This is not the case for baroclinic waves. Non-linear interactions of baroclinic modes project mainly on the barotropic mode (not shown). The barotropic mode thus acts also as an attractor in the "wave cascade".

c. Towards a realistic mean circulation

Intra-annual waves with large wavenumbers such as the ones described in this study are not unrealistic as they have been observed in the tropical ocean [START_REF] Farrar | Barotropic rossby waves radiating from tropical instability waves in the pacific ocean[END_REF][START_REF] Tuchen | Deep intraseasonal variability in the central equatorial atlantic[END_REF]Delpech et al. 2020b). [START_REF] Farrar | Barotropic rossby waves radiating from tropical instability waves in the pacific ocean[END_REF] 2011) are found to be stable to NTLI (e.g., Fig. 13). This could explain why they are seen propagating away from their generation site. On the contrary, waves unstable to NLTI have only sporadically been observed, their small wavelengths make them challenging to capture with in-situ data and as they are associated with a weak signature in sea surface height (typically less than 2 cm at low latitudes), it is also difficult to capture them with the current generation Nadir altimetry. The new generation swath satellite altimeter designed to observe submesoscale features (e.g., SWOT Fu and Ubelmann (2014)) may open new opportunities in this regard. The jet-like structures produced by the unstable intra-annual waves through NLTI are long waves with zonal wavelengths that can reach several thousands of kilometers (Fig. 16b) and with periods of several years (Fig. 16c). However, their time averaged velocities generally result in a weak mean circulation (e.g., Fig. 6a). This mechanism is therefore not sufficient to explain the observed mean deep circulation. Other mechanisms may intervene in the long term equilibration of these jet-like structures to produce the jet-structured mean circulation as has been observed in 10-year averages of Argo float displacements [START_REF] Ollitrault | Zonal intermediate currents in the equatorial atlantic ocean[END_REF][START_REF] Cravatte | Intermediate zonal jets in the tropical pacific ocean observed by argo floats[END_REF]. Equilibration processes may involve modification of the potential vorticity by the zonal jets combined with positive feedback mechanisms such as eddy-driven jet-sharpening [START_REF] Dritschel | Jet sharpening by turbulent mixing[END_REF][START_REF] Greatbatch | Evidence for the maintenance of slowly varying equatorial currents by intraseasonal variability[END_REF][START_REF] Ménesguen | Observations and mechanisms for the formation of deep equatorial and tropical circulation[END_REF]). In addition, the observed zonal jets are not purely barotropic and have a complex vertical structure [START_REF] Qiu | Generation of the north equatorial undercurrent jets by triad baroclinic rossby wave interactions[END_REF][START_REF] Cravatte | Subthermocline and intermediate zonal currents in the tropical pacific ocean: Paths and vertical structure[END_REF]. Therefore baroclinic dynamics could come into play and should be considered in future studies. Finally, bottom topography could have an influence on the dynamics of barotropic Rossby waves (such as causing scattering, refraction and trapping) through the topographic beta effect (e.g., [START_REF] Wang | Influence of mid-ocean ridges on rossby waves[END_REF][START_REF] Durland | Model of 33-day barotropic rossby waves in the north pacific[END_REF] and the modification of classical vertical modes [START_REF] Lacasce | The prevalence of oceanic surface modes[END_REF]. The topographic beta effect scales as the Coriolis parameter and the bottom slope. We therefore expect it to be of lesser importance at low latitudes, away from steep continental margins or seamounts. It is however interesting to note that the influence of topography could transform a stable large scale long period Rossby wave forced by the atmosphere, into a small scale long period Rossby wave, subject to NLTI [START_REF] Rhines | Topographic rossby waves in a rough-bottomed ocean[END_REF].

d. Summary and Perspectives

Deep intra-annual waves, which have been frequently observed at low latitudes in the ocean are prone to instability. Focusing on the barotropic mode, we demonstrated that the energy transfer mechanism of this instability is well described by a 3-mode truncation NLTI theory. The growth rate and characteristics of the instability are a function of the wavenumber, frequency and amplitude of the primary wave. In particular cases, the secondary waves have jet-like structures. In this study, we showed that intra-annual waves can destabilize into jet-like structures using idealized barotropic numerical simulations of an equatorial basin forced by an oscillatory wind-stress. We determined that the primary waves that are optimal to reproduce realistic jet-like structures are intra-annual zonally-propagating waves with period between 60 to 130 days, zonal wavelengths between 200 to 300 km and amplitude from 2 to 10 cm s -1 . The jet-like structures produced in the simulations and predicted by the theory do not capture all the properties of the time-mean deep circulation as it has been observed in low-latitude oceans (in particular the jets are too weak). Future work should address the potential role of three-dimensional dynamics in the generation of these jets and include baroclinic modes into NLTI theory. The effects of bottom topography and variable stratification should be investigated as they can modify the properties and propagation of Rossby waves. Future work will also involve running realistic numerical simulations to study the role of NLTI and other equilibration processes in the formation of the jet-structured mean circulation within a realistic equatorial wave field.

T 1. Experiment Table . Summary of all numerical experiments used for this study, where 0 , 0 , 0 , 0 , and 0 are the period, frequency, zonal wavelength, zonal wavenumber and amplitude of the forcing (Eq. 2).

The experiments are 2D and are either linear (L) or non-linear (NL). The letter of the experiment label indicates its position on the wavenumber-frequency spectrum (Fig. 1), and the number indicates the forcing amplitude. T 2. The waves involved in non-linear triad interactions. The primary wave is the forced wave, p , , thus are fixed by the forcing. As the experiment are resonantly forced for the barotropic mode, they satisfy the barotropic Rossby wave dispersion relation. The properties of the two secondary waves (q and p -= pq)

are determined from where the growth rate is maximum in spectral space (Eq. 11).

Forced wave (primary wave) Long secondary wave Short secondary wave

(m -1 ) = -2.3 × 10 -5 = -2.1 × 10 -6 -= -2.0 × 10 -5 (km) p = 250 km q = 4958 km p -= 238 km (m -1 ) = 0 = 1.4 × 10 -5 -= 1.7 × 10 -5 (km) p ∼ ∞ q = 429 km p -= 429 km (s -1 ) p = 9.8 × 10 -7 q = 2.6 × 10 -7 p -= 7.2 × 10 -7
T (days) p = 79 days q = 540 days p -= 110 days and (j-r) 1l, at times: (a,d,g,j) t = 100 days, (b,e,h,k) t = 400 days and (c,f,i,l) t = 800 days for (a-c, j-l) (meridional velocity); (d-f, m-o) (zonal velocity); (g-i, p-r). The red and green squares correspond to the locations where the one-dimensional spectra are computed (Fig. 9). . 49 Fig. 6. averaged over years 2 to 5 of the simulations from 10 • to 40 • of longitude (i.e. west of the forcing region) for simulations (a) 1d (b) 1l (c) and simulation with optimal primary wave characteristics, corresponding to experiment 1c from Table 1, but with a shifted wind-stress position ( 0 =3 • S) (d) Observed zonal velocity at 1000 m depth in the tropical Pacific from Argo floats YoMaHA database [START_REF] Lebedev | Velocity data assessed from trajectories of argo floats at parking level and at the sea surface[END_REF]) reproduced after (Delpech et al. 2020b). . 50 is expected (red marker in Fig. 5c) and (green curves) Spectrum of the zonal velocity at the location where the long secondary wave q is expected (green marker in Fig. 5f). , (b) zonal wavelength , (c) period, to the amplitude of the primary wave. The sensitivity is computed for primary waves with zonal wavelengths 150 km, 200 km, 300 km, and 400 km and for amplitudes ranging from 0.01 m s -1 to 0.3 m s -1 with a 0.01 m s -1 increment.

LIST OF FIGURES

The corresponding experiments and the secondary waves' meridional wavelengths computed using a wavelet analysis are indicated by the stars with their labels referring to Table 1. (lower panel) The variation coefficient, namely the ratio of the standard deviation of the properties for the different amplitude to the mean of this propreties are computed for the secondary wave meridional wavelength (d), zonal wavelength (e), and period (f). . . . . . . 60 (grey) for simulation 1d (see Table 1). Black straight lines correspond to the forced wavenumber ± 0 and 0 = 0. Black circles correspond to the barotropic Rossby wave dispersion relation for the forced frequency 0 (T=79 days) and a non-resonant frequency (T=180 days). If the dispersion relation intersects a location in wavenumbers space where there is a maximum in energy, the forced wave is resonantly forced. The example shows that only the period corresponding to the barotropic Rossby wave dispersion relation for 0 and 0 satisfies the resonant condition. F . 3. (a) Experiment spin-up time (computed as in Eq. 6) as a function of the forcing frequency ( 0 ).

(b) Amplitude of the wave response (computed as in Eq. 7) as a function of the strength of the forcing ( 0 ) for experiments a-l (Fig. 1). The best linear fit (line) corresponds to a line = + and is the correlation coefficient. F . 6. averaged over years 2 to 5 of the simulations from 10 • to 40 • of longitude (i.e. west of the forcing region) for simulations (a) 1d (b) 1l (c) and simulation with optimal primary wave characteristics, corresponding to experiment 1c from Table 1, but with a shifted wind-stress position ( 0 =3 • S) (d) Observed zonal velocity at 1000 m depth in the tropical Pacific from Argo floats YoMaHA database [START_REF] Lebedev | Velocity data assessed from trajectories of argo floats at parking level and at the sea surface[END_REF]) reproduced after (Delpech et al. 2020b). zonal wavelength , (c) period, to the amplitude of the primary wave. The sensitivity is computed for primary waves with zonal wavelengths 150 km, 200 km, 300 km, and 400 km and for amplitudes ranging from 0.01 m s -1 to 0.3 m s -1 with a 0.01 m s -1 increment. The corresponding experiments and the secondary waves' meridional wavelengths computed using a wavelet analysis are indicated by the stars with their labels referring to Table 1. (lower panel) The variation coefficient, namely the ratio of the standard deviation of the properties for the different amplitude to the mean of this propreties are computed for the secondary wave meridional wavelength (d), zonal wavelength (e), and period (f). and always have a growth rate equals to zero, illustrating their particular stability.

  |k| is the modulus of the wavevector (|k| = 2 + 2 ) and |C gk | the modulus of the wave group velocity vector (|C gk | = 2 + 2 ). Following Qiu et al. (2013), p represents the propagation time scale of the primary wave p (i.e. the time needed for a wavepacket to transit of one wavelength).

  evidenced barotropic intra-annual waves in the northern tropical Pacific ocean, radiating away from the Tropical Instablity Waves with a 30-day period and a 1000-km wavelength. Tuchen et al. (2018) described evidence of 50-day Yanai waves with a baroclinic structure at the equator. Delpech et al. (2020b) described observations of waves at 1000 m in the northern tropical Pacific with a 70-day period and a 500-km wavelength with properties consistent with barotropic Rossby waves. It is interesting to note that the waves observed by Farrar (

Fig. 1 .Fig. 2 .Fig. 4 .Fig. 5 .

 1245 Fig. 1. Wavenumber frequency diagram with the dispersion relations for equatorial Kelvin waves (green), Yanai waves (red), and baroclinic Rossby waves (blue). Only the first baroclinic mode is represented. The dispersion relation for barotropic Rossby waves with different meridional wavenumbers (0, 1 × 10 -6 , 2 × 10 -6 , 4 × 10 -6 m -1 from top to bottom) are in black. The frequency and wavenumbers of the forcing used in the experiments conducted in this study are represented by the yellow dots. Wave instability simulations conducted in previous studies are indicated with colored dots along with their references. The square boxes represent the spectral range of waves that have been observed at the surface (solid lines) and at depth (dashed lines). . . . . . . . . . . . . . . . . . 45 Fig. 2. Two-dimensional spectrum of the wind stress function ( , , ) = ( ) ( ) ( 0 + 0 -0 ) (grey) for simulation 1d (see Table 1). Black straight lines correspond to the forced wavenumber ± 0 and 0 = 0. Black circles correspond to the barotropic Rossby wave dispersion relation for the forced frequency 0 (T=79 days) and a non-resonant frequency (T=180 days). If the dispersion relation intersects a location in wavenumbers space where there is a maximum in energy, the forced wave is resonantly forced. The example shows that only the period corresponding to the barotropic Rossby wave dispersion relation for 0 and 0 satisfies the resonant condition. . . . . . . . . . . . . . 46 Fig. 3. (a) Experiment spin-up time (computed as in Eq. 6) as a function of the forcing frequency ( 0 ). (b) Amplitude of the wave response (computed as in Eq. 7) as a function of the strength of the forcing ( 0 ) for experiments a-l (Fig. 1). The best linear fit (line) corresponds to a line = + and is the correlation coefficient. . . . . . . . . . . . . . . 47 Fig. 4. Evolution of the meridional barotropic velocity ( ) in the linear simulations (a-c) 1d-L and (d-f) 1l-L at times (a,d) t = 100 days, (b,e) t = 400 days, (c,f) t = 800 days. Black lines indicate the ray path originating from 70 • longitude (the center of the forcing region) and -5 • , 0 • and 5 • latitude respectively and integrated from the initial time of the simulation to the time of the snapshot (Eq. 5). Position of the ray at each wave period is indicated by the dots. The zonal barotropic velocity ( ) is almost zero in the whole basin and is not shown. . . 48 Fig. 5. Evolution of the velocity and sea surface height (ssh) in the non-linear simulations (a-i) 1d

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. The growth rate for the 3-mode truncated non-linear triad interaction (Eq. 11) (a) as a function of the secondary wavenumbers and . The red arrow indicates the primary wave p and the green arrows the most unstable secondary waves q and pq. The circles represent the magnitude of the wavevectors. (b) The growth rate as a function of the secondary wavenumber and period q . The red dot indicates the position of the primary wave and the green dots the position of the secondary waves inspace. Black lines indicate the dispersion relations for barotropic Rossby waves with = 0 and = 2 × 10 -6 m -1 . The two secondary waves have the same growth rate. Note that we use the convention of positive

Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .Fig. 15 .Fig. 16 .

 10111213141516 Fig. 10. Wavelet scalogram averaged in the interior of the domain (i.e the region away from 5 • of the boundary) as a function of (a) time and period KE (Eq. 8); (b) time and zonal wavelength KE ; (c) time and meridional wavelength KE . The spin up period is indicated by the arrow. The time when the instability occurs is indicated with the dashed line, and corresponds to a shift in the spectra. . . . . . . . . . . . . . . . . . . . . . 54 Fig. 11. Two-dimensional spectra computed from the wavelet transforms following the method described in Section 3. (a,c) KE ; (b,d) KE ; (a,b) before the instability emerges ( = 150 days); and (c,d) after the instability emerges ( = 680 days). The red (green) dots and arrows indicate the theoretical frequency and wavenumbers of the primary (secondary) waves respectively. White lines (b,d) indicate the dispersion relation for barotropic Rossby waves with = 0 and = 2 × 10 -6 . . . . . . . . . . . . . . . . . . 55 Fig. 12. (a) Theoretical instability emergence ratio as a function of the amplitude of the primary wave p for four waves with periods and zonal wavelengths indicated in the legend. The instability emergence ratio (Eq. 13) is defined as the ratio between the time it takes a wavepacket to travel one wavelength and the emergence time of the secondary wave, illustrating how the instability is triggered in the simulations when this ratio exceeds 1 (i.e. for times to the right of the dashed line). (right panels) KE ( ) for simulations (b) 1a; (c) 1d; (d) 1f and (e) 1h. . . 56 Fig. 13. Theoretical amplitude threshold to reach p = 1 (Eq. 13) as a function of the primary wave zonal wavenumber ( ) and frequency . Above this amplitude, barotropic Rossby primary waves are unstable to NLTI. . . . . . . . . . . . . . . . . . . . 57 Fig. 14. The emergence time of the instability in the simulations as a function of the theoretical growth rate for the primary waves of the different simulations. Each dot represents one simulation, with size and color indicating the amplitude and zonal wavenumber of the primary wave respectively. The black curve is proportional to the inverse of the growth rate, with a proportionality constant of 5.5. . . . . . . . . . . . . . . . . . . 58 Fig. 15. Properties of the long secondary wave as a function of the primary wave's wavelength and period with a fixed amplitude of 0.07 m/s. (a) , (b) , (c) period, (d) time of emergence of the secondary wave. The dispersion relation for barotropic Rossby waves with different is indicated by solid lines. The hatched regions correspond to regions where no instabilities grow. . . . . . . . . . . . . . . . . . . . . . . . . . 59

Fig. 17 .Fig. 18 .

 1718 Fig. 17. Evolution of the energy spectrum in wavenumber space for simulation 1d at times (a) = 120 days, (b) = 1000 days, (c) = 1200 days. . . . . . . . . . . . . . . . 61 Fig. 18. The maximum growth rate (Eq. 11) for | | ranging between 0 and 5e-5 m -1 for different meridional primary wavenumbers (line styles) and amplitude 0 | | (line widths). Jetlike waves correspond to = 0 and always have a growth rate equals to zero, illustrating their particular stability. . . . . . . . . . . . . . . . . . . . . 62

  Two-dimensional spectrum of the wind stress function ( , , ) = ( ) ( ) ( 0 + 0 -0 )

F . 4 .

 4 Evolution of the meridional barotropic velocity ( ) in the linear simulations (a-c) 1d-L and (d-f) 1l-L at times (a,d) t = 100 days, (b,e) t = 400 days, (c,f) t = 800 days. Black lines indicate the ray path originating from 70 • longitude (the center of the forcing region) and -5 • , 0 • and 5 • latitude respectively and integrated from the initial time of the simulation to the time of the snapshot (Eq. 5). Position of the ray at each wave period is indicated by the dots. The zonal barotropic velocity ( ) is almost zero in the whole basin and is not shown. F . 5. Evolution of the velocity and sea surface height (ssh) in the non-linear simulations (a-i) 1d and (j-r) 1l, at times: (a,d,g,j) t = 100 days, (b,e,h,k) t = 400 days and (c,f,i,l) t = 800 days for (a-c, j-l) (meridional velocity); (d-f, m-o) (zonal velocity); (g-i, p-r). The red and green squares correspond to the locations where the one-dimensional spectra are computed (Fig. 9).

  The growth rate for the 3-mode truncated non-linear triad interaction (Eq. 11) (a) as a function of the secondary wavenumbers and . The red arrow indicates the primary wave p and the green arrows the most unstable secondary waves q and pq. The circles represent the magnitude of the wavevectors. (b) The growth rate as a function of the secondary wavenumber and period q . The red dot indicates the position of the primary wave and the green dots the position of the secondary waves inspace. Black lines indicate the dispersion relations for barotropic Rossby waves with = 0 and = 2 × 10 -6 m -1 . The two secondary waves have the same growth rate. Note that we use the convention of positive frequency for waves, such that the wavevector with > 0 (and thus a theoretical negative frequency) in panel (a) is represented with a . Wavelet scalogram averaged in the interior of the domain (i.e the region away from 5 • of the boundary) as a function of (a) time and period KE (Eq. 8); (b) time and zonal wavelength KE ; (c) time and meridional wavelength KE . The spin up period is indicated by the arrow. The time when the instability occurs is indicated with the dashed line, and corresponds to a shift in the spectra. . Two-dimensional spectra computed from the wavelet transforms following the method described in Section 3. (a,c) KE ; (b,d) KE ; (a,b) before the instability emerges ( = 150 days); and (c,d) after the instability emerges ( = 680 days). The red (green) dots and arrows indicate the theoretical frequency and wavenumbers of the primary (secondary) waves respectively. White lines (b,d) indicate the dispersion relation for . (a) Theoretical instability emergence ratio as a function of the amplitude of the primary wave p for four waves with periods and zonal wavelengths indicated in the legend. The instability emergence ratio (Eq. 13) is defined as the ratio between the time it takes a wavepacket to travel one wavelength and the emergence time of the secondary wave, illustrating how the instability is triggered in the simulations when this ratio exceeds 1 (i.e. for times to the right of the dashed line). (right panels) KE ( ) for simulations (b) 1a; (c) 1d; (d) 1f and (e. Theoretical amplitude threshold to reach p = 1 (Eq. 13) as a function of the primary wave zonal wavenumber ( ) and frequency . Above this amplitude, barotropic Rossby primary waves are unstable . The emergence time of the instability in the simulations as a function of the theoretical growth rate for the primary waves of the different simulations. Each dot represents one simulation, with size and color indicating the amplitude and zonal wavenumber of the primary wave respectively. The black curve is proportional to the inverse of the growth rate, with a proportionality constant of 5. Properties of the long secondary wave as a function of the primary wave's wavelength and period with a fixed amplitude of 0.07 m/s. (a) , (b) , (c) period, (d) time of emergence of the secondary wave. The dispersion relation for barotropic Rossby waves with different is indicated by solid lines. The hatched regions correspond to regions where no instabilities grow. . (upper panel) Sensitivity of the long secondary waves' properties: (a) meridional wavelength , (b)

F

  . Evolution of the energy spectrum in wavenumber space for simulation 1d at times (a) = 120 days, (b) = 1000 days, (c) = 1200 days. . 18. The maximum growth rate (Eq. 11) for | | ranging between 0 and 5e-5 m -1 for different meridional primary wavenumbers (line styles) and amplitude 0 | | (line widths). Jet-like waves correspond to = 0

Acknowledgments. The authors aknowledge the support of the CROCO community. CROCO is provided by http://www.croco-ocean.org. This work was granted access to the HPC facilities DATARMOR of the Pole de Calcul Intnsif pour la Mer at Ifremer, Brest (France). The authors are grateful to N. Ducousso and R. Benschila for their advices on the numerical code. This work benefited from fruitful discussions with C. Connaughton, E. Simmonet, B. Qiu, B. Delorme and L. Qu. The authors are also grateful to two anonymous reviewers, whose comments helped to improve the quality of the manuscript. This publication is supported by the CNRS/LEFE and the CNES/TOSCA Alti-ETAO and Saral-MANCA projects. APPENDIX

Analytical solutions for waves non-linear triad interactions

The aim of this section is to understand the dynamics of the forced wave and in particular its stability. We describe its stability (or instability) using a theoretical truncated non-linear model which allows to resolve analytically the first order of non linearity. We derive this problem in a quasi-geostrophic framework, which corresponds to the limit of the primitive equation model for small amplitudes (a reasonable assumption when working with waves).

The two dimensional quasi-geostrophic equation for the stream function is:

where is the Jacobian operator and is the inverse of the squared deformation radius. The deformation radius are defined as ,0 for the barotropic deformation radius and , for the th 29 baroclinic deformation radius.

A plane monochromatic wave, = 0 ( + -) , is solution of the equations if it satisfies the relation.

which is known as the Rossby wave dispersion relation. Plane monochromatic Rossby waves are thus exact solutions of the non-linear quasi-geostrophic equation (Eq. A1).

However, a superposition of waves, = = 0 ( ) ( ) , is not an exact solution of the Eq. A1.

+

( -

Because of the quadratic nature of non linear Jacobian term in Eq. A1, the waves interaction is triadic: the waves i, j and n interact together and their amplitude 0 depend on time.

The equation for the evolution of a wave of streamfunction k and wavenumber k is then given by:

where

is the triadic interaction term.

For simplicity and in order to keep the problem analytically solvable, we consider in a first time a reduced problem for which only one triad (three waves) is retained. We note p, q and pq the wave vectors of these three waves, such that p = q + pq. Following Eq. A5, we can derive an 30 equation for each of these waves.

where the star denotes the conjugate of a complex.

Eq. A7 becomes:

where Δ = pqp-.

The process we want to study in the framework of our numerical simulations is the instability of a primary wave of wavevector p. We linearize the Eq. A8 around the basic state p , with the perturbations q and p -.

Solutions of these equations are harmonics:

The instability occurs if Im( ) > 0 and Im( -) > 0. Im( ) is the growth rate of the instability.

The secondary waves q and p -are determined as the waves that give the maximum growth rate of instability.

The linear system in Eq. A9 then becomes, using

Non trivial solutions exist if the determinant of the matrix is equal to 0.

Which gives solutions for the growth rate:

Note that the two secondary waves q and p -have the same growth rate as Im( ) = Im( -). where 0 , 0 , 0 , 0 , and 0 are the period, frequency, zonal wavelength, zonal wavenumber and amplitude of the forcing (Eq. 2). The experiments are 2D and are either linear (L) or non-linear (NL). The letter of the experiment label indicates its position on the wavenumber-frequency spectrum (Fig. 1), and F . 9. (red curves) Spectrum of the meridional velocity at the location where the primary wave p is expected (red marker in Fig. 5c) and (green curves) Spectrum of the zonal velocity at the location where the long secondary wave q is expected (green marker in Fig. 5f). (a) The time wavelet transform and (b) the space wavelet transform in the zonal direction (solid line) and in the meridional direction (dashed line). The wavelet transforms are averaged over a 2 • × 2 • horizontal domain and a 20 day time interval. Vertical lines indicate the characteristics that are predicted for the primary and secondary waves (Table 2).

LIST OF TABLES