1.

## 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, 1997b; 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.

2.

## Definition of TIW indices

We use the NOAA high-resolution blended analysis of daily SST over the period from January 1st 1997 until October 20th 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 TIW1, 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.

((1))
$TI{W}_{1}\left(t\right)=\sum ±\mathit{SST}\prime \left(t,\mathit{nodes}\right)$
Fig. 1.

(a) Hovmöller diagram of daily SST anomalies between 180° and 100°W (0°-6°N latitudinal average). The red and blue horizontal lines represent the longitudinal positions of the nodes used to generate TIW indices. Red (blue) nodes/lines are positive (negative). (b) $TI{W}_{1}\left(t\right)=\sum ±\mathit{SST}\prime \left(t,\mathit{nodes}\right)$ time series. (c) Global wavelet spectrum of TIW1. See the mathematical definition of TIW1 and TIW2 indices in Section 2.

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 TIW2 in the same way as TIW1, 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).

((2))
$TI{W}_{2}\left(t\right)=\sum ±\mathit{SST}\prime \left(t,\mathit{nodes}+{6.25}^{\circ }\right)$

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

((3))
$Z=TI{W}_{1}+\mathit{iTI}{W}_{2}$
and the TIW amplitude is simply:
((4))
$|Z|=TI{W}_{12}=\sqrt{TI{W}_{1}^{2}+TI{W}_{2}^{2}}$

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 TIW1 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 TIW1 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 6th 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 TIW1 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 TIW1 timeseries also reveals a strong annual and interannual modulation of its enveloppe (i.e. TIWs activity/amplitude, cf. Fig. 2b). The TIW2 index is nearly a 9-day phase shift of TIW1 and has similar spatio-temporal and spectral features (e.g. propagation, amplitude… Not shown).

3.

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

3.1.

### 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\left(t\right)Z+\omega \left(t\right)$, 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\left(t\right)=-{\lambda }_{0}+{\lambda }_{A} \text{cos} \left({\omega }_{a}\left(t\right)+\varphi \right)$, where λ0 is the mean SST damping rate, λA the SST damping rate annual cycle and ϕ 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:

((5))
$\begin{array}{c}Z=TI{W}_{1}+\mathit{iTI}{W}_{2}\\ \frac{dZ}{dt}=\left[-\left({\gamma }_{0}+\frac{2i\pi }{T}\right)+\left({\gamma }_{A} \text{cos} \frac{2\pi \left(t-\varphi \right)}{{T}_{A}}\right)+\left({\gamma }_{N}N\mathit{ino}3\left(t\right)\right)\right]Z+\omega \left(t\right)\\ ↓ ↓ ↓ ↓\\ \mathrm{ }\mathit{TIW}\mathit{oscillator} \mathit{TIW}×AC \mathit{TIW}×\mathit{ENSO} \mathit{Noise}\end{array}$
where dZ/dt is the TIW amplitude tendency, ω(t) is a white noise forcing and Nino3(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\left(t\right)= \text{cos} \frac{2\pi \left(t-{\varphi }_{2}\right)}{{T}_{E}}$ with TE = 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 TA = 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. 2007), we chose to restrict the characterisation of this phenomenon to its dominant period (i.e. 36 days) as identified on the TIW1 spectrum (cf. Fig. 1c). γ0, γA and γN are respectively the TIWs mean, annual and interannual damping rates. ϕ 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.

Fig. 2.

(a) Regression of daily SST anomalies onto TIW1 index (1997–2015). Dots denote the 95% statistical significance level. Black filled triangles and circles show the positions of longitudinal nodes (positive for circles, i.e. red lines on Fig. 1a, and negative for triangles, i.e. blue lines on Fig. 1a). The bottom right panel (b) is the (TIW1, TIW2) phase diagram. (c) Auto correlation of TIW1-TIW2 in black and cross correlation between TIW1 and TIW2 in red.

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 γN only from the decorrelation time of observed TIWs. To do so we contrast auto correlations of (TIW1-TIW2) and cross correlations between TIW1 and TIW2 to evaluate the TIWs propagation e-folding times (or growth rates) during different seasons of climatological years (for γA) and different phases of ENSO (for γN) 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 γA and γN as follows:

Fig. 3.

(a) Annual cycle (monthly mean climatology) of the observed TIW amplitude (1997-2015) (black) and of the ensemble average of the 2nd order TIW amplitude analytical solution inferred with different reasonable values of γA0 and γN0 ($0.2\le |\frac{{\gamma }_{A}}{{\gamma }_{0}}|,|\frac{{\gamma }_{N}}{{\gamma }_{0}}|\le 0.8$ ) and different versions of the model, i.e. using only a annually-varying TIW damping rate (TIWxAC) in red, using only an interannually-varying TIW damping rate (TIWxNino) in green and using both (TIWxTOT) in magenta. (b) Same but for the monthly variance of TIW amplitude. (c) FFT single-sided spectrum of the observed TIW amplitude (i.e. |Z| = TIW12, black) and ensemble averaged of FFT single-sided amplitude spectra of the 2nd order approximation of the TIW amplitude analytical solution for the different versions of the model (same colour code as (a)). The dashed blue line indicates the averaged spectrum of 5000 generations of a random red noise process. The Coloured vertical bands indicate select time scales: ENSO (2.5-5 years) in green, the Annual Cycle (AC, 340-390 days) in pink, the Semi Annual Cycle (SAC, 160-200 days) in orange, and the combination mode (C-MODE) in purple. Period: Jan 1st 1997-Dec 31st 2015.

${\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}}};$δTTOT is the maximum amplitude of the total SST signal (climatology plus interannual anomalies), i.e. δTTOT = max(δTA+δTN) over the total time period. Based on this evaluation, we estimate the following rounded parameters: γA0= 0.30 and γN0= −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 3rd 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:

((6))
$\frac{dZ}{dt}=\left[-\left({\gamma }_{0}+\frac{2i\pi }{T}\right)+\left({\gamma }_{A} \text{cos} \frac{2\pi \left(t-\varphi \right)}{{T}_{A}}\right)+\left({\gamma }_{N}N\mathit{ino}3\left(t\right)\right)+\left({\gamma }_{{N}^{3}}N\mathit{ino}3{\left(t\right)}^{3}\right)\right]Z+\omega \left(t\right)$

Using a 3rd-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 (γNNino3(t) in Eq. 6).

3.2.

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

As long as |γA/γ0| < 1, |γN/γ0| < 1 and |γN3/γ0| < 1, the oscillating solutions 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:

((7))
$〈|{Z}^{2}|〉=K{e}^{-{\gamma }_{0}}\left[1+\frac{m\left(t\right)}{{\gamma }_{0}}+\frac{{m}^{2}\left(t\right)}{{\gamma }_{0}^{2}}\right]$
where $〈|{Z}^{2}|〉$ represents the slow variation of TIWs amplitude |Z|2 (i.e. the low-frequency modulation, diagnosed in the following as the 50-days running mean of the solution Z2) and $m\left(t\right)=\left({\gamma }_{A} \text{cos} \frac{2\pi \left(t-\varphi \right)}{{T}_{A}}\right)+\left({\gamma }_{N}N\mathit{ino}3\left(t\right)\right)+\left({\gamma }_{{N}^{3}}N\mathit{ino}3{\left(t\right)}^{3}\right)$

The 2nd 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 m2(t) term of the analytical solution:

((8))
$C-\mathit{MODE}=\frac{2}{{\gamma }_{0}}×\left({\gamma }_{A} \text{cos} \frac{2\pi \left(t-\varphi \right)}{{T}_{A}}\right)×\left({\gamma }_{N}N\mathit{ino}3\left(t\right)\right)$

As explained in Section 1, this combination tone mostly manifests itself at a couple frequencies fC-MODE = fAC ± fENSO (Fig. 3c, frequency bands highlighted in purple) and cannot be explained by any direct external forcing (Stuecker et al. 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 |γA0| and |γN0| 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.

3.3.

### 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, γA 0 and γN3 0), the model forced only by the annually varying TIW damping rate (noted TIWxAC, i.e. γN = 0 and γN3 = 0) and the model forced only by the interannually varying TIW damping rate (noted TIWxNino, i.e. γA = 0). All versions of the model are integrated from January 1st 1997 to August 31st 2016 using a 4th order Runge-Kutta method and the parameters values from Section 3.1.

Fig. 4.

(a) Average of (TIW1-TIW2) autocorrelations (in black) and (TIW1,TIW2) cross correlations (in red) logarithm for climatological years (i.e. neither El Niño nor La Niña years) between November and February. The black dotted line represents the best linear fit to these curves local maxima or in other words the estimation of TIWs propagation e-folding time during neutral years. Similarly, the green (magenta) dotted line represents the best linear estimation of TIWs propagation e-folding time for El Niño (La Niña) years. The red (blue) dotted line represents the best linear estimation of TIWs propagation e-folding time for extreme El Niño (La Niña) years. (b) The inverse of these e-folding times (in other words, TIWs growth rates) are represented as a function of the daily equatorial (2°S-4°N; 180°W-100°W) SST anomalies averaged during the corresponding periods and years (in red). The black curve is the same but for TIWs growth rates as a function of the meridional gradient of SST anomalies, i.e. the difference between the equatorial region average and the [7°N-13°N; 180°W-100°W] average. (c) Same as (a) but during different seasons of neutral years: Auto (black) and cross (red) correlations logarithm averaged during February-September. The corresponding e-folding time linear estimation is shown in the dotted black line. The green, respectively magenta, dotted line represents the best linear estimation of TIWs propagation e-folding time during the summer and fall (July-September), respectively spring and winter (February-April) of neutral years. (d) Same as (b) but the TIWs seasonal growth rates are shown as a function of the corresponding equatorial SST average in red and the meridional climatological SST gradient in black.

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 2nd order approximation of the analytical solution (red line). The full model (TIWxTOT i.e. γN 0, γA 0 and γN3 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 (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 |\frac{{\gamma }_{A}}{{\gamma }_{0}}|,|\frac{{\gamma }_{N}}{{\gamma }_{0}}|\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 TIW1 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.

4.

## 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 -∇·(vT′) (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:

((9))
$ND{H}_{\mathit{TIW}}\approx -\nabla .\left(v\prime T\prime \right)\approx \left(2\frac{{\gamma }_{A}{\gamma }_{N}}{{\gamma }_{0}^{2}}\text{cos} \left(\frac{2\pi \left(t-\varphi \right)}{{T}_{A}}\right)×\mathit{Nino}3\left(t\right)\right)+{\left(\frac{{\gamma }_{N}}{{\gamma }_{0}}\right)}^{2}N\mathit{ino}3{\left(t\right)}^{2}+{\left(\frac{{\gamma }_{{N}^{3}}}{{\gamma }_{0}}\right)}^{2}N\mathit{ino}3{\left(t\right)}^{6}+\left(2\frac{{\gamma }_{A}{\gamma }_{{N}^{3}}}{{\gamma }_{0}^{2}}\text{cos} \left(\frac{2\pi \left(t-\varphi \right)}{{T}_{A}}\right)×\mathit{Nino}3{\left(t\right)}^{3}\right)+2\frac{{\gamma }_{A}{\gamma }_{{N}^{3}}}{{\gamma }_{0}^{2}}N\mathit{ino}3{\left(t\right)}^{4}$

Fig. 5.

(a) Individual members (thin Coloured lines) and ensemble average (thick blue line) of the daily TIW amplitude numerical solutions using the full TIWxTOT version of the model. The black line shows the observed daily TIW amplitude. The thick red line is the 2nd order analytical solution of TIW amplitude. All time series are normalised and smoothed using a 50-days running mean. (b) Spectra of the observed and 2nd order approximation of the analytical solution of TIW amplitude from different sensitivity experiments. Thick black lines for observation, red for the experiments using only a annually-varying TIW damping rate (TIWxAC), green for the experiments using only an interannually-varying TIW damping rate (TIWxNino) and magenta for the experiments using both (TIWxTOT). The dashed blue line indicates the averaged spectrum of 5000 generations of a random red noise process.

Fig. 6.

(a) Spectra of the observed and ensemble average of TIW amplitude (Z = TIW12) simulated by different sensitivity experiments (the Coloured bands representing different typical timescales, see Fig. 5b caption). Thick black lines for observation, red for the experiments using only a annually-varying TIW damping rate (TIWxAC), green for the experiments using only an interannually-varying TIW damping rate (TIWxNino) and magenta for the experiments using both (TIWxTOT). The dashed blue line indicates the averaged spectrum of 5000 generations of a random red noise process. (b) Annual cycle of the observed (black line), the ensemble average (blue line) and the analytical solution (red) of TIW amplitude; the annual cycle of the model forcing is shown in green. (c) Scatter plot of simulated TIW amplitude (TIWxTOT ensemble average) versus the prescribed Niño3 forcing.

Fig. 7.

(a) Observed interannual part (i.e. cleared from the 1997–2015 monthly mean climatology) of TIW amplitude (black) and interannual part of the analytical solution 2nd order approximation of TIW amplitude (red) using the following set of parameters: γ0 = 1/23 days−1; γA0 = 0.3; γN0 = −0.5 and γN3 = 0.007; The correlation between these two timeseries is 0.65. The blue line shows the Niño3 index. (b) Scatter plot of the nonlinear relationship between the observed interannual TIW amplitude and ENSO (Niño3 index); Red (summer/fall only) and blue (winter/spring only) dots highlight the seasonality of this relationship. (c) Same as (b) for the the nonlinear relationship between the 2nd order approximation of the analytical solution of the interannual TIW amplitude and ENSO.

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 2f, 4f, 6f and 1 ± 3f. 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 -∇·(vT′), 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.