The recharge oscillator model (Jin, 1997) describes the dynamic relationship between the equatorial Pacific thermocline depth (or OHC) anomaly (H) and the eastern SST anomaly (T) and can be written as (Levine and McPhaden, 2015):

where λ is the damping rate, ω ₀ is the ENSO linear frequency, ξ is white noise (stochastic forcing) and σ is the noise amplitude. Eq. (2.1) indicates the combined effects of the relaxation of SST anomaly (negative feedback), advection, Ekman upwelling, thermocline positive feedback, and noise forcing on the SST anomaly. Eq. (2.2) suggests a description of the basin-wide equatorial oceanic adjustment. Note that the damping term of OHC is neglected here. This is because, on seasonal and longer time scales, changes in OHC are mainly governed by the geostrophic response to the wind stress forcing rather than by the damping itself (Burgers et al., 2005). Here, wind stress can be linearly represented by T and eq. (2.2) is obtained consequently. According to eq. (2.2), there is a phase lag between H and T, suggesting that H can be a predictor for SST at large lead times. Meanwhile, SST can also be a predictor of H changes (see also Figure 2 ).

Unlike H in eq. (2.2) (WWV or WWV_W), the reliable SST data is longer. This encourages us to construct a new predictor using SST data to represent H. Therefore, mathematically motivated by eq. (2.2) and for prediction purposes, H can be calculated as follows:

H(t)≈ω ₀ SST_H (t) and ω ₀ does not change dramatically with time compared with SST anomalies, the lead-lag correlation between H and T (r(τ)) can be obtained as follows:

According to eq. (2.4), r(τ) is independent of ω ₀. As the correlation relationship between H and T is independent of ω ₀, the SST-constructed index SST_H which approximately represents OHC, can be written as:

It can be simply called a temporally lagged SST index. For initial time t = 1, SST_H (1)=T(1). We will show that this index is a predictor of ENSO SST anomalies, especially after the 21^st century.

Two datasets are used in this study. One is the monthly Simple Ocean Data Assimilation (SODA; 0.5°×0.5°; Carton and Giese, 2008). The SODA 3.12.2 data (from 1981 to 2017) is integrated to interpolate the depth of the 20°C isotherm (Z ₂₀), which is used to calculate Z ₂₀This depth of Z ₂₀ is determined by interpolating the gridded subsurface temperature. The other one is the Extended Reconstructed Sea Surface Temperature, version 5 (ERSSTv5; 2°×2°; Huang et al., 2017) from 1980 to 2022, which is used to construct SST_H . All monthly data are used after removing the climatological seasonal cycle and linear trends. The indexes of WWV (warm water volume above 20°C isotherm between 5°S~5°N, 120°E~80°W; gray box in Figure 1 ) and WWVw (warm water volume above 20°C isotherm between 5°S~5°N, 120°E~155°W; Red box in Figure 1 ) from 1980-2022 are directly downloaded from the website https://www.pmel.noaa.gov/elnino/upper-ocean-heat-content-and-enso. Following Bretherton et al. (1999), the effective sample size, which takes into account the serial autocorrelation at lag one, has been used in the Student’s t-test. The effective sample size (N_e ) is defined as N_e =N×(1–r_x ×r_y )/(1+r_x ×r_y ), where N is the length of time series and r_x (r_y ) is the lag-one autocorrelation coefficient of the time series for variable x (y).

To further demonstrate the role of SST_H on ENSO prediction, we also analyze the monthly data of 36 CMIP6 in historical simulations from 1900 to 2014 ( Table 1 ; Eyring et al., 2016). All monthly mean data are used after removing the climatological seasonal cycle and quadratically trends.

The role of tropical Pacific OHC in ENSO predictability is illustrated in Figure 2 . According to the numerical solution of the simplest recharge oscillator model ( Figure 2A ; Burgers et al., 2005), OHC leads SST anomalies by about 6-9 months ( Figure 2B ), especially during the early calendar months of the year ( Figure 2C ), which is consistent with McPhaden (2003) using the observation data.

In recharge oscillator model, SST_H is in phase with H, as indicated by a correlation coefficient close to 1 (red line vs. blue line in Figure 2A ). SST_H leads SST anomalies by about 6-9 months, which is also the same with H ( Figure 2B ). This suggests that SST_H can effectively represent H in this simple recharge oscillator model.

SST_H , WWV, and WWVw are predictors of SST anomalies. Firstly, WWV or WWVw leads SST anomalies by several months (blue or green line vs. black line in Figure 3A ). For example, WWV or WWVw peaks before the SST anomalies reach their maximum in 1982. This point can be further found in Figure 3B . Consistent with Izumo et al. (2019), WWV leads SST by half a year, while WWVw has a longer lead time, especially for crossing the SPB (blue and green lines in Figure 3B ). Meanwhile, the correlation coefficient between SST_H and WWV is 0.5 (significant at a 95% confidence level), which indicates that SST_H can also be a precursor for SST anomalies. The increase of SST_H in early spring suggests a recharge state in the tropical Pacific and the El Niño event may occur in the future (red and black lines in Figure 3A ). SST_H always leads Niño3.4 SST by about 10 months (red line in Figure 3B ). Among the 50 ENSO events from 1950 to 2022, 27 years were identified to be caused by SST_H signals. This finding is expected as SST_H is based on the subsurface temperature information. According to the recharge oscillator theory of ENSO, the time when OHC anomalies reach their maximum precedes the development of SST anomalies. For a longer lead time, the correlation between SST_H and the Niño3.4 index is higher than that of WWV. As shown in Figure 4 , SST_H from January to March is a predictor for the following winter ENSO SST anomalies. Particularly, T = a _* SST_H +b, where T indicates winter SST anomalies and SST_H indicates the value of SST_H from January to March. Here, a=0.0402 and b=0.017 are trained by using observational data. The prediction skill of anomaly correlation coefficient (ACC) is 0.46, which is significant at a 99% confidence level. By using this index, the spring predictability barrier is weakened.

The role of SST_H in ENSO prediction is also identified in the CMIP6 datasets. In 36 CMIP6 models, SST_H consistently leads Niño3.4 SST by about 6-16 months (gray lines in Figure 3C ). The highest correlation is about 0.7. The multi-model ensemble (MME) mean of the 36 CMIP6 models shows that SST_H leads Niño3.4 SST by 10 months with the highest correlation (0.34; purple line in Figure 3C ). This implies that SST_H can also act as a precursor for ENSO events in CMIP6 models, consistent with observation.

SST_H is effective in crossing ENSO SPB. The primary characteristic of the ENSO SPB is the occurrence of a band with the maximum decline of monthly autocorrelation in spring. This is evident from the monthly autocorrelation of SST anomaly variability and its lag gradient. It indicates that ENSO forecasting experiences the most rapid loss of predictability in spring ( Figure 5A ). Subsequently, persistent SST anomalies show higher skills for winter SST anomalies. However, SST_H still suggests a higher skill in the early spring ( Figure 5B ). When SST_H leads Niño3.4 SST anomaly about 10~16 months, their correlation is higher. Particularly, the cross-correlation between winter SST and 12 months earlier of SST_H is about 0.5, indicating that it is a superior predictor for crossing ENSO SPB compared to Niño 3.4 itself. Compared with SST_H , spring WWV exhibits a higher correlation with winter Niño3.4 ( Figure 5C ). This lead correlation coefficient is about 0.6, suggesting that WWV cannot be reflected entirely by SST_H . According to eq. (2.2), our index only reflects the low frequency of WWV. SST_H index is derived by eq. (2.2), which is the simplest relationship between T and WWV. Only by this simple relationship, we can use SST data to represent WWV. In fact, WWV is related to T and the damping of thermocline depth itself although the damping term is small (Burgers et al., 2005). As such, the WWV exhibits a higher correlation magnitude compared with SST_H . Note here although the WWV index is more effective in crossing ENSO SPB, our index can be utilized to identify the role of low frequency component of the subsurface in the tropical Pacific in ENSO prediction before 1980. Additionally, WWVw exhibits a longer lead time in ENSO prediction compared to WWV, which is consistent with Izumo et al. (2019) ( Figure 5D ).

As mentioned in the previous sections, SST_H is a predictor for ENSO events, indicating the role of subsurface temperature in ENSO development and prediction. SST_H consistently leads SST anomalies by about 10 months, allowing it to cross ENSO SPB.

The effective role of the subsurface tropical Pacific, represented by SST_H , in ENSO development and prediction before 1980 is evident in Figure 6 . SST_H has been a predictor since 1950 (red line in Figure 6A ). For example, SST_H reached its maximum before an El Niño event occurred in 1958. This point can be further observed in Figure 6B . Since 1950, SST_H always leads SST anomalies by about 10 months (the correlation coefficient is about 0.4) so that can cross the SPB.

The effectiveness of SST_H remains relatively stable on the interdecadal time scale. The correlation between the Niño3.4 index and SST_H , WWV exhibits an interdecadal change ( Figures 6C, D ). The correlation coefficient is higher when WWV leads Niño3.4 by about 4~8 months before 2000 ( Figure 6C ). However, the lead time of WWV in relation to Niño3.4 shortens to 0~4 months after 2000, indicating a weakened role of WWV in ENSO prediction, consistent with the previous study (McPhaden, 2012). For SST_H , the lead time roughly maintains at 10 months since 1985 ( Figure 6D ), suggesting that it is a relatively stable index for the lead time.

SST_H is more effective in predicting ENSO with a lead time of about 10 months compared with WWV and WWVw after the 21^st century ( Figure 7 ). We divide the time period into two periods: 1980-1999 and 2000-2022. Before 2000, SST_H leads Niño3.4 SST anomalies by about 10 months (the correlation coefficient is about 0.6; N_e =28), while after 2000, the same lead time exhibits a lower correlation (the correlation coefficient is about 0.4; N_e =31). During the 1980-1999 period, WWV led Niño3.4 SST anomalies by about 6 months (the correlation coefficient is about 0.7; N_e =41). However, after 2000, this lead time is shortened to about 3 months, and the correlation decreases (lower than SST_H for lead time > 6 months; blue lines in Figure 7 ; N_e =63). Furthermore, the cross-correlation for WWVw decreases dramatically from 0.5 in 1980-1999 (N_e =41) to 0.15 in 2000-2022 (N_e =43), indicating reduced effectiveness of WWVw after 2000 (green lines in Figure 7 ). On the other hand, SST_H suggests a relatively smaller modulation compared to WWV and WWVw for crossing the SPB (red line in Figure 7 ).

SST_H is also more effective in crossing ENSO SPB after 2000. Before 2000, WWV was a predictor in the spring ( Figure 8C ). Although the in-phase relationship between WWV and SST is small (close to 0), spring WWV anomalies lead winter SST anomalies by 8-12 months (the correlation coefficient is about 0.7), indicating that it is a superior predictor for crossing ENSO SPB compared to Niño 3.4 itself ( Figure 8A vs Figure 8C ). Following spring, persistent SST anomalies exhibit higher skills in predicting winter SST anomalies ( Figure 8A ). Therefore, WWV is a constraint for forecasts starting early in the calendar year. Both SST_H and WWVw suggests a similar role in ENSO prediction, displaying higher skills for winter and following spring SST anomalies ( Figures 8B, D ). In the 21^st century, the SPB becomes more pronounced ( Figure 8E ). Spring WWV shows a smaller skill in predicting winter SST anomalies, and the lead time is notably shortened ( Figure 8G ), which is consistent with the findings of McPhaden (2012). WWVw also exhibits a dramatic drop in effectiveness ( Figure 8H ). While SST_H suggests a slightly lower skill in the early spring ( Figure 8F ), and the lead time remains consistent with that before 2000. The cross correlation between winter SST and 12 months earlier of SST_H is about 0.5, indicating that it is an effective predictor for crossing the SPB.

In this section, we will explain why SST_H becomes more effective after the 21^st century. Mathematically, if we assume the time series of T in the simplest case: T = sin(ω ₀ t) , according to eq. (2.2), the time series of H can be derived as H = cos(ω ₀ t). The lead-lag correlation (r_H,T (τ), τ>0 indicates H leads T for τ months) between H and T is:

According to eq. (4.1), the maximum cross-correlation occurs when τ equals 1/4 ENSO period (ENSO period is 2 π ω 0 ). When the ENSO period shortens, the lead time becomes smaller. The lead-lag cross-correlation relationship is controlled by the ENSO period (black lines in Figure 7 ), which can be theoretically derived from the recharge oscillator (Jin et al., 2021; their eq. (3.8)). A shortened ENSO period can lead to a decreased lead-lagged correlation. Note here other factors (e.g., the damping term of the WWV) can also lead to this lead-lag correlation as this theoretical solution is derived by the simplest recharge oscillator. Therefore, due to this model may be too simple, the relationship between lead-lag correlation and ENSO cycle is consistent with the observation before 2000 but not after 2000.

The shortened ENSO period after the 21^st century may contribute to the reduced effectiveness of WWV and WWVw. According to recharge oscillator theory, under the effect of westerly wind forcing, the accumulated warm water in the western Pacific flows eastward, finally affecting SST anomalies in the eastern Pacific (e.g., an El Niño event). At the same time, the zonally integrated Sverdrup transport, caused by wind stress, reduces the thermocline depth in the western Pacific (Wyrtki, 1985; Jin, 1997). During the period of 1980-1999, the positive correlation coefficient between thermocline depth and winter SST anomalies in the eastern Pacific can reach the west of 140°E with a lead time of 12 months ( Figure 9A ), indicating that the western Pacific is still “recharging” 12 months before an El Niño event happens. During this period, the region of OHC (e.g., WWV, 120°E~80°W) has a reasonable lead time of about 9 months since the major area of the Pacific exhibits positive correlation coefficients. However, the situation is different in the 21^st century. A negative correlation occurs in the western Pacific with a lead time of about 12 months ( Figure 9B ), indicating that the western Pacific is “discharging” at 12 months before an El Niño event occurs. This means that the OHC in the western Pacific is building up for a La Niña event at that lead time. It shows an accelerated recharge/discharge rate (a shortened ENSO period) after the 21^st century. The correlation coefficient of the west (east) of the dateline is negative (positive) at about 9 months lead, indicating that the area of WWV or WWVw is not appropriate. On the other hand, SST_H does not need to consider this situation as it is constructed using SST anomalies. The Niño3.4 area is able to represent the major feature of ENSO activity. This is why SST_H is more effective than WWV after the 21^st century.

It should be noted that SST_H , WWV and WWVw are all less effective after the 21^st century (dashed vs. solid lines in Figure 7 ). According to eq. (4.1), the lead time is directly related to the ENSO period. As the ENSO period is longer during 1980-1999 (McPhaden, 2012), the OHC-based indexes are more effective, meaning that ENSO predictability is higher. In this section, we have shown that due to the accelerated recharge/discharge after the 21^st century, WWV and WWVw have become less effective. However, as SST_H is constructed by SST anomalies, it remains relatively effective compared to WWV and WWVw after the 21^st century.