1 Introduction

Although it develops in the tropical Pacific Ocean, the El Niño-Southern Oscillation (ENSO) has far-reaching impacts on climate and society around the world and is the dominant source of interannual climate variability (Ropelewski and Halpert 1987; McPhaden et al. 2006). Effectively simulating and predicting the ENSO is therefore of great importance. However, Barnston et al. (2012) suggested that real-time ENSO prediction skills declined during 2002–2011 compared with those of the 1980s and 1990s, even though the models used then were less advanced. Climate model predictions largely depend on the initial conditions of the ocean. They are also closely related to the interdecadal variation in ENSO characteristics (e.g., McPhaden 2012; Zheng et al. 2016). For example, England et al. (2014) suggested that the pronounced strengthening in Pacific trade winds around 2000 induced a La Niña type of background state anomaly; that is, the sea surface temperature (SST) became warmer in the west and cooler in the east. Xiang et al. (2013) argued that such changes in the Pacific may favor the generation of the central Pacific (CP) type of El Niño, whose major warming center is in the CP. This indicates that the change in the air–sea coupled system in the tropical Pacific has made the CP type of ENSO occur more frequently since the beginning of the twenty-first century (Lee and McPhaden 2010).

In the CP region, the role of the thermocline in changing SST is limited because the mean thermocline depth (TCD) and mixed layer are relatively deep there. Rather, the zonal advection plays a more important role due to the large zonal gradient of the background SST (Kug et al. 2009; Yeh et al. 2014; Fang and Mu 2018; Fang and Zheng 2018). In the eastern Pacific (EP) region, the vertical processes related to the TCD are dominant (Yang et al. 2020). This means that the influence of the thermocline feedback and its related warm water volume (WWV) variation, which can be expressed as the equatorial mean TCD anomaly in the Pacific, weakens in the CP type of ENSO and has been weak in the period since 2000 (Kug et al. 2009). Since the CP type of ENSO is generally weaker than the EP type (Zheng et al. 2014; Fang et al. 2015), it is accordingly less predictable (McPhaden 2012; Zheng et al. 2016). To summarize, the change in the air–sea coupled system in the tropical Pacific has led to significant changes in the ENSO’s characteristics and predictability.

Another problem in ENSO prediction is the so-called spring predictability barrier (SPB). That is, regardless of when a prediction is made, the prediction skill will drop so seriously when it crosses the boreal spring as to make the following prediction almost meaningless (e.g., Webster and Yang 1992; McPhaden 2003; Mu et al. 2007a, b; Zheng and Zhu 2010a, b). As a result, if we cannot sufficiently resolve the SPB problem, ENSO prediction with a long lead time (n.b., without considering the specific starting month) will be impossible. In this work, four physically oriented variables are firstly extracted to demonstrate the close connection between the spring information and the following ENSO evolution. Then, the influence of the change in the air–sea coupled system in the tropical Pacific is also evaluated through comparison analyses of the periods before and after 2000.

2 Data and methods

2.1 Data

The observational and reanalysis datasets are obtained from the National Centers for Environmental Prediction Global Ocean Data Assimilation System (Behringer and Xue 2004). The datasets are analyzed for the period 1980–2018, and the mean seasonal cycles and linear trends are removed to calculate the anomalies in each field. The tropical TCD is defined as the depth at which the oceanic potential temperature is 20 °C. It should be mentioned that the focus of this article is a relative shorter period (i.e., 1980–2018), for which the CP El Niño events occurred more frequently during 2000–2018 than 1980–1999. For a long period, the CP type of El Niño occurs more frequently since 1978 than that during 1901–1978 (Wang et al. 2019).

2.2 Physically oriented variables

Based on the classical recharge oscillator theory and other efficient statistical prediction models, the WWV (called TCDa_M in the following equations for consistency) and the zonal wind stress anomaly in the western Pacific (called Tauxa_W) are effective at depicting ENSO development (Jin 1997; Clarke and van Gorder 2001). Therefore, these two variables are also adopted in this study.

Furthermore, Bunge and Clarke (2014) suggested that the TCD anomalies in the tropical Pacific have two main variability modes. One is the mode that shows one sign across the Pacific, which has been associated with the WWV. Another is the “tilt” mode, which shows opposite signs in the eastern and western equatorial Pacific. This tilt mode is closely related with the simultaneous Niño3.4 index, so it can reflect the persistence of the ENSO signal. In the meantime, Lai et al. (2018) also emphasized the effectiveness of the TCD anomalies in the western Pacific, rather than the equatorial mean, for ENSO prediction in the boreal spring. As a result, the third variable in this study is chosen to be the zonal gradient of the TCD anomalies in the equatorial Pacific, which is called TCDa_G.

Apart from the dominant zonal processes that influence ENSO evolution, meridional processes have attracted increasing attention in studying ENSO diversity and complexity. For example, Xie et al. (2018) indicated that the meridional air–sea coupling processes in the eastern equatorial Pacific have an important role in ENSO termination. They argued that southeasterly cross-equatorial wind anomalies in the boreal spring cause a moderate El Niño to dissipate rapidly, as they intensify ocean upwelling south of the equator. Consequently, the mean meridional wind stress anomalies over the eastern equatorial Pacific is chosen as the last variable to reflect the meridional processes mentioned above and are called Tauya_E.

In this study, the TCDa_M index is calculated as the mean TCD anomalies in the region (120° E–80° W, 2° S–2° N). The Tauxa_W index is the mean zonal wind stress anomaly in the region (120° E–160° W, 2° S–2° N). The TCDa_G index is the difference between the mean TCD anomalies in the regions (160° W–80° W, 2° S–2° N) and (120° E–160° W, 2° S–2° N), the zonal ranges of which are based on Bunge and Clarke (2014) and Lai et al. (2018). Tauya_E is defined as the mean meridional wind stress anomaly in the region (120° W–80° W, 2° S–2° N), the zonal range of which is based on Xie et al. (2018) and Hu and Fedorov (2018). The main results are not sensitive to a relatively wider meridional range (i.e., 5° S–5° N). The relatively narrow meridional range of the variables is selected to more precisely express the equatorial processes, since they are the core of ENSO mechanisms and the oceanic Rossby deformation radius is approximately 100–200 km in the tropics.

It should be stated that the choice of these variables is mainly based on the physical interpretation but not their direct linear correlation with the Niño3.4 (170° W–120° W, 5° S–5° N) index, since the Bjerknes feedback is so dominant that many high Niño3.4 index-correlated variables may not provide additional assistance to TCDa_M and Tauxa_W, i.e., the two widely used predictors for ENSO prediction with the support of the classical theories. It can be seen in Table 1 that TCDa_M and Tauxa_W highly correlate with the October–December (OND) mean Niño3.4 index, which is relevant with our expectations. The other two variables (TCDa_G and Tauya_E), however, just show little correlations, indicating that they may be negligible for the following ENSO evolution.

Table 1 Correlation coefficients (CC) between the four spring (March–April–May) mean variables (i.e., TCDa_M, Tauxa_W, TCDa_G and Tauya_E) and the OND mean Niño3.4 index
Full size table

However, when we calculate the partial correlation coefficient (PCC), which measures the degree of association between two variables, with the effect of other variables removed (Table 2), it shows that except for TCDa_G, the rest three variables can all play important roles for the following ENSO evolution. This means that Tauya_E could provide some new information differing from the TCDa_M and Tauxa_W indices for capturing the ENSO evolution. For the TCDa_G, although its contribution to the Nino3.4 index cannot pass the 5% significance level using a Student’s t test, we still retain it because it will be shown in Fig. 6 that TCDa_G is significantly important for the far eastern Pacific region during 1980–1999. And this region is crucial for distinguishing the two types of ENSO (Kao and Yu 2009; Kug et al. 2009). In addition, it suggests that TCDa_G and Tauxa_W are not totally independent, i.e., they are significantly correlated with each other. Physically, this relationship mainly reflects the response of the atmosphere to the SST variations over the central to eastern Pacific, i.e., one important link of the Bjerknes positive feedback (Jin 1997). However, since the main response region of the atmosphere does not coincide with the Tauxa_W, i.e., the whole far western Pacific, considering these two variables could bring more information than any one of them.

Table 2 Same as Table 1, but for the PCC, which measures the degree of association between two variables, with the effect of other variables removed
Full size table

2.3 Quaternary linear regression equation

In this study, a quaternary linear regression equation is constructed to measure the relationship between the spring air–sea information and the following ENSO property as follows:

$${Nino}_{p}^{OND}=a\times {TCDa\_M}_{o}^{Mar}+b\times {TCDa\_G}_{o}^{Mar}+c\times {Tauxa\_W}_{o}^{Mar}+d\times {Tauya\_E}_{o}^{Mar}+e,$$
(1)

where \({Nino}_{p}^{OND}\) is the Niño index at the end of the year (i.e., OND mean) and \({TCDa\_M}_{o}^{Mar}\), \({Tauxa\_W}_{o}^{Mar}\), \({TCDa\_G}_{o}^{Mar}\) and \({Tauya\_E}_{o}^{Mar}\) are the TCDa_M, Tauxa_W, TCDa_G and Tauya_E indices in March, respectively. \(a\), \(b\), \(c\), \(d\) and \(e\) are the regression coefficients. As a result, this equation can reflect both the sign and magnitude of ENSO in OND.

3 Results

3.1 Connection between the March information and the following ENSO evolution during the period 1980–2018

As demonstrated in Clarke and Zhang (2019), the seasonal phase locking of ENSO is fundamental to its development and prediction. To quantitatively illustrate this characteristic, the empirical orthogonal function (EOF) analysis is conducted to the seasonal variation (i.e., from May to the following April for each year) of the observed Niño3.4 index, that is, the typical index to denote the strength of ENSO. The results are shown in Fig. 1. It can be seen that the leading mode, which explains nearly 90% of the total variance, exhibits a typical ENSO evolution, i.e., it always initiates in spring, grows up to mature in winter and quickly decays in the following spring. This is consistent with the composite analyses of Rasmusson and Carpenter (1982). The corresponding principle component of the leading mode (i.e., PC1) and their reconstructed time series (i.e., N34_cal; PC1 multiplies the leading mode) are shown in Fig. 1b, c. It can be seen that the reconstructed Niño3.4 index is nearly identical to the original one, i.e., the correlation coefficient between these two time-series is 0.95. Furthermore, Fig. 1b shows that the PC1 is almost identical to the OND mean Niño3.4 index, i.e., the mature amplitude of ENSO. Specifically, their correlation coefficient is 0.99. This tells us that we can potentially provide a quite excellent ENSO prediction for a year ranging May–next April if the OND mean Niño3.4 index can be accurately obtained in advance, with the support of the steady seasonal phase locking of ENSO.

Fig. 1
figure 1

a The leading EOF mode of the seasonal Niño3.4 index variation, i.e., from May to the following April. The number at the top of this panel is the percentage of the variance explained by this mode. b Time series of the normalized principle component of the leading mode (i.e., PC1; black curve) and the OND mean Niño3.4 index (red curve). The correlation coefficient between these two time-series is 0.99. c Monthly time series of the original Niño3.4 index (N34_ori; black curve) and the reconstructed one (N34_cal; red curve), i.e., multiplying the leading mode with the PC1 of the specific year. Both time series range from May 1980 to April 2018. This analysis is similar with Clarke and Zhang (2019), who referred to the months April–March as their “El Niño year”

Full size image

Then, based on the method described in Sect. 2, we first construct an equation for the period 1980–2018 to illustrate the connection between the March information and the following ENSO evolution. Figure 2 shows the reconstructed OND mean Niño3.4 index by utilizing the four variables from March. The calculation suggests that the correlation coefficient between the reconstructed and observed indices is 0.87. This indicates that the four physically oriented variables selected in this study are reasonable. Additionally, the equation can accurately capture La Niña events, i.e., the blue dots are all in the third quadrant. For El Niño, it can capture most events but occasionally fails. This is in accordance with the asymmetry of the two phases of ENSO, i.e., the pattern of La Niña is relatively identical, while El Niño shows more diversity (Kug et al. 2009). Moreover, the coefficients of the equation suggest that the TCDa_M and Tauxa_W indices play dominant roles in ENSO development and are significant at the 5% significance level based on Student’s t test. Tauya_E is the main attenuating factor, illustrating the importance of the meridional processes in the eastern Pacific region, which were paid less attention in classical ENSO theory, to modulating the following ENSO evolution starting in March (Fang and Xie 2020). The TCDa_G index, however, does not pass the significance test, meaning that it might be not an effective factor on affecting the ENSO variations over the whole period 1980–2018. Recently, Wang et al. (2020) also found that simply using SST and zonal winds before April can successfully predict different types of El Niño, which is consistent with our viewpoint. Besides, the experiment with the preceding ENSO signal removed from the original variables is also conducted, the result of which is nearly identical with Fig. 2, except for a little quantitative difference for the regression coefficients (figure not shown). This means that the close relationship between the spring air–sea information and the following ENSO evolution is not significantly influenced by the preceding ENSO signals.

Fig. 2
figure 2

The upper panel shows the relationships between the reconstructed and observed OND mean Niño3.4 index during 1980–2018. The reconstructed indices are obtained by the quaternary linear regression equation utilizing the TCDa_M, TCDa_G, Tauxa_W and Tauya_E indices in March. The lower panel shows the coefficients (bars) and their 5% significance interval by Student’s t test (error bars) of the four variables, which are obtained by the product of the regression coefficient of each variable and its corresponding standard deviation. The red, orange and blue dots represent EP El Niño, CP El Niño and La Niña events, respectively, whereas the black dots represent neutral years. The correlation coefficient (R) between the reconstructed and observed Niño3.4 indices is also shown in the panel, which passes the 5% significance level based on Student’s t test

Full size image

Further analysis shows that the constructed quaternary equation is more effective for the period after 2000 than for the period before 2000. Specifically, the correlation coefficients between the reconstructed and observed Niño3.4 index are 0.92 and 0.84 for 2000–2018 and 1980–1999, respectively. This reflects the interdecadal variation of the ENSO evolution and its characteristics.

It should be stressed that the purpose of this study is not to construct an ENSO forecast model for practical service, but to investigate the relationship between the spring air–sea information and the following ENSO evolution. That is why we perform the linear regression for the whole period to demonstrate that there is an obvious relationship between the air–sea state in spring and the ENSO mature phase. Before performing further investigations, a leave-two-out cross-validation is adopted to verify the robustness of the linear regression method by constructing the equation based on random 37 years during 1980–2018 and using it to estimate the OND mean Niño3.4 index of the two excluded years (i.e., totally 741 samples; shown in Fig. 3). It shows that the correlation coefficient between the reconstructed and observed indices is 0.84, which passes the 5% significance level based on a Student’s t test and well verifies the efficiency of the linear equations.

Fig. 3
figure 3

Leave-two-out cross validation tests for the quaternary linear regression equation constructed during 1980–2018, which measures the relationship between March information and the OND mean Niño3.4 index. We constructed the equation based on random 37 years during 1980–2018 and using it to estimate the two excluded years, i.e., totally 741 samples. Each dot in the subplots represents a random sample. The correlation coefficients (R) between the reconstructed and observed Niño indices are shown in the panel, which passes the 5% significance level based on a Student’s t test. The 1:1 line is marked in dashed grey

Full size image

3.2 Influence of the air–sea coupled system change on the ENSO evolution from boreal spring

As stressed by many previous studies, there is a significant coupled system change of the tropical Pacific. Specifically, compared with the period 1980–1999, 2000–2018 shows a pronounced strengthening trade winds and a La Niña type of interdecadal anomalies; that is, the SST (TCD) becomes warmer (deeper) in the west and cooler (shallower) in the east. These mean state changes were reported to significantly influence the ENSO characteristics (e.g., Xiang et al. 2013; England et al. 2014), as discussed above. In addition, the relationship between the WWV, i.e., an important ENSO predictor based on the recharge oscillator theory that emphasizing the oceanic adjustment processes, and the Niño3.4 index experienced great change before and after 2000. During 1980–1999, the peak correlation (0.66, passing the 5% significance level based on Student’s t test) occurs with WWV leading Niño3.4 by 6 months. During 2000–2018, however, the peak correlation (0.57, passing the 5% significance level based on Student’s t test) occurs with WWV leading Niño3.4 by only 3 months, i.e., the lead time and the correlation both declined. This indicates that the recharge oscillator theory, which dominated the ENSO evolution before 2000, has played a less important role in recent decades (Kug et al. 2009; McPhaden 2012; Ren and Jin 2013; Kumar and Hu 2014a, b; Bunge and Clarke 2014).

Therefore, to utilize the WWV during 1980–1999 as much as possible, the late spring (i.e., May) information is adopted for this period. For 2000–2018, March information is still utilized, just as in the method used in Sect. 3.1 but with a shorter period. That is, separate linear regressions are performed for the two subperiods to study the relationship between the spring air–sea information and the ENSO mature phase. Figure 4 shows the results.

Fig. 4
figure 4

Same as Fig. 2, but separately considering the two periods, i.e., 1980–1999 and 2000–2018. Specifically, during 1980–1999 (2000–2018), the reconstructed OND mean Niño3.4 indices are obtained by the quaternary linear regression equation that uses the TCDa_M, TCDa_G, Tauxa_W and Tauya_E indices in May (March). The blue, red and yellow bars in the right panels represent the periods 1980–2018, 1980–1999 and 2000–2018, respectively. The correlation coefficient (R) between the reconstructed and observed Niño3.4 indices are also shown in the panel, which passes the 5% significance level based on Student’s t test

Full size image

Compared with the whole-period case (Fig. 2), this separate consideration further improves the model performance in capturing the ENSO evolution from the boreal spring. Specifically, the correlation coefficients between the reconstructed and observed OND mean Niño3.4 index during 1980–2018 are as high as 0.91. And all the ENSO events are accurately captured, i.e., the blue (red and orange) dots are in the third (first) quadrant. Based on 1000 Monte Carlo randomly selected samples, these improvements are all statistically significant, details for which can be found in Fig. S1 in the supporting information. The coefficients of the equations also indicate that the equations for the period 2000–2018 are similar to those for the period 1980–2018. However, during 1980–1999, the SST variations in the eastern Pacific were mainly related to the TCD variations, especially the TCDa_M index. This is in accordance with the classical ENSO theory and the relatively high ENSO prediction skills in this period (Barnston et al. 2012). The wind stress, on the contrary, is relatively less important. The experiments just based on TCDa_M and Tauxa_W, i.e., the two key variables in our method and classical ENSO theories, also agree with the main conclusions mentioned above (Fig. S2 in the supporting information).

Further analysis indicates that this better performance is basically due to the improvement of the former period, i.e., the correlation coefficient improves to 0.91 from 0.84 in the whole-period case. This verifies that with the dominance of the recharge oscillator theory in the years before 2000, efficiently utilizing the TCD variation can largely capture the following ENSO evolution. While in the years after 2000, correctly addressing the March information can result in a similarly high capability. Again, the leave-two-out cross-validations for the periods 1980–1999 and 2000–2018 both verify the robustness of the linear regression method (Fig. 5). They also show that the performance is much better when utilizing the May information than when using the March information during 1980–1999.

Fig. 5
figure 5

Same as Fig. 3 but conducting the leave-two-out cross validation tests for the quaternary linear regression equation constructed during 1980–1999 and 2000–2018. For the former period, the equations that utilize the March (a) and May information (b) are both tested. For the latter period, the equation that utilizes the March information is tested. The correlation coefficients (R) between the reconstructed and observed Niño indices are also shown in each panel. They all pass the 5% significance level based on a Student’s t test. The 1:1 line is marked in dashed grey

Full size image

To further illustrate these interdecadal differences, Fig. 6 illustrates the contributions of each variable and their combination to the OND mean SST variations. During 1980–1999, the SST variations explained by the four variables (Fig. 6i) exhibit a typical EP El Niño pattern, i.e., the maximum SST anomalies are in the EP region. This is in accordance with the fact that the EP type of ENSO is prevalent during this period. More specifically, the SST variations over the central to eastern Pacific are mainly related to the TCDa_M and TCDa_G indices, i.e., they can together explain ~ 1.5 °C. The Tauxa_W index plays a moderately positive role in the SST variations in the central Pacific. The Tauya_E index has almost no influence at all.

Fig. 6
figure 6

Contributions of each variable and their combination to the OND mean SST variations, which are obtained by the product of the regression coefficient and the corresponding standard deviation of the SST anomalies in each grid. The unit is °C. Only the coefficients that pass the 5% significance level based on Student’s t test are shown. The left and right panels are for 1980–1999 and 2000–2018, respectively. a, b Are for TCDa_M, c, d are for TCDa_G, e, f are for Tauxa_W, g, h are for Tauya_E, i, j are for their total contribution

Full size image

Nevertheless, the situation is quite different for the period 2000–2018. First, the total SST variation explained by the four variables (Fig. 6j) shows a pattern that is more like the mixed EP and CP El Niños, i.e., the main SST deviations (relatively weak) are in a wide range over the central to eastern Pacific. This also agrees with the fact that during this period, the CP type of ENSO, whose maximum SST anomalies are confined to the CP region and that is generally weaker than the EP type of ENSO, is more prevalent. Furthermore, the SST variations over the central to eastern Pacific are mainly related to the Tauxa_W and Tauya_E indices, while the importance of the TCD variation decreases greatly. The TCDa_G index has almost no influence at all. This indicates that during this period, the processes that classical ENSO theory does not fully cover are very important, e.g., the latent heat feedback related to the thermodynamics of the mixed layer and the meridional air–sea interaction in the eastern Pacific region, as suggested by previous studies on ENSO dynamics (e.g., Yu and Kim 2011; Zheng et al. 2014; Hu et al. 2017; Xie et al. 2018; Hu and Fedorov 2018). This also reminds us that if we forecast the ENSO events after 2000 by using the ENSO models that are constructed based on the air–sea interaction before 2000, i.e., emphasizing the dominant role played by the WWV index, there will be inevitable drawbacks.

It was documented that the ENSO spring persistence barrier, which is tightly related to the SPB, became stronger after 2000 compared to that in the previous decades (e.g., Fang et al. 2019; Liu et al. 2019). The main reason for this is the weakening of the memory of the ocean, as reflected in both SST and TCD. This finding is also consistent with our results.

4 Conclusions and discussion

Since the ENSO evolution shows a strong seasonal phase locking characteristic, i.e., the leading EOF mode of seasonal Niño3.4 variation can explain ~ 90% of the total variance, and the principle component of the leading mode is almost identical to the OND mean Niño3.4 index (i.e., the mature amplitude of ENSO), this means that a quite excellent ENSO prediction for a year ranging May-next April could be achieved if the OND mean Niño3.4 index can be accurately obtained in advance. So, in this study, four physically oriented variables are extracted and used to construct linear regression equations to study the connections between the spring information and the following ENSO evolution. The results indicate that utilizing the TCDa_M, TCDa_G, Tauxa_W, and Tauya_E indices in March can largely capture the mature phase of ENSO. However, whether the four variables are the optimal variables needs further investigation. In addition, this does not indicate that the four variables are totally independent, since all the variables in the model are tightly connected (Tables 1, 2). Further investigation indicates that the specific equation is modulated by the change in the air–sea coupled system. During 1980–1999, the SST variations over the central to eastern Pacific were mainly related to the TCD_M and Tauxa_W variations. This is in accordance with the fact that the EP type of ENSO was more prevalent and the classic recharge oscillator theory was effective on capturing the ENSO evolution during this period. While during 2000–2018, when zonal advection and thermodynamics were more important, it is essential to make good use of the wind stress information and correctly depict their controlled processes. In particular, the meridional processes in the eastern Pacific region, which were less considered in classical ENSO theory, were verified to be important in depicting the following ENSO evolution from March. In this work, the reference climatology is calculated based on the whole analysis period, i.e., 1980–2018. We also performed the analyses based on their own mean states when dealing with the 1980–1999 and 2000–2018 cases, which does not change the main conclusions. This means that the interdecadal change of the air–sea coupled system over the tropical Pacific mainly reflects on the specific relationship between the spring air–sea variables and the OND mean Niño3.4 index, rather than directly on the relevant amplitude of the variables.

These results not only verify the efficiency of the four variables, but also indicate that correctly accounting for the coupled system change and understanding the role of the dominant processes in spring in modulating the following ENSO evolution could reduce the impact of SPB and improve ENSO prediction (Hu et al. 2019). Nevertheless, the purpose of this study is not to construct an ENSO forecast model for practical service, but to investigate the relationship between the spring air–sea information and the following ENSO evolution based on the observational datasets. Through simple linear regression analyses and investigating the interdecadal variation of the specific relationship, we mainly emphasize importance of the physical processes, especially the meridional processes over the eastern equatorial Pacific, which was paid less attention by the classical ENSO theories. In reality, we cannot anticipate the interdecadal shift of the ENSO regime, thus it is a big challenge for the ENSO prediction models to make the prediction using different forecast schemes for different periods. We also intend to investigate how the models of the Coupled Model Intercomparison Project (CMIP) models perform on capturing these close relationships, since there are large biases on depicting the mean state and the seasonal phase locking characteristic of ENSO.

The mechanism of the change in the WWV lead time after 2000 has been explored by previous studies. For example, Horll et al. (2012) and Lübbecke and McPhaden (2014) speculated that the breakdown of the relationship between WWV and ENSO may relate to a shift towards more central Pacific versus eastern Pacific El Niños in recent decades. Hu et al. (2013, 2020) and Li et al. (2019, 2020) suggested it is associated with variability decrease, frequency increase, and westward shift of the atmosphere–ocean coupling center that coincided with the fact that the zonal advection feedback increased while the thermocline feedback declined after 2000. Hu et al. (2017) argued that the shortening of the lead time of WWV to ENSO since 2000 is due to their frequencies increment. And our results also supported these physical explanations.

It should be noted that the role played by the atmospheric variables in influencing the ENSO evolution is not through their own persistent development, but through their influence on the oceanic processes and the subsequent air–sea interactions, i.e., they play roles in triggering the following ENSO evolution from the boreal spring (e.g., Clement et al. 2011; Anderson et al. 2013). Additionally, the finding of this study that the importance of the atmospheric processes during the boreal spring is comparative with the oceanic processes suggests that elaborate consideration should be given to the accuracy of the initial atmospheric conditions and their dynamical consistency with oceanic conditions; this could be achieved by applying coupled data assimilation to climate models (e.g., Zheng and Zhu 2010a, b; Lu and Liu 2018). Also, the four identified variables here are confined to the equatorial Pacific, since we aim to reflect the close relationship between the spring air–sea information in the tropical Pacific and the ENSO signals in its mature phase. Recently, Chen et al. (2020) indicated that the impact of the extratropical atmosphere on the ENSO evolution has intensified during the past 2 decades. Furthermore, they proposed a novel approach to combine the tropical preconditions/ocean–atmosphere interaction with extratropical precursors to noticeably increase the ENSO prediction skill beyond the spring predictability barrier. This provides a research direction to further investigate the causality of the four spring variables, i.e., their physical relationships with the preceding extratropical Pacific variations.