## Introduction

The El Niño–Southern Oscillation (ENSO) is the largest source of interannual climate variability in the world (McPhaden et al. **1998**). It has been largely understood as a coupled ocean-atmospheric mode of the tropical Pacific essentially described by the simple delayed-oscillator/recharge-oscillator paradigm (e.g. Cane and Zebiak **1985**; Suarez and Schopf **1988**; Battisti and Hirst **1989**; Philander & Niño **1990**; Jin and Neelin **1993 Jin 1997a**, **1997**b; Neelin et al. **1998**; Wang and Picaut **2004**). However, while providing the basic understanding of the oscillation’s main dynamical and statistical features, this linear view of ENSO fails to account for its spatial diversity and temporal complexity, manifested for instance by large variations in amplitude, events duration, predictability, seasonal phase locking, a broad spectrum… (Timmermann et al. **2018**; Capotondi et al. **2015**; Stein et al. **2014**; Kug et al. **2009**; Ashok et al. **2007**). In particular, ENSO’s interaction with the annual cycle of the climate system has been recognised as a major source for its temporal complexity through the generation of “weak chaos” (Jin et al. **1994**; Tziperman et al. **1994**). Recently, the simple inclusion of seasonally varying parameters, such as the air-sea coupling strength, in the framework of the recharge oscillator has allowed this model to capture ENSO seasonal synchronisation characteristics and in particular the ENSO seasonal variance modulation (Stein et al. **2014**). Following this line of research, recent efforts have focussed on this phase locking and related ENSO climate modes known as ENSO-combination modes (Stuecker et al. **2013**, **2015**, **2017**). In short, if we represent the frequency of the annual harmonic as 1 cycle per year and an ENSO main period as 1/*f* years for instance, nonlinear interaction between the annual cycle and ENSO will produce ENSO responses with frequencies *n*+*mf* and *n*-*mf* where *n* and *m* are integers. Normally the response is dominated by the low order combinational near-annual frequencies 1*-f* and 1*+f* (i.e. *n =* 1, *m =* 1), which were in fact noted in earlier ENSO studies (e.g. Jin et al. **1994**, **1996**) and observations (Jiang et al. **1995**) nearly two decades ago. This represents an important ENSO dependent mode and a general route of ENSO frequency cascade, which have been suggested to contribute to potential predictability of climate variability associated with ENSO (Zhang et al. **2007**).

In addition, a number of studies indicates that higher frequency processes or stochastic atmospheric forcing, such as the Madden-Julian Oscillations (MJO) or Westerly Wind Bursts (WWB) in the Western Pacific, have significant effects on ENSO characteristics (e.g. Kleeman and Moore **1997**; Moore and Kleeman **1999**; Kessler and Kleeman **2000**; Perigaud and Cassou **2000**; Yu et al. **2003**; Lengaigne et al. **2004**; McPhaden **2004**). In turn, there are also evidences suggesting that MJO/WWB may be modulated by ENSO (Zhang and Gottschalck **2002**; Yu et al. **2003**; Eisenman et al. **2005**; Hendon et al. **2007**; Kug et al. **2008**; Tang and Yu **2008**; Kug et al. **2009**). Eisenman et al. (**2005**) showed that modulated WWB (i.e. multiplicative noise) become an effective source for destabilising the ENSO mode, as demonstrated theoretically by Jin et al. (**2007**). Numerical studies have shown that the ENSO state dependent MJO/WWB events stochastic forcing plays a critical role in the excitation, amplification and termination of El Niño events, in generating diverse ENSO behaviours and in casting difficulties for ENSO predictions.

On the other side of the Pacific basin, another intraseasonal oceanic process and its connection to ENSO and the annual cycle have been studied but have not received quite as much attention as the MJO/ENSO relationship. Tropical Instability Waves (TIWs) appear generally in the equatorial Eastern Pacific in June and persist until early in the following year as westward-propagating oscillations of the temperature front between the cold upwelled equatorial water and the warmer water to the north (Legeckis **1977**; Düing et al. **1975**). TIWs are generated by meridional shears instabilities related to the strong equatorial currents (Philander **1978**; Cox **1980**), in particular through the barotropic energy conversion associated with large meridional shears between the North Equatorial counter current and South Equatorial current (Im et al. **2012**) and between the Equatorial Undercurrent and South Equatorial Current (Luther and Johnson **1990**; Lyman et al. **2007**). They can also arise from the SST meridional gradient across the Cold Tongue through a baroclinic instability (Hansen and Paul **1984**; Wilson and Leetmaa **1988**). Their zonal wavelength, period, and phase speed are about 1000–2000 km, 20–40 days and 0.3–0.6 m/s respectively (Qiao and Weisberg **1995**). In addition, Wang and McPhaden (**1999**) used the Tropical Atmosphere Ocean buoy array to show that the seasonally varying TIWs were responsible for warming the equatorial Cold Tongue mostly through meridional heat advection. The intensified TIWs during the boreal summer/fall increase the tropical eastern Pacific SST due to the warm thermal advection by anomalous currents, with a rate estimated of up to 1 °C/month by the modelling study of Im et al. (**2012**); thus leading to an interactive feedback between seasonal and intraseasonal variations in this region. Wang and McPhaden (**2000**, **2001**) and An (2008) also highlighted the interannual modulation of TIWs by ENSO. During El Niño, TIWs are suppressed and induce an anomalous cooling, therefore acting as a negative feedback to ENSO. This feedback is stronger during La Niña than during El Niño and can partly explain ENSO asymmetry. In summary, TIWs are an integral part of the Eastern Pacific upper-ocean heat budget. They have an impact on the local SST seasonal cycle as well as on ENSO characteristics and they can be strongly modulated at seasonal and interannual timescales. Therefore, similarly to MJO/WWB, they can be seen as an oceanic-sourced state dependent stochastic forcing (transient) of the climate system.

In this study, we propose a simple linear model to study the covariance of the stochastically forced TIWs activity. Our approach differs from a dynamical theory and is rather based on simple parameterizations of the interactions between the dominant time scales. In particular, the inclusion of a time varying damping rate allows capturing the modulation of TIW variability by ENSO and the annual cycle. The rest of the paper will be structured as follows: in Section 2, we present the datasets and indices designed to represent TIWs spatio-temporal characteristics. Section 3 details our time-varying damping rate linear model, the main steps to approximate the analytical solution of the TIWs amplitude modulation and the results from different sensitivity experiments. The last section presents a summary and the main conclusions of this study.

## Definition of TIW indices

We use the NOAA high-resolution blended analysis of daily SST over the period from January 1^{st} 1997 until October 20^{th} 2016. This product combines observations from different platforms (satellites, ships, buoys) on a regular 1/4° degree global grid (Reynolds et al. **2007**). We first define a “TIW index”, noted *TIW*_{1}, as the simple equally spaced and weighted (but alternating sign) summation of SST anomalies at 8 fixed points along 0–6°N shown in Fig. 1a.

It should be pointed out that choosing different numbers and locations of base points (such 4, 6 or 10) and shifting them together in the zonal direction within the TIW active region give similar results. To capture TIWs westward propagation and thus the main TIW period, we define *TIW*_{2} in the same way as *TIW*_{1}, except the base points are all shifted by a fixed distance representing a 90° degree zonal phase shift (i.e. roughly 6.25° in longitude).

Thus, the TIWs propagation characteristics (amplitude and phase) are contained in the quantity *Z*:

This way of quantifying TIWs is original, yet extremely simple compared to prior methods. Historically, complex EOF (Lyman et al. **2007**) or demodulation (Kutsuwada and McPhaden **2002**) were commonly used to characterise TIWs and in general any propagating features for these methods allows evaluating both an amplitude and a phase information. Our approach doesn’t involve any mathematical transformation beforehand, is time independent and could, for instance, easily be implemented routinely to evaluate and/or parameterise TIWs in coarse-resolution models (i.e. which does not explicitely resolve them).

Figure 1a shows the time-longitude Hovmöller diagram of daily SST anomalies used to obtain such TIWs indices zoomed over the period 1997–2011. The corresponding extracted *TIW*_{1} index is shown on Fig. 1b and its global wavelet spectrum on Fig. 1c. The TIW phase diagram is shown in Fig. 2b and exhibits a regular wave propagation pattern. The regression of daily SST anomalies onto *TIW*_{1} shown in Fig. 2a displays the typical equatorial TIWs cusps at the right locations as compared with the spatial pattern of the daily SSTA’s 6^{th} mode of the Emprirical Orthogonal Functions decomposition, which accounts for TIW variability (Supplementary material Fig. S1a). This confirms the ability of our indices to capture both TIW propagation (i.e. phase and amplitude) and spatial structure. The TIW spectra obtained from both the direct *TIW*_{1} index and the principal component from the EOF decomposition exhibit a significant energy in a broad 20–40 days band as noted by previous studies (Figs. 1c and S1c). A preliminary visual inspection of *TIW*_{1} timeseries also reveals a strong annual and interannual modulation of its enveloppe (i.e. TIWs activity/amplitude, cf. Fig. 2b). The *TIW*_{2} index is nearly a 9-day phase shift of *TIW*_{1} and has similar spatio-temporal and spectral features (e.g. propagation, amplitude… Not shown).

## A linear model with a time-varying damping rate for TIW amplitude

### Model description and parameters

A stochastically-forced model of climate variability, in which slow changes of climate (i.e. low frequency modulation) are explained as the integral response to continuous white noise forcing, has first been introduced by Hasselmann (**1976**) and Frankignoul and Hasselmann (**1977**). It can be written as follows: $\frac{dZ}{dt}=-m(t)Z+\omega (t)$, where *Z* is the climate variable, *ω* the white noise forcing and *m*, the resistance of Z to changes produced by the forcing). Starting from the observation that most climate variables vary seasonally, de Elvira and Lemke (**1982**) extended this original stochastic climate model by including an annual feedback in *m(t)*. For instance, Stuecker et al. (**2017**) used this extended model to study the variability of the Indian Ocean Dipole. In their formulation, $m(t)=-{\lambda}_{0}+{\lambda}_{A}\hspace{0.17em}\text{cos}\hspace{0.17em}({\omega}_{a}(t)+\varphi )$, where *λ _{0}* is the mean SST damping rate,

*λ*the SST damping rate annual cycle and

_{A}*ϕ*its phase. Here, we test the hypothesis that the modulation of TIW amplitude is not only a response of the tropical Eastern Pacific annual cycle but also of ENSO and possibly of the nonlinear interaction between these two modes, therefore implying potential rectifying effects of TIWs low-frequency variability onto both ENSO and the annual cycle. The hypothesis of TIWs amplitude modulation being partly driven by a nonlinear interaction between these two modes originates from the observed strong seasonal modulation of TIWs amplitude variance (cf. Fig. 3). Similarly to ENSO synchronisation to the annual cycle, which manifests in the tendency of ENSO events to peak during boreal winter (Stein et al.

**2014**), both TIW amplitude (Fig. 3a, black line) and TIW amplitude’s variance (Fig. 3b, black line) exhibit a maximum during the boreal summer/fall. Therefore, our postulate is that TIWs amplitude is the result of a stochastically-forced wave with a seasonally and interannually varying damping rate. As an extension of de Elvira and Lemke (

**1982**) and Stuecker et al. (

**2017**), we add an interannually varying damping rate and formulate our model as follows:

*dZ*/

*dt*is the TIW amplitude tendency,

*ω*(

*t*) is a white noise forcing and

*Nino*3(

*t*) the interannual (ENSO) forcing. The strength of the noise impacts the amplitude of the solution but has no effect on its spectral properties. The primary objective of the present study is to assess TIWs amplitude’s variability, thus in the following we will only present standardised daily time series. Two daily ENSO forcings of different complexity are tested: (

*i*) the classic Niño3 index (SST anomalies averaged in the Eastern Pacific box 90°W-150°W; 5°S-5°N), (

*ii*) an ideal ENSO forcing, i.e. a simple sinusoidal interannual oscillation, $\mathit{Nino}3(t)=\hspace{0.17em}\text{cos}\hspace{0.17em}\frac{2\pi (t-{\varphi}_{2})}{{T}_{E}}$ with

*T*

_{E}= 1100 days, i.e. the ENSO main period (cf. Niño3 spectrum, i.e. red line on Supplementary material Fig. S2a) and

*ϕ*

_{2}the phase for the interannual damping rate.

*T*= 36 days and

*T*= 365 days are respectively the TIW and annual cycle (AC) periods. Although TIWs cover a wide range of periods (roughly between 15 and 40 days, cf. Fig. 1c, S1c and Lyman et al.

_{A}**2007**), we chose to restrict the characterisation of this phenomenon to its dominant period (i.e. 36 days) as identified on the

*TIW*spectrum (cf. Fig. 1c).

_{1}*γ*

_{0},

*γ*

*and*

_{A}*γ*are respectively the TIWs mean, annual and interannual damping rates.

_{N}*ϕ*is the phase for the annual damping rate and is chosen so that TIWs amplitude reaches a maximum in August/September and a minimum in March/April (

*ϕ=*120 days). Based on the TIWs indices auto and cross correlations (Fig. 2c), we estimate a mean TIW damping rate

*γ*

_{0}of 1/23 days

^{−1}.

Several options are available to estimate the model’s parameters. First of all, since our model’s fornulation is similar to a Linear Inverse Model (for a general description, see Penland and Matrosova (**1994**); Penland and Sardeshmukh (**1995**) and references therein), we could simply perform a multiple linear regression of TIW amplitude onto the Eastern Pacific SST seasonal cycle and the Niño3 index. This would lead to the estimation of the model’s optimal parameters. However, one can argue that this “tuning method” is too specific to the data and time period considered and leads to an artificially too optimised solution. In order to stay as general as possible, and in a similar fashion as for *γ*_{0}, we choose to estimate *γ _{A}* and

*γ*only from the decorrelation time of observed TIWs. To do so we contrast auto correlations of (

_{N}*TIW*

_{1}-

*TIW*

_{2}) and cross correlations between

*TIW*

_{1}and

*TIW*

_{2}to evaluate the TIWs propagation e-folding times (or growth rates) during different seasons of climatological years (for

*γ*) and different phases of ENSO (for

_{A}*γ*) respectively. We distinguish between El Niño, La Niña and climatological (or neutral) years following the NOAA Climate Prediction Center classification. ENSO seasons are from November to February and 1997/98, 2002/03, 2006/07, 2009/10, 2015/16 are considered El Niño events; while 1998/99, 1999/2000, 2000/01, 2007/08, 2010/11, 2011/12 are considered La Niña events. 1997/98, 2009/10 and 2015/16 are considered intense El Niño events. 1998/99, 1999/2000 and 2007/08 are considered intense La Niña events. 2001/02, 2004/05, 2005/06, 2008/09, 2012/13, 2013/14 and 2014/15 are classified as climatologial calendar years (i.e. neutreal ENSO conditions). The estimations of these different TIWs growth rates are presented in Fig. 4a and c and show, as expected, a longer persistance of TIWs activity in summer/fall and during La Niña than in winter/spring (Fig. 4a) and during El Niño (Fig. 4c). Following the notations introduced in Fig. 4b and d, we estimate

_{N}*γ*and

_{A}*γ*as follows:

_{N}${\gamma}_{A}=\frac{\delta {\gamma}_{A}}{2}\frac{\delta {T}_{A}}{\delta {T}_{\mathit{TOT}}}\mathrm{and}{\gamma}_{N}=\frac{\delta {\gamma}_{N}}{2}\frac{\delta {T}_{N}}{\delta {T}_{\mathit{TOT}}};$*δT _{TOT}* is the maximum amplitude of the total SST signal (climatology plus interannual anomalies), i.e.

*δT*= max(

_{TOT}*δT*) over the total time period. Based on this evaluation, we estimate the following rounded parameters:

_{A}+δT_{N}*γ*

_{A}/γ_{0}

*=*0.30 and

*γ*

_{N}/γ_{0}

*= −*0.50. Interestingly, while the TIWs seasonal growth rate changes rather linearly with the climatological equatorial SST (and SST meridional gradient, cf. Fig. 4b), the TIWs interannual growth rate exhibits a strong nonlinear relationship with the SST interannual anomalies (Fig. 4d). The nonlinear relationship between the interannual TIW growth rate and equatorial SST anomalies is best represented by a 3

^{rd}order polynom, which accounts for the asymmetrical nature of this nonlinearity (i.e. the ENSO asymmetrical feedback on TIW amplitude; An (

**2008**)). To include this nonlinear effect, a more refined version of the model is adopted hereafter; and instead of Equation (5), we use the following:

Using a 3^{rd}-order polynomial fit, in the least square sense, we estimate ${\gamma}_{{N}^{3}}=0.007$. It is worth noting that, while this effect might slightly enhance the asymmetry of TIWs activity between the ENSO phases, it will not be its main contributor. The TIWs activity’s asymmetry will mostly emerge from the interannually varying TIWs growth rate’s term (*γ _{N}Nino*3(t) in Eq. 6).

### Approximation of the TIW amplitude’s analytical solution and parameters sensitivity

As long as |*γ _{A}*/

*γ*| < 1, |

_{0}*γ*/

_{N}*γ*| < 1 and |

_{0}*γ*/

_{N}^{3}*γ*| < 1, the oscillating solutions

_{0}*Z*stay in a stable damped regime and we can use a second order Taylor series expansion (O(ε

^{2})) to approximate the analytical solution of TIWs amplitude |

*Z*|. After some mathematical steps detailed in the appendix A, it can be shown that the low-frequency modulation of TIWs amplitude has the following form:

*|Z|*

^{2}(i.e. the low-frequency modulation, diagnosed in the following as the 50-days running mean of the solution

*Z*

^{2}) and $m(t)=({\gamma}_{A}\hspace{0.17em}\text{cos}\hspace{0.17em}\frac{2\pi (t-\varphi )}{{T}_{A}})+({\gamma}_{N}N\mathit{ino}3(t))+({\gamma}_{{N}^{3}}N\mathit{ino}3{(t)}^{3})$

The 2^{nd} order approximation of our simple model’s solution exhibits similar spectral properties as the observed TIW amplitude (cf. Fig. 3c black and magenta lines). In particular, this model is able to produce the variability associated with the combination mode (C-MODE) between the annual cycle and ENSO. This variability originates from the nonlinear interaction between ENSO and the annual cycle found within the *m ^{2}*(

*t*) term of the analytical solution:

As explained in Section 1, this combination tone mostly manifests itself at a couple frequencies *f _{C-MODE}* =

*f*±

_{AC}*f*(Fig. 3c, frequency bands highlighted in purple) and cannot be explained by any direct external forcing (Stuecker et al.

_{ENSO}**2013**;

**2017**).

Additionally, we performed a parameter sensivity analysis presented in Supplementary material Fig. S3 (cf. supplemental material). It shows the sensitivity of some features of TIW amplitude’s analytical solution to the |*γ _{A}/γ*

_{0}| and |

*γ*

_{N}/γ_{0}| ratio values. This not only corroborates the parameters estimation presented in the previous section based on a physical analysis of the observed TIWs features (and not from an optimal tuning method through linear regressions), but also illustrates the robustness of our model and its analytical solution within the parameters domain.

### Sensitivity of numerical and analytical solutions to forcing

Now that the TIW model can reproduce realistically the main spectral and temporal features of the observed TIW amplitude variability, we explore the sensitivity of the analytical and numerical solutions to different forcing terms. To do so, we compare the spectral and temporal properties of the TIW amplitude variability generated by a 50-members ensemble of three different versions of the linear model: the full model forced by both the annually and interannually TIW varying damping rates (noted TIWxTOT, i.e. all terms included, *γ _{N} ≠* 0

*, γ*0 and

_{A}≠*γ*0), the model forced only by the annually varying TIW damping rate (noted TIWxAC, i.e.

_{N}^{3}≠*γ*0 and

_{N}=*γ*0) and the model forced only by the interannually varying TIW damping rate (noted TIWxNino, i.e.

_{N}^{3}=*γ*0). All versions of the model are integrated from January 1

_{A}=^{st}1997 to August 31

^{st}2016 using a 4

^{th}order Runge-Kutta method and the parameters values from Section 3.1.

Figure 5a shows the timeseries of the observed (black line) and simulated TIW amplitude (blue line for the ensemble average and thin Coloured lines for all individual members integrated using the TIWxTOT version) along with the 2^{nd} order approximation of the analytical solution (red line). The full model (TIWxTOT i.e. *γ _{N} ≠* 0

*, γ*0 and

_{A}≠*γ*0) simulates realistically the TIW amplitude variability, as demonstrated by a good correlation between the observed and the analytical solution, respectively ensemble averaged, TIW amplitude (

_{N}^{3}≠*r*= 0.72, respectively 0.75; cf Fig. 5a). Figure 5b (respectively Fig. 6a) summarises the spectral properties of the analytical solution (respectively ensemble-averaged solutions) of the TIW amplitude simulated by different forcing experiments of the linear model. Additionally, Fig. 3c presents the analytical solutions power spectra from different versions of the model (i.e. TIWxTOT, TIWxAC and TIWxNino) and for a wide range of parameters values (i.e. $0.2\le \left|\frac{{\gamma}_{A}}{{\gamma}_{0}}\right|,\left|\frac{{\gamma}_{N}}{{\gamma}_{0}}\right|\le 0.8$). The TIW amplitude generated by the fully forced version of the model exhibits a relatively similar spectrum (magenta line) as the one from the observed TIW amplitude (thick black line), significant in particular for the key annual, ENSO and C-MODE frequencies. Interestingly, we note the presence of a significant peak of near decadal variability around 7–10 years

^{−1}in the spectrum of both the observed and simulated TIW amplitude (as well as in the spectrum of Niño3 index, cf. red line in Supplementary material Fig. S1a). However, although barely significant as compared to a random red noise process, this frequency peak is stronger in the spectrum of the TIW amplitude generated by the full forcing experiment (TIWxTOT) and less energetic in the ENSO-only case (see comparison between green and magenta spectra in Figs. 3c, 5b and 6a). This suggests that this near-decadal TIW amplitude variability might be amplified by rectification processes emerging from the nonlinear interaction between the annual cycle, ENSO and the seasonal and interannual modulations of TIW amplitude.

The asymmetric feedback between the ENSO state and TIW amplitude first evidenced by An (**2008**) is also captured by this simple model with an enhanced, respectively suppressed TIW amplitude during La Niña, respectively El Niño, as shown in Fig. 6c. Interestingly, this simple model also seems to capture the nonlinear interaction between the Eastern Pacific annual cycle and TIWs amplitude’s annual modulation, as shown in Fig. 6b by a shift in phase of ≈40 days and a broadening from the analytical solution to the forcing’s seasonal cycle. This indicates a potential nonlinear rectification of the tropical Eastern Pacific upper-ocean seasonal heat budget by the modulation of TIWs amplitude.

As evidenced on Supplementary material Fig. S2a (red line), the direct ENSO forcing displays a significant peak of variability in the C-MODE frequency band. Therefore, there is a chance that the C-MODE peak in the simulated TIW amplitude emerges from this prescribed ENSO forcing rather than from the nonlinear combination between the ENSO and annual frequencies. The use of interannual forcing of different complexity and in particular the comparison between TIW amplitude solutions inferred from an ideal (i.e. purely sinusoidal by construction and thus containing no C-MODE variability) and the prescribed ENSO forcing confirm that most of the C-MODE variability emerges from the nonlinear interaction between the Annual Cycle and ENSO (cf. Figs. 3c and 6a as compared to Supplementary material Fig. S4c which shows the spectra of the different versions of the model integrated using an ideal ENSO forcing).

Additionally, Supplementary material Fig. S2b shows the spectrum of the observed and simulated TIW, i.e. the *TIW*_{1} index (real part of Z in the nonlinear model). The richness in TIW spectrum (20–40 days^{−1}) only emerges when both the annually and inter-annually varying TIW damping rates are included in the model and cannot be accounted for without the nonlinear interaction between the Eastern Pacific Annual Cycle and ENSO. Although beyond the scope of this paper, the complex variability of the TIW high-frequency process arises partly from nonlinear dynamics in the Eastern Pacific. TIW amplitude modulation can potentially influence the local SST variability at different timescales, which in turn can alter the typical patterns of the annual and interannual atmospheric variability in this key region where sits the Inter-Tropical Convergence Zone.

## Conclusion and discussion

In this paper, we developed a simple TIWs stochastic linear model forced by ENSO and the annual cycle to assess the deterministic low frequency modulation of TIWs variance. In particular, we study the effect of the Cold Tongue SST annual cycle and ENSO on the modulation of TIW activity. Both the analytical and numerical solutions of our simple model depict realistically the basic features of TIWs amplitude variability; namely the seasonal increase in TIWs variance during summers and falls and the nonlinear relationship with the ENSO phase characterised by a suppression, respectively increase of TIW activity during El Niño, respectively La Niña. The novelty result evidenced through this innovative framework is that a significant fraction (diagnosed by a significant spectral peak in the C-MODE frequency band as compared to a red noise process) of TIWs amplitude modulation emerges from the deterministic nonlinear interaction between ENSO and the annual cycle and is therefore predictable. This reflects a similar combination mode dynamic as the one described by Stuecker et al. (**2013**; **2017**) but with a different seasonality (enhanced interaction during the northern hemisphere summer/fall). Therefore, this study revealed a new ENSO combination mode and a new route of ENSO frequency cascade through both the oceanic and atmospheric responses to TIW-induced changes in the upper-ocean heat transport. The deterministic variability of this state-dependent oceanic stochastic forcing will potentially trigger different oceanic and atmospheric teleconnections in the tropical climate system and deserves more in-depth investigation.

Additionally, the observed interannual anomalies of the 50-days running TIW amplitude are significantly correlated with those from the model’s 50-member ensemble average (not shown) and the analytical solution (Fig. 7a). Moreover, this variability is negatively correlated with ENSO. Since, in our simple TIW model, ENSO modulates TIW activity, it can be easily proved that the interannual anomalies of TIW variance will induce a proportional and significant variation in TIW-induced nonlinear dynamic heating (NDH) owing to the strong TIW heat transport meridional convergence -∇·(*v*′*T*′) (mainly due to the strong SST meridional gradient). This will be the object of an upcoming article but preliminary results reveal a high correlation between the TIWs-induced nonlinear meridional heat flux (-*v’T’*) and the interannual anomalies of both the observed (*r* = 0.57) and theoretical (*r* = 0.78) TIW amplitude (not shown). Thus, our simple model provides a simple analytical theory to describe the asymmetric feedback between the ENSO state and TIW activity, which was first empirically introduced by An (**2008**) into a simple ENSO model to assess the TIW contribution to the observed El Niño/La Niña asymmetry. We can now give an analytical formulation for the ENSO-induced TIW NDH feedback onto ENSO as follows:

This formulation not only accounts for the response to the ENSO-Annual Cycle interaction but also for potentially significant higher order nonlinear interactions, i.e. following the notation introduced in Section 1, at frequencies 2*f*, 4*f*, 6*f* and 1 ± 3*f*. Despite the fact that our model is linear, this approximate analytical solution, obtained by a perturbation method, demonstrates clearly that the TIW feedback onto ENSO is not only nonlinear but also strongly modulated by the seasonal cycle (cf. Fig. 7b and c). This nonlinearity comes from that fact that this feedback operates through the TIW nonlinear heat transport convergence -∇·(*v*′*T*′), i.e. NDH processes. We plan to conduct further investigations motivated by this conceptual advancement, to assess the interactions between ENSO and TIW activity.

Furthermore, the simple theoretical framework presented here can be generalised to assess the modulation by the annual cycle, ENSO and their combined effect of any so-called state-dependent transients of the tropical climate system. A stochastically forced linear model of high frequency processes with periods much shorter than ENSO, with an annually and interannually varying damping rate, provides a good approximate to infer nonlinear deterministic rectifications through the covariance fields. A follow-up application of this new framework is the study of the deterministic modulation of MJO/WWB activity due to ENSO, the annual cycle and their nonlinear combination. This will be the topic of an oncoming research article.