To characterize the oceanic circulation associated with TIWs and IKWs, we use the daily sea level (SLA) and geostrophic currents (U,V) anomalies from the multi-satellite Data Unification and Altimeter Combination System datasets available through the Copernicus Marine Environment Monitoring Service (CMEMS) and the Copernicus Climate Change Service (C3S) (Taburet et al., 2019) at a spatial resolution of 0.25°× 0.25° over the region (10°S–10°N, 130°E–80°W) and the period 1993–2018.

To document the IKWs atmospheric forcing, we calculated the zonal wind stress (τ_x) considering the drag coefficient formulation of Trenberth et al. (1989) from the 10-m zonal winds from the Fifth Generation ECMWF Atmospheric Reanalysis of the Global Climate (ERA5) (Hersbach et al., 2018) at a 0.25°× 0.25° horizontal resolution over the same region and period than the sea level. Finally, to estimate the meridional oceanic advection associated with the TIWs-IKWs interaction mode, we use the daily sea surface temperatures (SST) from the NOAA 1/4° daily Optimum Interpolation Sea Surface Temperature version 2.1 product (Huang et al., 2020) at a 0.25°× 0.25° resolution over the same period and domain. All data were first linearly detrended. Then, we define the intraseasonal anomalies for all variables as the departure of the instantaneous daily data from the monthly mean following the method by Lin et al. (2000), which is somewhat equivalent to applying a high-pass Lanczos filter with a frequency cut-off at (1/90) days-1 (cf. Dewitte et al., 2011).

Complex Empirical Orthogonal Functions (CEOF) are an efficient tool to disentangle the spatiotemporal characteristics of the intraseasonal variability related to TIWs and KWs as this analysis provides both amplitude and phase information, and is therefore well-suited to capture propagating features (Barnett, 1983; Stein et al., 2011; Boucharel et al., 2013, 2016).

To analyze the TIWs characteristics, we use the definition of the TIWs amplitude and index from Boucharel and Jin (2020). We take a simple equally spaced and weighted nodes (positive and negative SLA anomalies) index TIW₁. As stated in Boucharel and Jin (2020), the number of nodes and their longitudinal location within the (0°–6°N) band can be chosen arbitrarily. In this study, we chose 6 (SLA) nodes within the (3°N–7°N and 165°–107°W) band that essentially fit the TIWs cusps from the related CEOF mode’s spatial pattern (see section “Satellite Data”).

Where SLA refers to the sea level intraseasonal anomalies. Then, to grasp the TIWs propagation, we define TIW₂ similarly to TIW₁ except the nodes are now shifted by 90 degrees representing half the TIWs wavelength.

In this paper, we apply this methodology to derive the TIWs activity from the sea level.

To infer the IKWs contribution to sea level anomalies and interpret the results from the CEOF decomposition, we project the SLA and zonal wind stress anomalies onto the Kelvin/Rossby meridional wave structures in the long wavelength approximation (see Cane and Sarachik, 1976). We assume like in former studies based on satellite sea level (Boulanger and Menkes, 1995; Perigaud and Dewitte, 1996) that the oceanic vertical structure along the equator can be accounted for by just the first baroclinic mode, so that, total sea level (h) and zonal current (u) anomalies can be expressed as a linear combination of meridional mode functions accounting for Kelvin and Rossby waves. This writes as follows:

Where R0u(y),R0h(y) and Rnu(y),Rnh(y) are the normalized Kelvin and the n-th normalized Rossby waves meridional structures, n denotes the meridional modes being 0 for Kelvin waves and ≥ 1 for the Rossby modes. The r_n are the corresponding wave coefficients calculated at each grid point along the equator. The projection of sea level and zonal current anomalies onto the meridional modes is performed between the northern Y_N and southern limits Y_S (10°N–10°S in this study). This latitudinal range allows capturing the poleward exponential decrease of Kelvin and Rossby waves and ensuring an approximate orthogonality property between meridional modes.

Where the ψ_n(y) are the Hermite functions, calculated using a value for phase speed that varies zonally. Its value is set to 3 m/s at 170°W, which corresponds to an average value used in the literature for the first baroclinic mode (Picaut and Sombardier, 1993; Dewitte et al., 1999). We then approximate C(x) as C(x)=C(x=170W°)Zeq¯Z170W°¯ where Zeq is the mean depth of the 20°C isotherm at 170°W and between 2°S and 2°N. Zeq was derived from the global ocean data assimilation system (GODAS, Behringer and Xue, 2004) using the 20°C isotherm depth as a proxy and for the period 1980–2018. The total SLA as well as the SLA reconstructed from the individual CEOF modes are, unless stated otherwise, projected onto the meridional modes functions using the method by Boulanger and Menkes (1995) (see their Appendix B), which does not require the knowledge of zonal current anomalies and is comparable to assuming that Rossby wave coefficients are negligible from a certain rank (see also Perigaud and Dewitte, 1996).

The zonal wind stress is also projected onto the theoretical meridional mode structures of the first baroclinic mode, and its contribution to the Kelvin wave is inferred from the following formula: τ_x (x, y, t) = Fo (x, t) R_n (y), and

Where h, u and τ_x come from observations described in section” Satellite Data.”

Results from the CEOF decomposition of the sea level intraseasonal anomalies are presented in Figure 1. The first two modes grasp 21% of the explained variance of the total sea level anomalies over the spatial domain shown in Figure 1, almost equally distributed. Since higher-order modes explain individually less than 5.4% of the variance, in the following we focused on the first two modes only. The first mode (CEOF1), characterized by a strong amplitude between 3–7°N and 120–165°W (Figure 1A), mostly captures the variability associated with TIWs (with an explained variance of 11%) although a relative maximum amplitude is also observed along the equator reminiscent of IKWs contribution to sea level anomalies. The real part of the spatial CEOF pattern (Supplementary Figure 2A) exhibits TIWs cusps similarly positioned as in Xue et al. (2020). The decreasing phase from east to west in this region shown on the CEOF spatial phase pattern (Supplementary Figure 3A) characterizes the TIWs westward propagation. The second mode (CEOF2), which explains a similar amount of variance, exhibits a comparable phase (Supplementary Figure 3B) and amplitude pattern although the amplitude along the equator is about twice as large as that of the first mode while the amplitude of the TIWs cups is smaller. Note also the exact same location of the TIWs cups but with opposite value of amplitude of the real part between CEOF1 and CEOF2 (Supplementary Figure 2), which results from the orthogonality of the two modes’s patterns and suggests that they capture the same sequence of TIWs propagations. Their summed-up contribution explains 42% of the variance of the total sea level anomalies over the region (2°S–2°N, 170–120°W). The summed-up contribution from CEOF1 and CEOF2 over the TIWs active region (3–7°N and 165–107°W) was also found to be highly correlated with both TIWs index and TIWs amplitude from Boucharel and Jin (2020) (r = 0.5 and 0.7 respectively, not shown) which all together indicates that CEOF1 + 2 capture most of TIWs variability. In the equatorial band, the zonally elongated pattern symmetrical around the equator, more pronounced for CEOF2 than for CEOF1, is indicative of the IKWs activity and suggests that the two modes capture the non-linear interaction between IKWs and TIWs. The phase patterns of the two modes along the equator (Supplementary Figure 3) confirm the eastward propagation at a phase speed estimated to 2.6 m/s consistent with former studies based on altimetry data (Cravatte et al., 2003; Dewitte et al., 2003). Their projections onto the Kelvin wave structure indicate that CEOF2 explains 24% of the variance of the Kelvin wave contribution to the full sea level anomalies average over (2°S–2°N, 170–120°W) while CEOF1 only explains 12.9% (see Table 1). The maximum correlation between sea level anomalies reconstructed from CEOF1 and CEOF2 is reached (r = 0.3 between 2°S–2°N and 170–120°W) at a positive lag of 14 days [Principal Component (PC) 1 lags PC2], corresponding to about a quarter of the main period of IKWs, which confirms that the two modes capture the same sequences of IKWs. In order to diagnose the seasonality of the PC timeseries, their climatology and climatological variance are shown in Figures 1C,D, respectively, which indicates they both have a peak amplitude during the boreal summer and fall, consistent with the expectation that they grasp the same seasonally synchronized interaction mode between IKWs and TIWs. Lastly, the correlation between the 3-month running variance of PC1 and PC2 reaches 0.84 which confirms that the two modes grasp the same physical mode of interaction between IKWs and TIWs.

In order to document further the relationship between IKWs and TIWs, we perform a lead-lag composite analysis of strong events identified from the SLA timeseries reconstructed from CEOF2 and averaged in the equatorial region (2°S–2°N, 180–80°W). This choice is motivated by the fact that (i) CEOF2 grasps twice as much variance along the equator (see Table 1) and therefore can be interpreted as the portion of TIWs variability that relates to IKWs, and (ii) PC2 leads PC1 slightly (by ∼15 days). 25 downwelling and 28 upwelling events are identified for which the equatorial SLA reconstructed from CEOF2 is above or below 2 standard deviations, respectively (Figure 2A). Figure 2 shows the composite evolution of the dominant EOF mode of equatorial zonal wind stress intraseasonal anomalies in the equatorial band (Figure 2B, see also its spatial pattern and temporal evolution in Supplementary Figures 4, 5 respectively), the equatorial Kelvin wave contribution to sea level anomalies, and the reconstructed sea level anomalies from CEOF2 (Figures 2C). In particular, Figure 2B,C indicate that the positive equatorial SLA events in the central Pacific characterize strong downwelling IKWs triggered by westerly wind bursts that occurred on average ∼10 days prior in the Western Pacific. These westerly wind anomalies spread toward the central-western Pacific and maintain the eastward propagation of IKWs, which take approximately 30 days to cross the entire basin, consistently with the IKWs phase speed. Similar composite patterns but with a reverse sign are observed during upwelling events sparked by easterly wind bursts (Supplementary Figure 6). Interestingly, most of the events used to conduct this lag-composite analysis can also be found in the upwelling IKWs either preceding or following the composite downwelling IKWs presented in Figure 2. This means that the basin scale CEOF analysis, which inherently encompasses most IKWs propagation features into a single mode does not allow disentangling quantitatively the IKWs phases but rather represents the succession of upwelling and downwelling pulses as an integrated pattern of IKWs variability. This is why we observe active TIWs in the northeastern tropical Pacific throughout the entire IKWs propagation cycle across the basin (Figure 2C and Supplementary Figure 6B). Their presence in the eastern Pacific, whereas the IKWs has yet to reach the area, has likely been triggered by the preceding IKWs (∼60 days prior).

Although the CEOF analysis is not designed to grasp asymmetry in the variability, still, we can observe that the amplitude of the CEOF2 SLA as well as the SST and meridional geostrophic velocity anomalies composite (superimposed on Figure 2D and Supplementary Figure 6C) in the region of TIWs activity is larger during the upwelling than the downwelling IKWs propagation cycle. In order to further quantify such asymmetrical response in TIWs amplitude to the IKWs phase, we present in Figure 3, a bar plot between total SLA averaged in the equatorial region (2°S–2°N, 150–80°W) and TIWs amplitude calculated from SLA using Boucharel and Jin’s (2020) method for the strong upwelling and downwelling events identified previously. Figure 3 indicates that upwelling IKWs tend to promote a larger TIWs amplitude than downwelling IKWs by ∼15%. Note that this is significant at the ∼90% level based on a Student’s t-test. This is consistent with the theoretical and modeling study by HT16 who also evidenced an asymmetry in the amplitude of the TIWs activity between the passage of upwelling and downwelling IKWs, although in their case the TIWs response is somewhat stronger (in their experiments, the TIWs kinetic energy is decreased (increased) by 38% (42%) for downwelling (upwelling) IKWs compared to the background control). In addition, in their idealized setting, HT16 estimated a ∼60 days lagged TIWs response to the passage of IKWs, while we find a faster response (∼10–15 days). These differences can be attributed to the distinct methodological approaches between the two studies. In particular, while we diagnose the lag in SLA amplitude between the equatorial and TIWs regions, HT16 use the rate of volume averaged TIWs kinetic energy, which has a longer time scale adjustment due to its quadratic formulation and may encompass the contribution of high-order baroclinic modes due to the depth integration. In addition, unlike in HT16’s idealized settings, the different Kelvin waves phases are not well differentiated in our CEOF observational analysis as most of the time they form a sequence a 2–3 pulses of waves (i.e., wave trains). So, it is expected that the lagged relationship between TIWs and IKWs in the data differs from that in their modeling study, because of the cumulative effect of consecutive IKWs on TIWs activity.

The composite patterns of sea level, SST and meridional geostrophic velocity anomalies at the peak of the passing IKWs (Figure 4) in the region of TIWs activity not only confirms the somewhat asymmetrical response in the TIWs amplitude but also showcase an in-phase signal between temperature and meridional currents anomalies likely leading to a TIWs-induced meridional heat flux convergence (Xue et al., 2020). This now raises the question of the contribution of these coupled IKWs-TIWs individual events to the Cold Tongue heat budget and in particular, to which extent they compare to the strong TIWs-induced non-linear dynamical heating evidenced at seasonal and interannual time scales in previous studies (e.g., Wang and McPhaden, 1999; Menkes et al., 2006; An, 2008; Xue et al., 2020, 2021).

While most focus has been drawn to TIWs-induced heat flux convergence associated with interannual fluctuations in the background state (i.e., ENSO) in earlier studies, our results call for the evaluation the TIWs-induced non-linear dynamical heating (NDH) at intraseasonal time scales. The seasonally varying TIWs has been shown to warm the equatorial cold tongue through meridional heat advection (Wang and McPhaden, 1999). Similarly, many studies have highlighted the TIWs effect at interannual time scales as a potential source for ENSO asymmetry (An, 2008; Boucharel and Jin, 2020; Xue et al., 2020, 2021) due to the rectified effect of anomalous temperature transport by anomalous currents referred to as oceanic eddy heat flux (e.g., Baturin and Niiler, 1997; Menkes et al., 2006; Jochum and Murtugudde, 2006; Jochum et al., 2007, 2008; Graham, 2014; Xue et al., 2020). In this section, we estimate the eddy heat flux spatial derivative composite of the specific events that correspond to periods of strong activity of the IKWs-TIWs interaction mode. This allows diagnosing the NDH of the eastern equatorial Pacific triggered by the passage of strong downwelling and upwelling IKWs. The intraseasonal NDH is approximated here by NDH_TIWs = v′∂yT′ following Xue et al. (2020).

Figure 5 shows the TIWs-induced intraseasonal meridional NDH composite associated with the upwelling (a) and downwelling (b) phase of the IKWs-TIWs interaction mode. They both exhibit a warming pattern located just north of the equator along with a cooling pattern centered around 4°N indicating a rather similar TIWs-induced entrainment of warm equatorial water into the Cold Tongue region between the two IKWs phases. As an integrated measure of the warming effect of the IKWs-TIWs interaction mode on the eastern equatorial Pacific, we present in Figure 5C the intraseasonal NDH longitudinally averaged between 150°W and 110°W (the same region as in Xue et al., 2020). It confirms the rather similar magnitude of NDH_TIWs (∼0.8°C/month) during both the downwelling and upwelling IKWs phases, albeit slightly stronger during upwelling events due to the stronger TIWs activity (Figure 3). Note however that there is a significant spread (dotted line in Figure 5C) amongst the individual events, which reflects a large diversity in the response of NDH_TIWs to individual IKWs or IKWs trains. The spread tends to be larger among the upwelling events along the equator (see dashed turquoise thin line in Figure 5C). What is also noticeable from Figure 5C is that NDH_TIWs falls within the range of warming induced by the interannual eddy heat flux convergence as indicated by the latitudinal profiles of the El Niño and La Niña meridional NDH composites during their developing phase (i.e., July–November average) over the same period (1993–2018). Note that the latter is mostly induced by the prolonged La Niña event (1998–2001) that followed the strong Eastern Pacific El Niño of 1997/98 (not shown). Although both triggered by TIWs-induced NDH, the Cold Tongue warming presumably originates from different dynamical mechanisms.