Interdecadal wind stress variability over the tropical Pacific causes ENSO diversity in an intermediate coupled model

The role of interdecadal wind stress variability in the genesis of ENSO diversity is examined by using an intermediate coupled model (ICM) in the tropical Pacific; two types of experiments are performed, one with the original ICM, and the other with interdecadal wind stress (τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document}) effect being explicitly represented. The τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document} component is derived from NCEP/NCAR reanalysis dataset as follows. First, the ensemble empirical mode decomposition (EEMD) is used to extract the interdecadal component of wind stress anomalies on about a 10–40 yr timescale. Next, an idealized interdecadal cycle of τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document} is reconstructed by a principal oscillation pattern (POP) analysis based on the EEMD-extracted interdecadal wind component. A 110-yr model integration is then performed by explicitly incorporating the reconstructed τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document} cycle into the ICM. Compared with the regular interannual oscillation in the original ICM, the simulated ENSO events become highly irregular with interdecadal variations in the amplitude and asymmetry when the τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document} effect is included. Especially, the model reproduces two types of El Niño with different spatial distribution and temporal evolution of SST anomalies, namely Eastern-Pacific (EP) and Central-Pacific (CP) types. Further attribution analyses are performed to understand the modulating effects of τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document} in the tropical Pacific using the ocean component of the ICM, forced by the added τinterde\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\uptau }_{interde}$$\end{document} effect. Two different roles of the Interdecadal Pacific Oscillation (IPO) in modulating different types of El Niño are illustrated. On the one hand, the warm phase of IPO favors for the emergence of EP-El Niño events, in association with the initial warming signals occurring in the eastern equatorial Pacific, which are absent in the original ICM. On the other hand, the cold phase of IPO tends to shift the El Niño (i.e., the single type of El Niño in the original ICM) to an CP type, with which the SST anomalies propagate eastward along the equator but cannot extend into the eastern boundary. This simple modeling study highlights the significant contributions of interdecadal wind variability to the genesis of ENSO irregularity and diversity theoretically.


Introduction
The El Niño-Southern Oscillation (ENSO) phenomenon describes basin-wide anomalous warming (El Niño) or cooling (La Niña) of sea surface temperature (SST) in the tropical Pacific, accompanied by far-reaching atmospheric circulation anomalies. ENSO can affect global weather and climate variations of great relevance to society (Ropelewski and Halpert 1987;Philander 1990;Diaz et al. 2001;McPhaden et al. 2006;Yang et al. 2018). Therefore, understanding and predicting ENSO have always been subjected to extensive attention from the scientific community.
ENSO events exhibit significant diversity, including their onset time, amplitude, duration, spatial pattern, and temporal evolution of SST anomalies in the tropical Pacific (Zhang et al. , 2022Xu and Chan 2001;Okumura and Deser 2010;Capotondi 2013;Capotondi et al. 2015;Timmerman et al. 2018;Hu et al. 2020;Feng et al. 2020;Cai et al. 2021;Chen and Li 2021). Although these differences in ENSO evolutions have been known for many years, the reasons for ENSO diversity remain controversial, having attracted a renewed interest since the late 1990s, due to the frequent occurrence of different types of El Niño. One new type of El Niño, referred to as Central-Pacific (CP) El Niño (Kao and Yu 2009), is characterized by maximum warm SST anomalies located in the central equatorial Pacific rather than in the eastern Pacific, compared with the canonical Eastern-Pacific (EP) type. Different terms are also used to describe the CP type, such as dateline El Niño (Larkin and Harrison 2005), El Niño Modoki (Ashok et al. 2007), and warm pool El Niño ), respectively. These two types of El Niño exhibit different impacts on global weather and climate, because atmospheric teleconnection patterns strongly depend on the SST warming center in the equatorial Pacific (Kim et al. 2009;Weng et al. 2009;Wang et al. 2013;Xu et al. 2018). In addition, there are also obvious differences in the forecasting skill between EP and CP-El Niño events, with a general lower skill for the CP-El Niño (Hendon et al. 2009;Jeong et al. 2012;Tian and Duan 2016;Tang et al. 2018;Ren et al. 2019;Tao et al. 2020). Thus, it is important to understand the underlying mechanisms responsible for ENSO diversity so as to improve their predictions.
A large number of research has been conducted to explain the causes of ENSO diversity, although there is still no consensus on this issue. Some studies evaluate the impact of global warming on ENSO characteristics based on various long-term proxy datasets and climate modeling evidence from a perspective of future projections (Li et al. 2011;McGregor et al. 2013;Kim et al. 2014;Cai et al. 2014Cai et al. , 2015Cai et al. , 2021. In terms of the increased occurrence of CP-El Niño events, Ashok et al. (2007) and Yeh et al. (2009) suggested that there is a flattening of the equatorial background thermocline under global warming, which favors the oceanic vertical processes from the thermocline to influence SST more dominantly in the central Pacific. Others emphasize that the modulation of ENSO behavior can arise from natural variability within the climate system (Wittenberg 2009;Newman et al. 2011;Borlace et al. 2013;Russell and Gnanadesikan 2014). For instance, it is shown that ENSO properties can be modulated by multiple processes, including tropical instability waves, freshwater flux, and ocean biological process (Lengaigne et al. 2007;Busalacchi 2008, 2009;Zhang et al. 2019;Tian et al. 2019;Gao et al. 2020;Zhi et al. 2020). By analyzing the physical processes controlling the two types of El Niño, Kug et al. (2009) indicated that zonal advective feedback plays a crucial role for the CP-El Niño, while the EP-El Niño is mainly attributed to the thermocline feedback. Consistently, Hu et al. (2012) suggested that stronger (weaker) and more eastward extended (westward confined) westerly wind preceding El Niño can result in active (weak) air-sea interaction and strong (weak) thermocline feedback, favoring the development of an EP (CP)-El Niño. In addition, Chen et al. (2015) and Fedorov et al. (2015) concluded that westerly wind bursts acting on different ocean thermal states can potentially induce diverse ENSO events. Focusing on the extratropical effect, Yu et al. (2010Yu et al. ( , 2015 emphasized the importance of subtropical Pacific precursors in the generation of CP-El Niño, somewhat similar to the Pacific Meridional Mode (PMM; Chiang and Vimont 2004). They proposed that the SST anomalies associated with CP-El Niño initially appear first in the northeastern subtropical Pacific and then spread toward the central equatorial Pacific, where the local air-sea interactions are triggered to produce an SST warming condition.
Low-frequency climatic variability could also influence ENSO diversity. For instance, Yu et al. (2015) indicated that the phase change in the Atlantic multidecadal oscillation (AMO; Schlesinger and Ramankutty 1994) causes a strengthening of PMM so as to favor the development of more CP-El Niño events since the early 1990s. The Pacific Decadal Oscillation (PDO; Mantua et al. 1997), the leading dominant mode of decadal variability in the North Pacific, can lead to a systematic shift of ENSO behaviors by changing the background state of the tropical Pacific. The influence of these midlatitude decadal variations can be conveyed to the tropics through oceanic water mass and heat transport associated with the subtropical cell and atmospheric teleconnection (Gu and Philander 1997;Zhang et al. 1998;Rothstein et al. 1998;Barnett et al. 1999;Kleeman et al. 1999;Pierce et al. 2000;Nonaka et al. 2002;Ding et al. 2015). However, some studies argued that the decadal variability within the tropical Pacific could modulate ENSO diversity without influences from the midlatitudes (Choi et al. 2009;Ogata et al. 2013;Di Lorenzo et al. 2015). For example, Wang and An (2002) demonstrated that the decadal changes of background equatorial winds could qualitatively reproduce the observed decadal changes of ENSO behavior through numerical experiments with the Zebiak-Cane model (Zebiak and Cane 1987). In addition, the Interdecadal Pacific Oscillation (IPO; Power et al. 1999), exhibiting structures similar to an ENSO mature pattern in the tropical Pacific, can provide different decadal conditions in favor of differently dominated dynamic processes to affect ENSO behavior. In particular, the climate regime shift in late 1990s related to the IPO phase transition may cause ENSO flavors to change from a dominant EP type to a dominant CP type (McPhaden et al. 2011;Xiang et al 2012;Chung and Li 2013;Zhang et al. 2022).
At present, the observational record is still too short to obtain a reliable statistical relationship between the interdecadal variability and ENSO diversity. Coupled models can therefore provide a powerful tool to conduct long-term integrations for examining their relationship. Nevertheless, it is difficult to isolate the relative contributions of interdecadal variability from the tropical Pacific and midlatitude regions when using fully coupled general circulation models (CGCMs). In this study, we adopt an intermediate coupled model (ICM) to explore the influence of interdecadal wind stress variability on the ways ENSO diversity can be generated. Compared to the comprehensive CGCMs, there are several advantages by using this type of simplified models. First, the model region of the ICM is restricted to the tropical Pacific Ocean, thereby excluding the effects out of the tropical Pacific. Second, the simulated ENSO oscillation is quite regular, without obvious interdecadal variability signal internally induced for both SST and surface wind in the ICM. This model can be used to explicitly incorporate the interdecadal wind component whose effects on ENSO diversity can be clearly illustrated. Third, the ICM-based simulations are computationally efficient to run in time integration. With this ICM, we perform long-term simulations by explicitly including the interdecadal wind stress component to reveal its effects on ENSO diversity. Based on historical datasets, the interdecadal wind stress component is reconstructed by statistical methods of ensemble empirical mode decomposition (EEMD) and principal oscillation pattern (POP). This paper is organized as follows. The datasets and methods used in this study are introduced in Sect. 2. The model and experiments are described in Sect. 3. Section 4 examines ENSO diversity in the model simulations when the interdecadal wind effect is explicitly included. Section 5 addresses the generation of ENSO diversity using oceanonly experiments forced by the added interdecadal wind component, which is attributed to the modulation of the IPO. The summary and discussion are given in Sect. 6.

Datasets
Two datasets are used in this study. The monthly surface wind stress field is from the National Centers for Environmental Prediction (NCEP)-National Center for Atmospheric Research (NCAR) reanalysis (Kalnay et al. 1996). The NCEP/NCAR dataset is on a global T62 Gaussian grid with 192 × 94 points, and the period from 1948 to 2020 is chosen in order to represent interdecadal variability as much as possible. The monthly SST field is the extended reconstructed sea surface temperature, version 5 (ERSSTv5; Huang et al. 2017) from the National Oceanic and Atmospheric Administration (NOAA), with a 2° × 2° horizontal resolution from 1950 to 2020.

Ensemble empirical mode decomposition (EEMD)
The surface wind stress anomalies exhibit variability on different timescales. To isolate wind stress signals on different timescales, we adopt the method of ensemble empirical mode decomposition (EEMD; Wu and Huang 2009). The EEMD is a noise-assisted data analysis method, an improved version from the original empirical mode decomposition (EMD; Huang et al. 1998), which overcomes the mode mixing problem. The principle of EEMD is to add white noise into the signal (original data), with which the spectrum of white noise is evenly distributed so that the signal can be automatically distributed to the appropriate reference scale. The added white noise series cancel out each other in the final mean of an ensemble of trials with the corresponding intrinsic mode functions (IMFs). Then the mean IMFs stay within the natural dyadic filter windows and significantly reduce the chance of mode mixing and thus preserve the dyadic property. In brief, the EEMD can be used to decompose a complex signal into a set of IMFs with different frequencies and a residual long-term trend. More details of this method can be found in Wu and Huang (2009). In this study, the EEMD is applied to surface wind stress in the tropical Pacific from the NCEP/NCAR reanalysis dataset during 1948-2020. The added white noise has an amplitude 0.2 times that of the standard deviation of the original data and the ensemble number is 1000. Figure 1a illustrates an example for EEMD-based results for zonal wind stress component at 180°on the equator. The C1-C8 are the first to eighth direct output components ("C1-C8" in  Fig. 1b). In this work, our focus is placed on the interdecadal component of wind stress, referred to as the interdecdal wind stress (denoted as τ interde ) in the following, which are used to construct an idealized interdecadal wind cycle through the principal oscillation pattern (POP) analysis as follows.

Principal oscillation pattern (POP) and reconstructed interdecadal wind cycle
To characterize interdecadal variability, we use the method of principal oscillation pattern (POP; Hasselmann 1988; Von Storch et al. 1988) to reconstruct an idealized interdecadal wind stress cycle. The POP analysis is a technique to extract periodically-varying POP patterns and corresponding temporal coefficients from a multi-component system whose dynamics are too complex to be easily described (Von Storch et al. 1988). A time series X(t) is assumed to be generated by a first-order multivariate autoregressive process: where ξ is Gaussian white noise matrix and constant matrix B can be obtained with lag-1 and lag-0 covariance matrices of X(t). The linear system's normal modes, or called POPs, are eigenvectors of the matrix B, which are usually complex due to the asymmetry of B.
Since we are interested in the interdecadal timescale in this study, the POP analysis is performed on the EEMDextracted interdecadal wind stress ( τ interde ) component in Sect. 2.2. Three significant POP pairs are obtained, accounting for 53.5% of the total variance. Their corresponding periods are 211-month, 319-month and 407month, respectively. The second POP pair with the period of 319-month (~ 26.6-yr) is selected for reconstructing an (1) idealized interdecadal wind stress cycle as follows. Firstly, this complex POP mode pair is used to decompose the wind stress fields into a two-dimensional space spanned by the real and imaginary parts of the POP mode. The evolution of this POP mode is manifested in the form of alternating real ( P r ) and imaginary ( P i ) parts between a quarter cycle. Meanwhile, the temporal coefficients of imaginary POP mode [ Z i (t) ] exceed the real POP mode coefficients [ Z r (t) ] by a quarter cycle. Then, it is assumed that the evolution of POP temporal coefficients conforms to a sinusoidal function, with a period of 319-month and an amplitude of 1, so that the POP cycle is allowed to repeat itself undamped. The real and imaginary coefficients are written as: Thus, the reconstructed wind stress series, X (t) , can be written as Examples of decompositions into temporal components and reconstructions of zonal wind stress at 180° on the equator using ensemble empirical mode decomposition (EEMD) method. Shown in a are for the original time series, its intrinsic mode functions (IMFs; C1-C8) and the residual. Shown in b are reconstructed results from the corresponding components in a, where ''interdecadal'' stands for the combination of C6-C8, which represents the interdecadal signal, a focus in this study Figure 2 shows an example of the reconstructed interdecadal wind stress cycle during 1-50 years. Totally, this reconstructed interdecadal wind stress ( τ interde ) cycle is 110 years long and will be explicitly added into an ICM to perform sensitivity experiments. background stratification and partial nonlinear contributions in the momentum equations so that it can realistically simulate ocean currents in the tropical Pacific. The SST anomaly model is embedded into the ocean dynamical model to describe the evolution of interannual SST variability in the surface mixed layer. One distinguishing feature of the IOCAS ICM is a non-local empirical parameterization for the subsurface entrainment temperature ( T e ), which is developed by using an inverse approach and a relationship between T e and sea level anomalies. The atmospheric component of the IOCAS ICM is a statistical model for wind stress interannual anomalies ( τ inter ), which is constructed by using singular value decomposition (SVD) analysis, representing the interannual anomalies of wind stress in response to given SST anomalies.

The IOCAS ICM
During the coupling procedure, the component models exchange anomaly information every day. The dynamical IOM produces oceanic anomaly fields including currents and sea level fields. Then, the currents are used to calculate the SST anomaly based on its anomaly equation, in which T e anomalies are updated from the sea level anomalies. Given the SST anomalies, the τ inter model is used to calculate the wind stress interannual anomalies, which are then used to force the ocean dynamical model. Further details of the IOCAS ICM can be found in Zhang and Gao (2016a, b).

Experiments
Two simulations are performed using the IOCAS ICM. For the control run, the original IOCAS ICM is integrated for 100 years. Then, another simulation is conducted to investigate the influence of interdecadal wind variability on ENSO; the reconstructed interdecadal wind stress cycle by EEMD and POP is added onto τ inter in the ICM (denoted as τ interde run). The τ interde run is integrated for 110 years, corresponding to the reconstructed interdecadal cycle length. A schematic diagram of the experimental procedure is given in Fig. 3. In the τ interde run, the interdecadal wind effect is purposely added when the wind stress anomalies are used to force the ocean model. Comparisons between these two simulations are made to demonstrate the effects of interdecadal wind variability on ENSO.

ENSO diversity due to the interde effect
As mentioned in the introduction, when the original IOCAS ICM is integrated without the interdecadal wind effects, the model tends to produce regular interannual oscillations. Figure 4 displays the time series of Niño3.4 indices (SST anomalies averaged within 5°S-5°N, 120°W-170°W) and the power spectra of Niño3.4 from the observations and model simulations, respectively. Compared to the observations, the control run gives rise to a quite regular ENSO oscillation, with a dominant 4-yr period (Fig. 4e). Note that the regularity of ENSO is mainly due to the fact that the statistical atmospheric τ inter model used in the ICM is linear, which excludes stochastic atmospheric forcing. Additionally, the 11-yr running mean of Niño3.4 standard deviation and skewness are also shown in Fig. 4 to examine the interdecadal variability of ENSO amplitude and asymmetry. As shown in Fig. 4h and k, there is no obvious interdecadal variation in the simulated ENSO amplitude and asymmetry in the control run.
When the interdecadal wind forcing is explicitly included into the ICM, the simulated ENSO oscillation in the τ interde run becomes irregular compared with that in the control run (Fig. 4c). From Fig. 4e and f, the spectral analysis of Niño3.4 shows that the dominant interannual period of 4-yr is not significantly affected by the interdecadal wind forcing, only with a reduced variance amplitude. This is mainly due to the strong coupling between atmosphere and ocean during ENSO in the IOCAS ICM. Nevertheless, the spectral analysis of the 11-yr running mean of Niño3.4 shows interdecadal variance with a period of 15-30 yr in the τ interde run, while the control run does not show any power on interdecadal time scales. The amplitude of ENSO events changes considerably and exhibits obvious interdecadal changes after including the interdecadal wind forcing (Fig. 4i). The asymmetry between El Niño and La Niña also shows an interdecadal change in the τ interde run, as witnessed by the skewness of Niño3.4 (Fig. 4l). These results suggest that the interdecadal wind variability can have significant modulating effects on the simulated ENSO regularity and amplitude in the IOCAS ICM. , an SST anomaly model, an empirical anomaly model for T e , and a statistical atmospheric model to calculate the interannual anomalies of surface wind stress ( τ inter ). In the control run, the IOCAS ICM is integrated with the τ inter model. In the τ interde run, the IOCAS ICM is integrated with the τ inter model and the prescribed interdecadal wind stress cycle explicitly superposed (i.e., τ inter + τ interde ) To explore the influence of the interdecadal wind variability on ENSO diversity, we classify El Niño events in the τ interde run. Firstly, an El Niño event is selected when the Niño3.4 is greater than 0.5 °C during November-January (NDJ) mature phase. As such, a total of 22 El Niño events are selected, respectively. Furthermore, we classify the 22 El Niño events as EP or CP type by comparing SST anomalies at the Niño3.4 and Niño1 + 2 (10°S-0°, 80°W-90°W). El Niño event with a Niño1 + 2 SSTA being comparable to the Niño3.4 SSTA is considered as an EP-El Niño. El Niño event with a Niño1 + 2 SSTA being smaller than the Niño3.4 SSTA is considered as an CP-El Niño, which means its warming is largely detached from the coast of South America. Following this criterion, 8 El Niño events are classified as EP-El Niño and 11 as CP-El Niño, respectively. It is noted that 3 El Niño events are shown to have mixed features of the EP and CP types and so they are not placed in either of these two types. The classified Niño3.4 and Niño1 + 2 indices for EP and CP-El Niño events in the τ interde run are illustrated in Fig. 5. All the El Niño events reach their mature phase during the boreal winter, showing the phase-locking feature of both EP and CP-El Niños. For the EP-El Niño, the peak intensity of SST anomalies in the Niño1 + 2 regions is as large as that in the Niño3.4 region, while the Niño1 + 2 of CP-El Niño is obviously weaker than the Niño3.4 during the peak season. Previously, the comparison between Niño3 and Niño4 indices is widely used to separate the two types of El Niño. We also made the comparison between Niño3 and Niño4 indices and found that the above classification results are not affected by using the Niño3 and Niño4 indices. Note that the differences in spatial pattern of EP and CP-El Niño in the IOCAS ICM are mainly the magnitude of SST anomalies in the eastern equatorial Pacific near the coast (Fig. 6). Therefore, we still use the comparison of the Niño3.4 and Niño1 + 2 regions to illustrate the results. The composite SST anomalies of the EP and CP-El Niño during the NDJ (Fig. 6) show obvious differences in their spatial distribution and amplitude, generally consistent with previous studies (Ashok et al;Kao and Yu 2009;Kug et al. 2009). In the case of EP-El Niño, large SST anomalies are widely located in the central-eastern equatorial Pacific, extending from the coast of South America to the dateline. As for the CP-El Niño, major SST anomalies are mostly confined west of the 120°W, while those in the 90°W-120°W longitudinal range are substantially weak. Furthermore, the maximum SST anomalies with the EP-El Niño can reach 3 °C while those with the CP-El Niño have about 2 °C, consistent with observations indicating that EP-El Niño is usually stronger than CP-El Niño.
In addition to the differences in the spatial characteristic and amplitude of SST anomalies, the temporal evolutions of EP and CP-El Niño are also quite different. Figure 7 displays the space-time evolutions of composite wind stress and SST anomalies during the development of EP-El Niño and CP-El Niño. For the EP-El Niño, the warm SST anomalies first appear in the eastern coast of the equatorial Pacific during the early spring and then propagate westward to grow into a basin-wide warming in the tropical Pacific (Fig. 7a), accompanied by anomalous westerlies in the central equatorial Pacific. However, the CP-El Niño exhibits quite different evolution in the SST anomalies and westerly wind anomalies. As shown in Fig. 7b, before the onset of CP-El Niño, there exists an obvious pattern Composite horizontal distributions and 5°S-5°N averaged equatorial profiles of SST anomalies during November-January (NDJ) for a EP-El Niño and b CP-El Niño in the τ interde run of warm SST anomalies (which extend from the northeastern Pacific to the central equatorial Pacific) and cold SST anomalies which are located in the eastern equatorial Pacific; this spatial pattern is somewhat similar to the Pacific Meridional Mode (PMM). The subtropical warm SST anomalies persist and gradually spread southwestward to the central equatorial Pacific in April. Then, the CP-El Niño event begins to develop and exhibits a slightly eastward propagation, but large SST anomalies have not extended all the way into the eastern coast region. Corresponding to the initial PMM-like pattern of SST anomalies, an anomaly east-west pattern of zonal wind stress is seen in the central-western equatorial Pacific. On the one hand, the westerly anomalies in the western Pacific help to link the subtropical SST anomalies to the equatorial Pacific and contribute to the subsequent growth of El Niño. On the other hand, the anomalous easterlies can suppress the SST warming in the eastern Pacific by anomalous upwelling in spite of their gradual weakening. To summarize, the development of EP-El Niño is mainly attributed to the air-sea interactions in the equatorial Pacific, exhibiting typical features of the conventional El Niño. However, the onset of CP-El Niño is closely associated with SST variability in the subtropical Pacific.
An SST budget analysis is performed to contrast the contributions of different physical processes to the development of EP and CP-El Niño events. The mixed layer SST anomaly equation of the IOCAS ICM is written as: where overbars and primes indicate monthly climatology and anomaly, respectively. T ′ and T ′ e are anomalies of SST and the subsurface entrainment temperature ( T e ); H is the depth of the mixed layer; H 2 is a constant (125 m); M(x) is the Heaviside step function; T ′ is the parameterized surface heat flux with = (100d) −1 ; the other variables are conventional.
According to the Eq. (5), Fig. 8 shows the composite evolutions of SST tendency, the oceanic advection, the horizontal and vertical diffusion, and the surface heat flux (5) averaged in the Niño3.4 and Niño1 + 2 regions during the EP and CP-El Niño, respectively. The surface heat flux always has a damping effect on the SST anomalies. For the EP-El Niño, the warming tendency of Niño3.4 SST is strongest in September, while the Niño1 + 2 SST shows strong tendency in April. This indicates that the warming signals of EP-El Niño initially appear in the eastern Pacific and extend westward continually. Another positive tendency of Niño1 + 2 SST in November corresponds to the peak phase of EP-El Niño. In the CP-El Niño composites, the largest warming tendency of Niño3.4 occurs in April, corresponding to the extension of the subtropical warm anomalies into the central equatorial Pacific (Fig. 7b).
Another smaller warming tendency is seen to occur in September which indicates a local growth of SST anomalies but with less magnitude during the peak season. For the Niño1 + 2 regions, there is also a large warming tendency in April, corresponding to the weakening of the eastern cold SST anomalies during early spring in Fig. 7b. After May, there is little increase of SST anomalies in the Niño1 + 2 regions of CP-El Niño. The oceanic advection primarily contributes to the SST warming of both EP and CP-El Niño. In addition, large positive vertical diffusion effects appear in the Niño1 + 2 regions especially during Fig. 8 Composite evolution of different terms in the SST budget analyses averaged in the Niño3.4 (left panels) and Niño1 + 2 regions (right panels) for a, b EP-El Niño and c, d CP-El Niño in the τ interde run. The black line indicates the SST tendency term, the red line is for the heat flux, the green and purple lines are for the vertical and horizontal diffusion, and the blue line is for the oceanic advection Fig. 9 Same as Fig. 8 but for the different oceanic advection terms. In each panel, the blue line indicates the zonal advection, the green line is for meridional advection, the red line is for the vertical advection, and the black line is the sum of the three oceanic advection terms the EP-El Niño, which is mainly associated with the positive T e anomalies. The three ocean advection terms are shown in Fig. 9. For the EP-El Niño, the zonal advection is dominant in both Niño3.4 and Niño1 + 2 regions. In the Niño3.4 region, the largest contribution is from the anomalous zonal advection of mean temperature gradient ( −u � T∕ x ), related to the zonal advective feedback (Fig. 10a). Interestingly, the contribution of the zonal mean current associated with anomalous temperature gradient ( −u T � ∕ x ) is also obvious, which is dominant in the Niño1 + 2 regions (Fig. 10b). This may relate to the modulation of the IPO as shown in Sect. 5. During the positive IPO phase, for example, SST warm anomalies occur in the eastern Pacific, which yield a positive temperature gradient and thus contributes to SST warm tendency through the zonal mean westward current. Among the vertical advection terms, the mean upwelling and anomalous temperature gradient ( −w T � ∕ z , i.e., the thermocline feedback) is relatively large as shown in Fig. 10.
For the CP-El Niño, the zonal advection mainly contributes to the SST warming in early phase during January-June (Fig. 9c). Specifically, the −u � T∕ x term is positive but with less magnitude relative to EP-El Niño, while the nonlinear term of zonal advection ( −u � T � ∕ x ) is comparable to the −u � T∕ x term. The positive −u � T � ∕ x effect during Jan-May is due to the anomaly eastward current and negative SST gradient in the central-western Pacific associated with the PMM-like pattern during the onset of CP-El Niño (Fig. 7b). Afterwards, the negative −u � T � ∕ x counteracts the positive contribution of the −u � T∕ x term so that the zonal advection is almost negligible (Figs. 9c and 10c). The thermocline feedback term ( −w T � ∕ z ) in the Niño1 + 2 regions is much weaker during CP-El Niño. In addition, the meridional advection is also large in both types of El Niño; it is dominated by the mean current and the anomalous temperature gradient, which tend to advect the warm SST on the equator poleward and broaden the SST anomalies (figures not shown).
To sum up, the heat budget analysis indicates that the zonal advective feedback and thermocline feedback significantly contribute to growth of SST anomalies in EP-El Niño. In particular, the contribution of the zonal mean current and anomalous temperature gradient ( −u T � ∕ x ) also plays a critical role, which may reflect the modulation of the IPO. As for the CP-El Niño, the zonal advective feedback is dominant, while the nonlinear term of zonal advection shows great influence during the early phase. This is mainly due to the fact that the CP-El Niño events in the τ interde run are closely related to the PMM-like mode.

ENSO diversity associated with Interdecadal Pacific oscillation (IPO)
As demonstrated in Sect. 4, when the interdecadal wind effects are explicitly included in the IOCAS ICM, the simulated ENSO oscillation becomes irregular in terms of amplitude and asymmetry. In particular, the model reproduces two types of El Niño events with different spatial patterns and temporal evolutions. We therefore naturally ask: how the ENSO diversity is generated by the interdecadal wind stress in the IOCAS ICM? To address this question, an additional ocean-only experiment is conducted by using the ocean component of the ICM (i.e., IOM). For Fig. 10 Composite evolution of different zonal and vertical advection terms averaged in the Niño3.4 (left panels) and Niño1 + 2 regions (right panels) for a, b EP-El Niño and c, d CP-El Niño. The blue lines correspond to each zonal advection term and the red lines correspond to each vertical advection term. The dotted lines indicate con-tributions by the mean current and anomalous temperature gradient. The solid lines are for the anomalous current and mean temperature gradient. The dashed lines are for the anomalous current and anomalous temperature gradient this ocean-only experiment with air-sea interactions being switched off, the IOM is directly forced by the reconstructed interdecadal wind stress anomaly cycle from POP. Figure 11 displays the first leading EOF mode and PC time series of the 11-yr running mean SST anomalies obtained from the ocean-only run (i.e., IOM run). It shows that when the interdecadal wind component forces the IOM, an interdecadal mode is generated, which is characterized by the Interdecadal Pacific Oscillation (IPO; Di Lorenzo et al. 2015). The IPO features a broad triangular shape of SST anomaly in the tropical Pacific, resembling an ENSO pattern but with a wider meridional scale (Fig. 11a). The correlation analysis is applied to examine the relationship between the IPO and ENSO amplitude and asymmetry. It is found that the simultaneous correlation between the PC time series of the IPO and the 11-yr running mean of Niño3.4 standard deviation is -0.64. This suggests that ENSO amplitude can be significantly modulated by the low-frequency IPO mode. In addition, the 11-yr running mean of skewness also shows a significant relationship with the IPO, which has a correlation coefficient of 0.85. This in-phase relationship between the skewness and IPO indicates that the warm conditions related to the IPO act to enhance El Niño events, while the cold interdecadal states favor La Niña events. In Fig. 11, we perform the same analysis with the observations during 1950-2020 to examine the interdecadal changes of ENSO amplitude and skewness with respect to the IPO phase. Overall patterns of the IPO in the tropical Pacific are captured well in the IOM run ( Figs. 11a and b). As seen, the IPO has undergone two typical phase transitions occurring around in 1976/77 and 1999/2000 respectively. In observations, there shows little relationship between the IPO and ENSO amplitude and skewness throughout the entire period of 1950-2020. The simultaneous correlation between the PC time series of IPO and the 11-yr running mean of Niño3.4 standard deviation (skewness) is -0.08 (0.04). Nevertheless, we found that when considering time periods before 1980, the relationship between the IPO and ENSO statistics is significantly enhanced, with the correlation between the IPO and amplitude (skewness) being −0.72 (0.22). This suggests that the interdecadal relationship between the IPO and ENSO statics identified in the IOM run can be also found in observations. The modeled high correlation between the IPO and ENSO properties is mostly due to the fact that only single interdecadal signal (i.e., ~ 26.6-yr) is considered in the IOM run and also that the model only resolves the oceanic processes in the tropical Pacific without midlatitude-tropical connections. In nature as observed, on the other hand, ENSO can be modulated by multi-time scale and multi-regional variability, which may lead to the lower correlation compared to the model simulations.
Next, we explore whether the two types of El Niño events are also related to the IPO. As shown in Sect. 4 with the τ interde run, we have identified 8 EP-El Niño and 11 CP-El Niño events, respectively. Here, these events are furtherly separated according to the warm/cold phases of the IPO in the IOM-only simulation. Table 1 shows the numbers of two types of El Niño events occurring during the warm and cold phases of the IPO. It shows that there are 7 EP-El Niño and 2 CP-El Niño events during the warm phase of the IPO, while there are 1 EP-El Niño and 9 CP El Niño during the cold IPO phase. This implies that the warm IPO phase is favorable for the EP-El Niño, while the cold IPO phase favors the occurrence of CP-El Niño. This result is Fig. 11 a, b The first leading EOF patterns of the 11-yr running mean SST anomalies in the tropical Pacific, c, d the 11-yr running mean of Niño3.4 SST standard deviation (black line with left y axis), and e, f the 11-yr moving window of Niño3.4 skewness (black line with left y axis). In c-f, the red line (with right y axis) represents the PC time series corresponding to EOF1 in a and b. Results in left panels are from the IOM ocean-only run which is forced by the reconstructed interdecadal wind stress anomalies. Results in right panels are from observations during 1950-2020 consistent with the observational analysis of Chung and Li (2013), who examined the interdecadal change of the mean state and two types of El Niño around 1999. They clearly demonstrated that when an El Niño-like (La Niña-like) mean state occurs, the interannual SST variability is dominated by EP-El Niño (CP-El Niño).
Physically, it is conceivable that the interdecadal changes in the tropical Pacific associated with the IPO provide a changing background state on which ENSO evolves, which exerts effects on ENSO properties. During the warm phase of the IPO, for example, there are positive SST anomalies in the eastern Pacific (Fig. 11a), which contribute to the emergence of initial warming signals in the eastern equatorial Pacific. When such warm SST anomalies lead to anomalous westerlies along the equator so as to establish the Bjerknes positive feedback, an EP-El Niño event is more likely to develop consequently. Regarding the CP-El Niño, we have to recall for the El Niño events simulated by the original IOCAS ICM (i.e., the control run). Figure 12 shows the composite evolutions of SST anomalies during El Niño events in the control run, and during CP-El Niño events in the τ interde run, respectively. By comparing Fig. 12a and b, it can be clearly seen that the temporal evolution of CP-El Niño identified in the τ interde run is largely similar to that of El Niño in the control run, with one exception being that the CP-El Niño in the τ interde run does not reach the farther eastern boundary regions. Note that a PMM-like pattern exists during the onset period of El Niño events in the control run. The initial PMM-like pattern exhibits warm SST anomalies in the subtropics and cold anomalies in the eastern equatorial Pacific. When such PMM-like condition occurs during the cold phase of IPO, the corresponding patterns of IPO can strengthen the cold SST anomalies in the eastern Pacific and thus anomalous easterlies, which hinder the eastward propagation of warm anomalies to the eastern coast, thereby favoring the development of an CP-El Niño event in the τ interde run. In this sense, the El Niño events in the original IOCAS ICM are changed into CP-types in the τ interde run through the modulation of the IPO. Meanwhile, the IPO can lead to the emergence of new EP-type of El Niño events, which does not exist in the original IOCAS ICM when the interdecadal wind effect is not presented.

Summary and discussion
The role of interdecadal wind variability in the genesis of ENSO diversity is examined by using the IOCAS ICM in which wind stress interdecadal anomalies are explicitly included. The IOCAS ICM by itself does not contain interdecadal signals and thus produces a quite regular interannual oscillation. When the interdecadal wind effects are explicitly included, the simulated interannual oscillation becomes irregular with the amplitude and skewness showing obvious interdecadal changes. In particular, the model reproduces two types of El Niño exhibiting different spatial distributions of SST anomalies. For the EP-El Niño, large SST warm anomalies are located in the central-eastern equatorial Pacific extending from the coast of South America to the dateline. In the case of the CP-El Niño, the SST anomalies are mainly confined west of the 120°W, while those in the eastern coastal regions are substantially weak. In addition, the temporal evolutions of EP and CP-El Niño events also show distinct behaviors. SST anomaly signals during the onset of EP-El Niño are primarily seen in South America coast and subsequently shows a westward extension to the central equatorial Pacific, bearing typical characteristics of conventional ENSO events (Rasmusson and Carpenter 1982). However, the CP-El Niño events tend to show initial SST anomalies in the subtropical Pacific, which gradually spread into the central equatorial Pacific and then propagate eastward along the equator but cannot reach the eastern boundary regions. The SST heat budget analysis indicates that for both EP and CP-El Niños, the ocean dynamical advections contribute to the development of SST anomalies. In the EP-El Niño, both the zonal advective feedback and thermocline feedback are important for the SST growth. Particularly, in the Niño1 + 2 regions, the zonal mean current associated with anomalous temperature gradient is a dominant process, while the thermocline feedback is secondary. This may reflect the modulation of the IPO phase. On the other hand, in the CP-El Niño, the well-known zonal advective feedback significantly contributes to the growth of SST anomalies. Meanwhile, the nonlinear term of zonal advection also plays an important role, which is related to the initial PMM-like mode.
The causes of the ENSO diversity in the IOCAS ICM are explained by the modulation of the IPO, which is realized by the added interdecadal wind forcing in the ocean-only experiment. The IPO shows an ENSO-like SST anomaly pattern but is broader in the meridional direction. On the one hand, the IPO contributes to the emergence of the EP-El Niño events in the τ interde run, which are not reproduced in the control run of the IOCAS ICM. The SST warm anomalies during the positive phase of IPO in the eastern Pacific favor the emergence of initial anomaly signals occurring in the eastern equatorial Pacific, which then induce anomalous westerlies along the equator to develop into an EP-El Niño event through the Bjerknes positive feedback. On the other hand, the IPO tends to change the El Niño simulated in the control run into an CP-type of event in the τ interde run. The IPO can strengthen the eastern SST cold anomalies of the PMM-like pattern during the onset period of the CP-El Niño, which hinder the eastward propagation of SST anomalies to the eastern boundary. To summarize, the interdecadal wind stress effects cause the ENSO diversity in the IOCAS ICM through the modulation of the IPO.
In this study, we emphasize the role of the interdecadal wind stress variability in the genesis of ENSO irregularity and diversity in the IOCAS ICM, which is accomplished by explicitly adding interdecadal wind stress effect into the model. The interdecadal wind stress anomalies are prescribed as an external forcing with no interactive coupling within the ICM model system. With the prescribed interdecadal wind stress, the model shows consistent changes of the wind stress and SST. However, one should note that the model results do not indicate a simple linear cause-andeffect relationship between interdecadal wind variability and ENSO changes. This is because that the wind stress forcing is part of the air-sea coupling system in nature, in which the atmosphere and ocean anomalies are always occurring coherently and consistently due to the coupling (It is not able to identify a simple relationship between the interdecadal wind variations and ENSO changes in the observations as explicitly represented in the model). Therefore, an interactive way is further needed to take into account the air-sea coupling, so as to better represent and understand the relationship between interdecadal variability and ENSO. It is also hopeful that the approach adopted in this paper can be used to improve the ENSO representation and prediction of the IOCAS ICM in terms of interdecadal wind effects. In addition, the results in this work are based on the NCEP/ NCAR reanalysis. Kumar and Hu (2012) indicated that the surface wind stress and convection anomalies associated with ENSO SST in the NECP/NCAR are substantially weaker than other later developed reanalysis. The possible impact of biases in NCEP/NCAR should be identified through comparison by analyzing multiple reanalysis datasets. Furthermore, the results obtained in this study could be model-dependent. The IOCAS ICM consists of an intermediate complexity ocean model based on a baroclinic modal decomposition, and an SVD-based statistical atmospheric model to calculate τ inter . As mentioned earlier, the simulated El Niño events in the control run with the IOCAS ICM tend to show an obvious extension of SST anomalies from the northeastern subtropical Pacific during the onset period. This is quite different from another commonly used intermediate Zebiak-Cane model (Zebiak and Cane 1987). The Zebiak-Cane model is composed of a Gill-type steady state linear atmospheric model and a reduced-gravity oceanic model. The ENSO events in the Zebiak-Cane model often start from the eastern equatorial Pacific and extend westward to cover the central-eastern Pacific. As such, due to the differences between ENSO behaviors simulated in different models, the influences of the interdecadal wind could also be different. Therefore, a comparative modeling study using differently formulated models is necessary to fully understand the modulating effects of interdecadal wind variability on ENSO events.