A Hybrid ENSO Prediction System Based on the FIO−CPS and XGBoost Algorithm

Abstract

Accurate prediction of the El Niño–Southern Oscillation (ENSO) is crucial for climate change research and disaster prevention and mitigation. In recent decades, the prediction skill for ENSO has improved significantly; however, accurate forecasting at a lead time of more than six months remains challenging. By using a machine learning method called eXtreme Gradient Boosting (XGBoost), we corrected the ENSO predicted results from the First Institute of Oceanography Climate Prediction System version 2.0 (FIO−CPS v2.0) based on the satellite remote sensing sea surface temperature data, and then developed a dynamic and statistical hybrid prediction model, named FIO−CPS−HY. The latest 15 years (2007–2021) of independent testing results showed that the average anomaly correlation coefficient (ACC) and root mean square error (RMSE) of the Niño3.4 index from FIO−CPS v2.0 to FIO−CPS−HY for 7− to 13−month lead times could be increased by 57.80% (from 0.40 to 0.63) and reduced by 24.79% (from 0.86 °C to 0.65 °C), respectively. The real−time predictions from FIO−CPS−HY indicated that the sea surface state of the Niño3.4 area would likely be in neutral conditions in 2023. Although FIO−CPS−HY still has some biases in real−time prediction, this study provides possible ideas and methods to enhance short−term climate prediction ability and shows the potential of integration between machine learning and numerical models in climate research and applications.

1. Introduction

The El Niño–Southern Oscillation (ENSO) is the strongest signal of interannual variability in the coupled ocean–atmosphere system [1,2]; it not only significantly affects tropical climate and weather, but also influences extratropical areas through processes such as atmospheric teleconnection and ocean circulation [1,3,4,5,6]. These processes subsequently affect global terrestrial and marine ecosystems [7,8,9], as well as agricultural production and food security [10,11]. Therefore, accurate ENSO prediction is not only important for short−term climate prediction studies, but also crucial for disaster prevention and mitigation [12,13].

ENSO prediction has been extensively examined over the past four decades [6,13,14,15]. Traditionally, the tools for ENSO prediction can be divided into statistical and dynamical models. Statistical models are mainly based on using historical data to find the relationship between ENSO and predictors. These models include multiple linear regression, principal oscillation pattern [16], analog analysis [17], canonical correlation analysis [18,19], climatology and persistence modeling [20], singular spectrum analysis [21], Markov modeling [22], and linear inverse modeling [23]. Current statistical models can achieve an anomaly correlation coefficient (ACC) greater than 0.7 for the Niño3.4 index at a 6−month lead time [24]. These statistical models have played a steady role in ENSO predictions with the advantages of lower computational cost and ease of operation [25].

Dynamical prediction models are based on numerical models that describe changes in the ocean–atmosphere and their interactions to achieve ENSO predictions. Since the first air–sea coupled ENSO prediction model was established [26,27], these models have been improved significantly. There are three main types of dynamical prediction models, namely simple air–sea coupled models [28], intermediate coupled models [29], and coupled general circulation models or climate models [15,30]. With increasing knowledge of physical processes, development of assimilation prediction technology, and increase in available satellite observation data, the prediction skill of dynamical models for ENSO has realized significant progress [31,32], and the dynamical model has currently become the most important ENSO prediction tool [13]. Currently, there are numerous institutes that perform real−time ENSO predictions, such as the International Research Institute for Climate and Society (IRI) at Columbia University, the European Centre for Medium−Range Weather Forecasts (ECMWF), and the Beijing Climate Center of China. The IRI released the multi−model mean predictions’ ACC of the monthly Niño3.4 index (https://iri.columbia.edu/our−expertise/climate/forecasts/enso/current/#, accessed on 20 March 2023.), which can be greater than 0.8 and near 0.6 at 6− and 12−month lead times, respectively [33,34]. The skill of the latest generation prediction system, Seasonal Forecasting System 5 (SEAS5), developed by ECMWF, can exceed 0.85 at a 6−month lead time [35]. The China multi−model ensemble prediction system (CMME) provides monthly releases of ENSO real−time predictions (http://ncclcs.ncc−cma.net/Website/index.php?ChannelID=254&WCHID=252, accessed on 20 March 2023), which shows that the multi−model mean prediction skill is over 0.8 at a 6−month lead time [36].

Although current models have achieved significant progress in ENSO prediction, there are still many challenges. For example, the prediction skill at one−year lead times is still relatively low, and the ACC values of many models are still below 0.6 [33,34,37]. This indicates that the model cannot provide effective predictions. The spread of the prediction results from multiple models is large, which yields significant prediction uncertainty [6,34,38]. Moreover, both dynamical and statistical models generally suffer from the spring predictability barrier (SPB) in ENSO prediction [15,39,40,41,42]. In addition, some researchers have found that the prediction skill after 2000 is much lower than during 1982−2000 by using multi−model ensemble hindcasts [43,44,45]. All of these challenges severely limit a longer and more accurate ENSO prediction.

In recent years, a new type of statistical method represented by machine learning, also including deep learning, has shown significant potential in geoscience applications [46,47,48,49]. Relevant research on ENSO predictions has been performed via two main methods: data−driven direct prediction and model result correction [37,50,51,52,53]. In terms of direct prediction, Ham et al. [37] combined convolutional neural networks (CNNs) and transfer learning to predict ENSO at long lead times. Based on 100−year−long SST observations and CMIP5 climate model data, they achieved the first breakthrough in an ENSO prediction skill over 0.5 at a 17−month lead time. Zhou and Zhang [54] combined the principal oscillation pattern (POP) method with a CNN−LSTM network to develop a POP−Net ENSO prediction model with a prediction skill over 0.5 at a 17−month lead time. Patil et al. [55] predicted ENSO using CNNs based on more than 100 years of observed SST anomalies and vertically averaged SST anomalies. Their results show a skillful prediction for long lead times (18− to 23−month lead times) and over 0.62 in January at an 18−month lead time. In terms of model result correction, most studies have focused on multi−model ensemble results. For example, Zhang et al. [50] used the Bayesian model averaging method to perform corrections on multiple statistical and dynamical model results. They showed that the corrected prediction skill was higher than any single model value. Li [56] utilized multiple decision tree algorithms to correct the multi−model prediction results; the prediction skill was equal to or higher than the results of the dynamical and statistical models. Furthermore, Li concluded that selecting dynamical and statistical model results as input features is more important for short and long lead times, respectively. All of these studies indicate the powerful role and fairly good application prospects of machine learning in ENSO prediction. However, these studies are based on multi−model ensemble results, and correcting short−term climate predictions based on single model predictions remains challenging due to insufficient data.

Recently, the First Institute of Oceanography, Ministry of Natural Resources developed the second−generation short−term climate prediction system FIO−CPS v2.0 (First Institute of Oceanography Climate Prediction System version 2.0), which is based on the developed coupled wave Earth system model, FIO−ESM v2.0 [30,57]. Although the ENSO hindcast skill of FIO−CPS v2.0 exceeds 0.75 at a 6−month lead time, the value of the ACC decreases to 0.45 at a 12−month lead time. Furthermore, as previously mentioned, the prediction skill of FIO−CPS v2.0 during 2002−2021 is significantly lower than that during 1982−2001. Meanwhile, the SPB phenomenon still exists in ENSO predictions. To improve the ENSO prediction skills of FIO−CPS v2.0, we attempt to utilize machine learning to correct the hindcast results based on FIO−CPS v2.0 dynamical predictions, design a hybrid prediction model with FIO−CPS v2.0 and machine learning, referred to as FIO−CPS−HY, and then verify its correction skills.

The remainder of this paper is organized as follows. Section 2 introduces the method, data, and idea of correcting the results of FIO−CPS v2.0. Section 3 evaluates the corrected results based on FIO−CPS−HY. Section 4 discusses FIO−CPS−HY input features and the impact of observation data on the prediction results. Finally, Section 5 summarizes the results and provides real−time ENSO prediction.

2. Datasets and Methods

2.1. FIO−CPS v2.0 Bias Correction Ideas

The predicted bias of FIO−CPS v2.0 has two main sources: initial bias due to data assimilation [58,59] and model bias due to various assumptions and approximations introduced by the imperfect description of dynamical processes and parameterization schemes [60,61], also known as intrinsic bias. We considered both the initial and intrinsic biases and corrected predicted SST anomalies over the Niño3.4 area using the machine learning model.

Figure 1 shows a schematic diagram of the ENSO hybrid prediction model, FIO−CPS−HY (First Institute of Oceanography Climate Prediction System Hybrid), based on FIO−CPS v2.0 and machine learning. FIO−CPS−HY first characterizes the intrinsic bias (CO₂ concentration introduced to highlight the effect of external forcing) and initial bias of FIO−CPS v2.0. It obtains a time series of the intrinsic bias anomaly, the initial bias anomaly of FIO−CPS v2.0, and the global mean CO₂ concentration anomaly. Then, the two sources of biases are combined with the prediction dimension of FIO−CPS v2.0 and adapted to 1− to 13−month lead times using transfer learning [62] to correct the original hindcast SST anomaly. Finally, the corrected SST anomaly and Niño3.4 index can be obtained. Specifically, for the intrinsic bias of FIO−CPS v2.0, the original biases of the SST in the FIO−ESM v2.0 historical simulation and future scenario experiment, i.e., SSP2−4.5, are used to characterize the intrinsic bias of the model (in FIO−CPS v2.0, external forcing data, including solar radiation, CO₂ concentration, and other forcing data after 2015, were employed for the future scenario projection data of SSP2−4.5). The XGBoost (eXtreme Gradient Boosting) model in machine learning is used to predict and obtain the final intrinsic bias. We then adapt the bias values required for the target months of the FIO−CPS v2.0 hindcast and real−time predictions. The global mean CO₂ concentration required by FIO−ESM v2.0 is directly adapted to the target month of FIO−CPS v2.0. For the initial bias, the SST of the assimilation outputs is compared with the observations to characterize the initial bias. As the prediction system corresponds to the same initial field before predicting from the same month, the initial bias anomaly is relatively small and constant. Thus, the assimilated bias anomaly is adapted as the initial bias of the 13 target months starting from the same month. Finally, based on the two adapted sources of biases and the original SST anomaly predicted by FIO−CPS v2.0, a hybrid prediction model is developed to improve ENSO predictions.

2.2. Datasets

The data used in this study included observation data, FIO−ESM v2.0 simulations and projections, CO₂ concentrations, and assimilation and prediction data from FIO−CPS v2.0 (

Table 1. Data used in the construction of the FIO−CPS−HY hybrid prediction model.

Dataset		Horizontal Resolution	Temporal Resolution	Time Span
SST	ERSST v5	2° × 2°	Monthly	01/1854–01/2023
	OISST v2.1	0.25° × 0.25°	Daily	12/01/1981–12/31/2022
	FIO−ESM v2.0 Historical simulation	1.1° × (0.27–0.54°)	Monthly	01/1854–12/2014
	SSP2−4.5			01/2015–01/2023
	Assimilation			12/01/1981–11/30/2021
	FIO−CPS v2.0			01/1982–12/2022
Global Mean CO2 Concentration		——		01/1854–01/2023

). The observation data include (1) the fifth−generation Extended Reconstructed Sea Surface Temperature (ERSST v5) monthly results provided by the U.S. Oceanic and Atmospheric Administration (NOAA), with a horizontal resolution of 2° × 2° and monthly temporal resolution, spanning from 01/1854 to 01/2023 [63]. As the observational period is relatively long, we used it for intrinsic bias prediction in the prediction system. The observation data also include (2) daily results from the Optimum Interpolation Sea Surface Temperature (OISST v2.1) provided by NOAA, which incorporates satellites and other different platforms (ships, buoys, and Argo floats), with a horizontal resolution of 0.25° × 0.25° and a daily temporal resolution, spanning from 12/01/1981 to 12/31/2022 [64]. We first converted the OISST data into monthly average data. They then could be used for the initial bias calculation and transfer correction in FIO−CPS v2.0.

The SSTs from FIO−ESM v2.0 in historical simulation and SSP2−4.5 projection experiments were used. We selected the first set of experiments, i.e., r1i1p1f1, which was used to develop the short−term climate prediction system, FIO−CPS v2.0, in combination with assimilation prediction technology, with a horizontal resolution of 1.1° × (0.27–0.54°) and a monthly temporal resolution. We chose them during 01/1854−01/2023, with the historical simulated SST from 01/1854 to 12/2014 and SSP2−4.5 projected SST from 01/2015 to 01/2023. Additionally, the global mean CO₂ concentration data were provided by CMIP6 (Coupled Model Intercomparison Project phase 6).

The assimilated results of the SST, also known as the initial field of the FIO−CPS v2.0, have the same horizontal and temporal resolution for the SSTs in FIO−ESM v2.0; the time span was 1 month in advance of the starting time of FIO−CPS v2.0 for predictions, from 12/1981 to 11/2021, including 480 months. The prediction data were the historical hindcast Niño3.4 SST from FIO−CPS v2.0 with 10 ensemble members. Here, we calculated the ensemble mean results as the prediction values. The starting prediction period of the hindcast experiment was from 01/1982 to 12/2021. We selected the prediction results for the subsequent 13 lead months for the starting prediction for each, such that the sample size was 480 for each lead month.

We adopted the Niño3.4 index to characterize the ENSO. The SST anomaly, two sources of biases’ anomaly, and global mean CO₂ concentration anomaly were calculated with respect to their respective climatology in 1991–2020. We obtained the Niño3.4 index from the three−month running mean value of the SST anomaly [65].

2.3. Machine Learning Model

XGBoost (Extreme Gradient Boosting) is an integrated machine learning algorithm based on a gradient boosting regression tree (GBDT), which integrates many weaker learners to form a strong learner [66].

XGBoost is an additive model consisting of $k$ base models, assuming that the tree model for sample i, $f_t\left(x_i\right)$, to be trained in the t−th iteration is written as follows:

$${\widehat{y}}_i^{\left(t\right)}={\sum }_{k=1}^t{f}_k\left(x_i\right)={\widehat{y}}_i^{\left(t-1\right)}+f_t\left(x_i\right),$$

(1)

where${\widehat{y}}_i^{\left(t\right)}$denotes the prediction result of sample i after the t−th iteration, ${\widehat{y}}_i^{\left(t-1\right)}$denotes the prediction result of the previous t – 1 trees, and $f_t\left(x_i\right)$ denotes the model of sample i in the t−th tree. The objective function is defined as follows:

$$Ob{j}^{\left(t\right)}={\sum }_{i=1}^{n}l\left(y_i,{\widehat{y}}_i^{}\right)+{\sum }_{i=1}^t\Omega \left(f_t\right),$$

(2)

where${\sum }_{i=1}^{n}l\left(y_i,{\widehat{y}}_i^{}\right)$ is the training loss for all samples, and${\sum }_{i=1}^t\Omega \left(f_t\right)$ is the complexity summation of all t trees, added to the objective function as a regularization term to prevent overfitting of the model. The objective function is expanded by a second−order Taylor series at $f_t=0$ and simplified to obtain the final objective function as follows:

$$Ob{j}^{\left(t\right)}\cong {\sum }_{i=1}^n\left[g_i{f}_t\left(x_i\right)+\frac{1}{2}{h}_i{f}_t^2\left(x_i\right)\right]+\Omega \left(f_t\right),$$

(3)

where$g_i$ and $h_i$ are the first− and second−order derivatives of the loss function, respectively. Only $g_i$ and $h_i$ at each step are required; then, the objective function is optimized to obtain $f\left(x\right)$at each step. Finally, an overall model is obtained according to the additive model.

XGBoost has achieved strong power in major machine learning competitions [67] and has emerged in several applications in geosciences [68,69,70], especially climate prediction. For example, Zhang et al. developed a new model to predict droughts using a past drought index, meteorological measures, and climate signals from 32 stations from 1961 to 2016 in Shaanxi province, China. They selected three machine learning models including XGBoost and compared their validations for predicting the Standardized Precipitation Evapotranspiration Index (SPEI) 1–6 months in advance, and the XGBoost model more accurately predicted SPEI with a lead time of 1–6 months than other models [71]. Qian et al. developed three models including XGBoost to perform seasonal forecasts of the surface air temperature in winter in North America, and the XGBoost model hindcasts showed better forecast performances [72]. Huang et al. constructed 12 machine learning models to predict and compare daily and monthly values of solar radiation, and a stacking model using the best of these algorithms was developed to predict solar radiation. They concluded that the stacking model and the XGBoost model are the best models for predicting solar radiation [73]. Yang et al. applied many machine learning models to predict wave energy flux in the 31 geographical coordinate points in the Bohai sea separately, based on the data of the mean wave period, the significant wave height, and water depth. Their results showed that the gradient boosting machine, XGBoost, and stacked ensemble models yielded the best performance out of all algorithms [74].

Compared to the deep learning model (e.g. CNN), XGBoost can better process table data, and has stronger interpretability. In addition, it has merits such as easy parameter tuning, input data invariance, reliability, flexibility, and accuracy in many cases. And the actual performance of XGBoost is top notch in most regression and classification issues. Therefore, considering the sample size of the short−term climate prediction system is only 40 years in this study, a machine learning model that requires small training datasets is preferred over complex deep learning. Therefore, we selected XGBoost to conduct ENSO correction experiments.

2.4. Evaluation Metrics

Two statistical metrics, RMSE and ACC, were used to evaluate the mean bias and linear correlation between the prediction and observation, expressed as follows:

$$RMSE=\sqrt{\frac{{\sum }_{t=1}^N{\left(P_t-O_t\right)}^2}{N}} \text{and}$$

(4)

$$ACC=\frac{{\sum }_{t=1}^N\left(P_t-\overline{P}\right)\left(O_t-\overline{O}\right)}{\sqrt{{\sum }_{t=1}^N{\left(P_t-\overline{P}\right)}^2}\sqrt{{\sum }_{t=1}^N{\left(O_t-\overline{O}\right)}^2}},$$

(5)

where $N$ represents the sample size; $P_t$ and $O_t$ represent the Niño3.4 index predicted and observed in the t−th month, respectively; and $\overline{P}$ and $\overline{O}$ represent the mean of the prediction and observation in $N$ months, respectively.

2.5. FIO−CPS−HY’s Establishment

2.5.1. Intrinsic Bias Prediction

The XGBoost model was used to train and predict the intrinsic bias of the prediction system based on the outputs of historical simulations and future scenario experimental results from FIO−ESM v2.0. According to the scatter diagram of the SST and SST bias from the FIO−ESM v2.0 simulations and projections (Figure 2), along with the increase in the FIO−ESM v2.0 SST, the SST bias also gradually changed from negative to positive, which initially reflected a basic linear relationship between the two. Therefore, based on the FIO−ESM v2.0 SST as input, the global mean CO₂ concentration, important external forcing data affecting SST variation, was additionally introduced as a feature to achieve a more accurate regression prediction of the intrinsic bias.

In the XGBoost model, the real intrinsic bias (i.e., the bias of FIO−ESM v2.0) was regarded as ground truth. The input was the FIO−ESM v2.0 SST and global mean CO₂ concentration during the same month. The output was the predicted intrinsic bias of the given month. The dataset was randomly divided into a training set (including 1669 months) and a test set (including 360 months). In the training process in XGBoost, the GridSearch method based on 5−fold cross−validation was for tuning the model parameters. We mainly tuned nine parameters of XGBoost. Among them, n_estimators is the number of weak evaluators, and the larger the value of n_estimators, the stronger the learning ability of the model, but the easier the overfitting of the model; learn_rate is the step size of training the model, and when iterating the decision tree, it controls the iteration rate to prevent overfitting; max_depth is the maximum depth of the tree, and increasing its value will make the model stronger but more likely to overfit; min_child_weight is the minimum sample weight required on a leaf node, and when the value is large, it can prevent the model from learning local special samples to prevent the model from overfitting; gamma is the penalty term of model complexity, and the larger the value, the more conservative the algorithm will be; subsample and colsample_bytree can control the proportion of samples taken from the training data; and reg_alpha and reg_lambda are regularization penalty coefficients and can reduce model complexity.

The final model parameters are as follows: n_estimators = 500, learn_rate = 0.1, max_depth = 24, min_child_weight = 5, gamma = 0.1, subsample = 0.6, colsample_bytree = 0.6, reg_alpha = 0.01, and reg_lambda = 2.

As listed in

Table 2. RMSE and CC of the bias predicted by XGBoost and the real intrinsic bias.

Dataset	Statistical Metrics
	RMSE	CC
Training Set	0.1505	0.9937
Test Set	0.3863	0.9605
Hindcast Period	0.2010	0.9900

, the RMSE values of the predicted and real intrinsic bias for the training set, test set, and hindcast results (01/1982 to 12/2022) were all less than 0.39 °C, and correlation coefficient (CC) values were all above 0.96. Specifically, the RMSE of the hindcast results for subsequent FIO−CPS v2.0 SST correction was less than 0.21 °C, and the CC value exceeded 0.99. The time series of the real and predicted intrinsic bias for the hindcast period also indicate that the results of XGBoost fit the real intrinsic bias well (Figure 3), such that the predicted results were used as the intrinsic bias component of the subsequent transfer correction.

2.5.2. Initial Bias

The time series of the initial bias from FIO−CPS v2.0 was calculated as the difference between the assimilation SST and OISST, as shown in Figure 4. These values were all between −1 °C and 1 °C after assimilation, substantially reduced compared with the intrinsic bias (Figure 3). However, the initial bias was still not negligible. It was used as the initial bias of the subsequent transfer correction.

2.5.3. Transfer Correction of FIO−CPS v2.0

Based on the above analysis, we obtained the predicted intrinsic bias of FIO−CPS v2.0 based on the simulated and projected SST bias from FIO−ESM v2.0. In essence, the similarity between FIO−ESM v2.0 and FIO−CPS v2.0 is similar to the sample−based transfer in transfer learning. Therefore, the initial bias anomaly of FIO−CPS v2.0, predicted intrinsic bias anomaly, and global mean CO₂ concentration anomaly values could be adapted to the predicted results of the prediction system. For the initial bias anomaly of the starting prediction month, it was adapted to 13 target months starting from that month. For example, the initial bias anomaly of 12/1981 was adapted to the results of the subsequent 13 months starting from 01/1982. In addition, the intrinsic bias anomaly and global mean CO₂ concentration anomaly were directly adapted to the target month at the same time by the prediction system; e.g., the intrinsic bias anomaly and global mean CO₂ concentration anomaly values in 12/1982 were directly adapted to the results of 12/1982 predicted by the prediction system for 13 different lead times.

We constructed XGBoost models of the transfer correction for 1− to 13−month lead times in order: XGBoost_1, XGBoost_2, ..., XGBoost_13. For each correction model, the ground truth was obtained from OISST. The input had four feature factors, including the original SST anomaly in the Niño3.4 region predicted by FIO−CPS v2.0 in the same month, intrinsic bias anomaly of the prediction system, anomaly of the global mean CO₂ concentration, and initial bias anomaly of the prediction system. The output was the corrected prediction SST anomaly of each month. We divided the hindcast dataset in chronological order, with the training set’s prediction from 1982 to 2006 and the test set’s prediction from 2007 to 2021, such that we could independently evaluate the prediction skills of models in the later period. The parameter−tuning process also used the GridSearch method based on 5−fold cross−validation to optimize the parameters of each corrected model. Finally, 13 model parameter schemes were obtained.

The Niño3.4 index was calculated after applying the trained model to the test set results. A hybrid prediction model based on the combination of FIO−CPS v2.0 and XGBoost was established, denoted as FIO−CPS−HY.

3. Results

The prediction metrics (including RMSE and ACC) of FIO−CPS v2.0 and FIO−CPS−HY for the Niño3.4 index in the training set, test set, and entire hindcast period were evaluated (Figure 5 and Figure 6). The RMSE of FIO−CPS−HY in the training set had decreased to varying degrees for 13 lead times and obtained higher improvement for 7− to 13−month lead times. At a 13−month lead time, the RMSE of FIO−CPS−HY was 0.54 °C, which decreased by 26.47% compared with the values of 0.73 °C in FIO−CPS v2.0. The average RMSE for 7− to 13−month lead times from FIO−CPS v2.0 to FIO−CPS−HY could be reduced by 38.81% (from 0.63 °C to 0.38 °C) (see Figure 5a). Meanwhile, the ACC was more than 0.78 for all lead times. The improvement rate basically increased with the lead times; at 13−month lead times, the ACC of FIO−CPS−HY was 0.79, which increased by 41.76% compared with the values of 0.56 in FIO−CPS v2.0. An average ACC from FIO−CPS v2.0 to FIO−CPS−HY for 7− to 13−month lead times could be increased by 29.39% (from 0.69 to 0.90) (see Figure 5c). This shows that the hybrid prediction models were effectively trained.

As the independent period, the prediction skill of the test set is crucial to evaluate the generalization ability of our hybrid model. Except for the slightly worse than FIO−CPS v2.0 values (lower than 3.70%) at 2− to 3−month lead times, the RMSE from FIO−CPS v2.0 to FIO−CPS−HY has decreased. At a 13−month lead time, the RMSE decreased by 32.19% (from 0.95 °C to 0.64 °C), and the average RMSE for 7− to 13−month lead times could be reduced by 24.79% (from 0.86 °C to 0.65 °C) (see Figure 5b). For ACC results of FIO−CPS−HY, the improvement rate significantly increased with the lead times; an average ACC from FIO−CPS v2.0 to FIO−CPS−HY for 7− to 13−month lead times could be increased by 57.80% (from 0.40 to 0.63) (see Figure 5d). Representing an increase of more than 163.85%, the value of ACC reached more than 0.63 from 0.24 at a 13−month lead time. Furthermore, it is noteworthy that the prediction skills of FIO−CPS v2.0 in the test set were much worse than those in the training set, but the corrected results of FIO−CPS−HY have been able almost to approach and even exceed those of the earlier training period, especially for 12−to 13−month lead times.

Figure 5. RMSE and ACC of the Niño3.4 index in FIO−CPS v2.0 and FIO−CPS−HY among the OISST, training set, and test set prediction results. (a) RMSE in the training set, (b) RMSE in the test set, (c) ACC in the training set, and (d) ACC in the test set. The blue and red lines represent the RMSE and ACC of FIO−CPS v2.0 and FIO−CPS−HY for 1− to 13−month lead times, respectively. The red bar represents the percent improvement in FIO−CPS−HY relative to FIO−CPS v2.0, calculated via the absolute difference between two models divided by the RMSE or ACC of FIO−CPS v2.0 (units: %).

Figure 5. RMSE and ACC of the Niño3.4 index in FIO−CPS v2.0 and FIO−CPS−HY among the OISST, training set, and test set prediction results. (a) RMSE in the training set, (b) RMSE in the test set, (c) ACC in the training set, and (d) ACC in the test set. The blue and red lines represent the RMSE and ACC of FIO−CPS v2.0 and FIO−CPS−HY for 1− to 13−month lead times, respectively. The red bar represents the percent improvement in FIO−CPS−HY relative to FIO−CPS v2.0, calculated via the absolute difference between two models divided by the RMSE or ACC of FIO−CPS v2.0 (units: %).

Here, we could obtain the whole hindcast Niño3.4 index time series, and the difference in prediction skill between FIO−CPS v2.0 and FIO−CPS−HY is shown in Figure 6. The hindcast period (prediction from 01/1982 to 12/2021) shows that the RMSE of FIO−CPS−HY decreased to within 0.58 °C from 0.82 °C and decreased by 29.29% at the 13−month lead time, with the average RMSE for 7− to 13−month lead times could be reduced by 30.28% (from 0. 73 °C to 0.51 °C) (see Figure 6a). The ACC increased to more than 0.73 from 0.41, increasing by 79.15%; even an average ACC from FIO−CPS v2.0 to FIO−CPS−HY for 7− to 13−month lead times could be increased by 45.96% (from 0.55 to 0.81) (see Figure 6b). So, the hindcast results of FIO−CPS−HY could provide higher−quality data for the subsequent ENSO prediction research.

Figure 6. RMSE and ACC of the Niño3.4 index in FIO−CPS v2.0 and FIO−CPS−HY between the OISST and the whole hindcast results (prediction from 01/1982 to 12/2021). (a) RMSE and (b) ACC. The blue and red lines represent the results of FIO−CPS v2.0 and FIO−CPS−HY, respectively. The red bar represents the percent improvement by FIO−CPS−HY relative to FIO−CPS v2.0, calculated via the absolute difference between the two models divided by the RMSE or ACC of FIO−CPS v2.0.

Figure 6. RMSE and ACC of the Niño3.4 index in FIO−CPS v2.0 and FIO−CPS−HY between the OISST and the whole hindcast results (prediction from 01/1982 to 12/2021). (a) RMSE and (b) ACC. The blue and red lines represent the results of FIO−CPS v2.0 and FIO−CPS−HY, respectively. The red bar represents the percent improvement by FIO−CPS−HY relative to FIO−CPS v2.0, calculated via the absolute difference between the two models divided by the RMSE or ACC of FIO−CPS v2.0.

Figure 7 shows the RMSE and ACC of the Niño3.4 index of FIO−CPS v2.0 and FIO−CPS−HY as a function of the predicted target month and lead time in the test set. Compared with FIO−CPS v2.0, the prediction skills of FIO−CPS−HY significantly improved, especially in boreal summer and autumn. The significant positive difference in the ACC as the target month ranged from July to November indicates that the predictability in summer and autumn has been improved in FIO−CPS−HY (see Figure 7f). Furthermore, we found that the longer the lead times, the greater the improvement rate of FIO−CPS−HY; the ACC has increased by more than 0.46 in summer and more than 0.14 for other target months from 11−month lead times. Meanwhile, it was obvious that the SPB phenomenon has been alleviated more with the extension of lead times. In addition, the pattern of the RMSE reduction was basically similar to the pattern of the ACC improvement; the RMSE had a much sharper decline with longer lead times, decreasing by 0.21 °C from 11−month lead times (see Figure 7e). In general, FIO−CPS−HY produced a better ENSO predictability, especially in boreal summer and autumn; although SPB still existed both in FIO−CPS v2.0 and FIO−CPS−HY, FIO−CPS−HY improved it clearly and had a better persistence.

To explore the improvement of FIO−CPS−HY in more detail, the scatter diagrams of the Niño3.4 index for all lead months in the test set based on the FIO−CPS v2.0 and FIO−CPS−HY outputs are shown in Figure 8. We found that the spatial contour of scatters of FIO−CPS−HY was clearly smaller than that of FIO−CPS v2.0, especially for 8−month to 13−month lead times, and the scatterplot’s center was closer to the diagonal. Although we also noted that the corrected results of FIO−CPS−HY had a relatively poor effect on values of FIO−CPS v2.0 that are too low relative to OISST, most of the corrections were closer to the observation for longer lead times.

4. Discussion

Although the FIO−CPS−HY hybrid prediction model improved the prediction ability for ENSO, three issues require further discussion: (1) What are the effects of different input features in the correction model on the correction results? In other words, this is the issue of feature importance. (2) OISST was used for the initial bias and correction model of the prediction system while ERSST was used for the intrinsic bias prediction model. To what extent do the observation data, with different types and lengths, affect the correction results? In other words, this is the issue of the selection of training observation data. (3) The prediction skill of FIO−CPS−HY significantly increases for 7− to 13−month lead times. However, the increase in the prediction skill of FIO−CPS−HY is relatively smaller for 1− to 6−month lead times. The reasons need to be discussed. Based on these three issues, we provide further discussion.

4.1. Feature Importance

Feature importance is a measure of the magnitude of each feature input in a decision tree model. The greater the performance of a feature for decision tree split point improvement, the greater the importance of that feature. The importance values of the output feature are shown in Figure 9. We found that the hindcast SST anomaly of FIO−CPS v2.0 was the most important for 1− to 9−month lead times, with values all above 0.4, and the values also exceeded 0.3 from 10−month lead time. This is because our correction task was based on the hindcast SST anomaly results of FIO−CPS v2.0, followed by the introduction of other features. With the extension of the lead times, the importance of the hindcast SST anomaly gradually decreased; the sum of three other statistical features gradually increased (the sum of the values of the four important features is 1), especially for the intrinsic bias anomaly, which gradually increased with the lead time and was more than 0.4 from 10−month lead time. As proposed by Li [56], other statistical factors with a stronger correction role should be considered at long lead times. Meanwhile, the contribution of the global mean CO₂ concentration also played a slightly larger role than the initial bias for most lead months. This implies that SST anomaly correction with an increase in the lead time requires intrinsic bias and external forcings. From here, we propose that improvements to the physical processes, parameterization schemes, and accuracy of the response to external forcings are necessary. The initial bias was the least important component for most lead times because the overall assimilated initial bias was substantially reduced compared to the intrinsic bias. It played a relatively small role in the correction models but still had a certain contribution.

4.2. Comparison of Observation Data

In this study, ERSST v5 during 01/1854–01/2023 was used to produce the intrinsic bias prediction and introduce this intrinsic bias into the final correction model. OISST v2.1 during 12/1981–12/2022 was used to calculate the initial bias and evaluate the prediction system. To show the impact of two different observations on the final correction results, three additional sets of intrinsic bias prediction experiments using different observations and four sets of correction experiments using different observations were designed for comparison with the existing experiment above. As listed in

Table 3. Hybrid prediction experiments with different observations.

Observations for Initial Bias	Observations for Intrinsic Bias	Observations for Transfer Correction
Init_OI_198112	Intri_OI_198201	TC_OI_198201
	Intri_EROI_185401
	Intri_ER_185401
Init_ER_198112	Intri_ER_198201	TC_ER_198201
	Intri_ER_185401

For observations for initial bias, Init_OI_198112 indicates that OISST was used as the observation during 12/1981–11/2021; Init_ER_198112 indicates that ERSST was used as the observation during 01/1982–11/2021. For observations for intrinsic bias, Intri_OI_198201 indicates that OISST was used as the observation during 01/1982–01/2023; Intri_ER_198201 indicates that ERSST was used as the observation during 01/1982–01/2023; Intri_ER_185401 indicates that ERSST was used as the observation during 01/1854–01/2023; Intri_EROI_185401 indicates that ERSST was used as the observation for 01/1854–12/1981 and OISST was used as the observation for 01/1982–01/2023. For observations for transfer correction, TC_OI_198201 indicates that OISST was used as the observation during 01/1982–12/2022; TC_ER_198201 indicates that ERSST was used as the observation during 01/1982–12/2022. The experiment combination, i.e., Init_OI_198112 + Intri_ER_185401 + TC_OI_198201, represents the existing experiment mentioned above in this paper.

, a total of five experiments containing the initial bias, intrinsic bias, and transfer correction were designed based on the four intrinsic bias prediction experiments.

To evaluate the performances of different experiments at long lead times, we compare the average RMSE and ACC for 7− to 13−month lead times in

Table 4. Average RMSE and ACC of the Niño3.4 index in the test set for 7− to 13−month lead times using the hybrid prediction model with different observation data.

Statistical Metrics	Observations
	Init_OI_198112			Init_ER_198112
	Intri_OI_198201	Intri_EROI_185401	Intri_ER_185401	Intri_ER_198201	Intri_ER_185401
	TC_OI_198201			TC_ER_198201
RMSE	0.6849	0.6500	0.6537	0.6961	0.6930
ACC	0.5842	0.6338	0.6297	0.6442	0.6598

. The results show that all the experiments’ RMSE and ACC changed across a small range. The maximum difference of the RMSE was within 0.05 °C, and the ACC was between 0.58 and 0.66. In addition to an experiment combination, i.e., Init_OI_198112 + Intri_OI_198201 + TC_OI_198201, with slightly poor prediction skills, other different observation data selections have a small impact on the hybrid prediction results, also reflecting the stability and robustness of our hybrid model. The ACC results of ERSST were the best observation for the initial bias, intrinsic bias, and prediction system, but the improvement was not significant compared with our existing experimental results. Meanwhile, the existing experiment has a minimum RMSE of only 0.65 °C. Considering that the OISST has been widely used in the evaluation of ENSO prediction skill [24,33,34], as well as integrated satellite data with higher spatial and temporal resolutions, we selected the existing experiment as the final result.

4.3. Prediction Skill of FIO−CPS−HY

From the results above, we can find that the prediction skill of FIO−CPS−HY significantly increases for 7− to 13−month lead times. However, the increase in the prediction skill of FIO−CPS−HY is relatively smaller for 1− to 6−month lead times. Some reasons need to be discussed.

As shown in Figure 5a,c, The RMSE and ACC in the training set for 1− to 6−month lead times from FIO−CPS v2.0 showed good prediction performance. The average RMSE and ACC for 1− to 6−month lead times from FIO−CPS v2.0 was less than 0.40 °C and more than 0.88. Therefore, it is easy to understand that the statistics’ improvement for 1− to 6−month lead times from FIO−CPS v2.0 to FIO−CPS−HY was relatively small. The prediction error increases with the increase in lead time; in order to reveal the error characteristic changes in FIO−CPS v2.0 in more detail, Figure 10 shows the scatter diagrams of the Niño3.4 index for all lead months in the training set based on the FIO−CPS v2.0 and FIO−CPS−HY outputs. It is clear that the scatter points gradually deviate from the diagonal with the increase in lead time. Meanwhile, when we conducted ENSO correction experiments for 13 lead months separately, we adopted the same parameter range as large as possible to tune parameters, so as to ensure the fairness of parameter tuning and avoid the individual influence of different parameter−tuning schemes on the results of one model. Therefore, FIO−CPS−HY has more room for correction to achieve better results, on the premise of distribution characteristics with larger deviations in longer lead months.

As mentioned above, the ENSO prediction skill after 2000 is much lower than before. The average RMSE and ACC for 1− to 6−month lead times from FIO−CPS v2.0 were less than 0.46 °C and more than 0.80 (see Figure 5b,d). Following the rules of the training set, it seems that FIO−CPS−HY should have more room for correction in the test set. However, the scatter point data distribution characteristics of the test set have changed compared with the training set (see Figure 8). The coefficient of the best−fit linear regression line of FIO−CPS v2.0 became smaller, and the overall rotation was clockwise compared with the training set, which makes it harder for FIO−CPS−HY to correct FIO−CPS v2.0 values near the diagonal, especially on values of FIO−CPS v2.0 too low for OISST. Therefore, the statistics’ improvement for 1− to 6−month lead times from FIO−CPS v2.0 to FIO−CPS−HY was smaller than that for the training set. For 7− to 13−month lead times, the scatter point data distribution characteristics of the test set changed significantly compared with the training set, and the coefficient of the best−fit linear regression line of FIO−CPS v2.0 became much smaller. However, most of the scatter points are distributed in places where the value of FIO−CPS v2.0 is too high for OISST, which gives FIO−CPS−HY more room to correct the large deviation. This is the reason why the prediction skill of FIO−CPS−HY significantly increases for 7− to 13−month lead times.

In summary, the scatterplot distribution characteristics’ changes for FIO−CPS v2.0 and OISST make it easier for FIO−CPS−HY to correct the results for 7− to 13−month lead times but more difficult for FIO−CPS−HY to correct the results for 1− to 6−month lead times. This reflects that the prediction skills and bias characteristics of FIO−CPS v2.0 have changed in the recent 20 years, which suggests that we should strive to improve the prediction skills and carry out further dynamical model bias analysis and research in the future.

5. Conclusions

Based on the latest generation of the FIO−CPS v2.0 short−term climate prediction system and the XGBoost machine learning method, the hybrid prediction model FIO−CPS−HY was developed for ENSO prediction. Considering the initial bias and intrinsic bias of ENSO predictions, FIO−CPS−HY could find the nonlinear relationship among the FIO−CPS v2.0 hindcast SST, two−part bias, and observation while correcting the prediction results of the dynamical model.

The prediction skill of the Niño3.4 index in 2007–2021 significantly improved for 7− to 13−month lead times. The average ACC and RMSE from FIO−CPS v2.0 to FIO−CPS−HY for 7− to 13−month lead times could be increased by 57.80% (from 0.40 to 0.63) and reduced by 24.79% (from 0.86 °C to 0.65 °C). At the same time, a spring predictability barrier exists in FIO−CPS v2.0 and FIO−CPS−HY, but this problem is less pronounced in FIO−CPS−HY. Moreover, the prediction skill of FIO−CPS−HY significantly increases in boreal summer and autumn from 11−month lead times.

Based on the FIO−CPS−HY, we present the corrected results of the Niño3.4 index predicted from 09/2022 to 02/2023 in Figure 11. The latest 6 months’ predictions of FIO−CPS−HY show that the Niño3.4 index has dropped to a certain extent compared with FIO−CPS v2.0 and will be most likely to increase until March–April–May 2023 and then decrease. In general, the real−time predictions from FIO−CPS−HY indicate that the sea surface state of the Niño3.4 area would likely be in neutral conditions in 2023.