El Niño–Southern Oscillation Prediction Based on the Global Atmospheric Oscillation in CMIP6 Models

Abstract

In this work, the preindustrial control (piControl) and Historical experiments results from climatic Earth system models participating in the sixth phase of the Coupled Model Intercomparison Project (CMIP6) are analyzed for their ability to predict the El Niño–Southern Oscillation (ENSO). Using the principal component method, it is shown that the Global Atmospheric Oscillation (GAO), of which the ENSO is an element, is the main mode of interannual variability of planetary anomalies of surface air temperature (SAT) and atmospheric sea level pressure (SLP) in the ensemble of 50 CMIP6 models. It turns out that the CMIP6 ensemble of models reproduces the planetary structure of the GAO and its west–east dynamics with a period of approximately 3.7 years. The models showed that the GAO combines ENSO teleconnections with the tropics of the Indian and Atlantic Oceans, and with temperate and high latitudes. To predict strong El Niño and La Niña events, we used a predictor index (PGAO) obtained earlier from observation data and reanalyses. The predictive ability of the PGAO is based on the west–east propagation of planetary structures of SAT and SLP anomalies characteristic of the GAO. Those CMIP6 models have been found that reproduce well the west–east spread of the GAO, with El Niño and La Niña being phases of this process. Thanks to this, these events can be predicted with approximately a year’s lead time, thereby overcoming the so-called spring predictability barrier (SPB) of the ENSO. Thus, the influence of global anomalies of SAT and SLP on the ENSO is shown, taking into account that it may increase the reliability of the early forecast of El Niño and La Niña events.

1. Introduction

The El Niño–Southern Oscillation (ENSO) is the strongest interannual signal in the global climate system and a leading source of seasonal climate predictions worldwide [1]. ENSO prediction relies on dynamical and statistical models to estimate the predictability and skill of ENSO predictions, which vary seasonally [2,3] and from decade to decade [4]. Interdecadal variations in ENSO predictability are largely related to slowly changing ocean conditions, with short-term atmospheric oscillations often limiting ENSO predictability on seasonal timescales. Despite the constantly improving models and their initialization schemes, the quality of ENSO prediction remains low for forecasts passing the so-called spring predictability barrier (SPB) [5,6,7,8]. Therefore, the ENSO forecast system based solely on the physical processes of the Pacific Basin limits forecasts of El Niño and La Niña events to less than one year [9].

A hierarchy of ENSO forecasting schemes is developed in [10,11] that includes statistical schemes and physical models. Statistical schemes are generally based on linear statistical methods and can be divided into models that use low-frequency oscillations of the atmosphere (sea level pressure or surface wind) or the ocean (surface temperature or upper ocean heat content) as predictors. Physical models consist of coupled atmosphere–ocean models of varying complexity, ranging from simplified shallow water models to coupled general circulation models.

With the development of coupled models, improvements in initialization schemes, and new theoretical studies, some progress has been made in ENSO forecasting on time scales from months to seasons [12]. However, there are two fundamental directions: improving practical forecasting skills and theoretical study of the inner limit of predictability. The former includes progress in coupled models, data assimilation, ensemble forecasting, etc., and the latter focuses on efforts to study the optimal error growth and estimate the inner limit of ENSO predictability. However, despite this progress and the existence of a number of models that predict ENSO dynamics, the prediction of this phenomenon is not always successful. An example is the failed El Niño forecasts of 2014–2015 [13].

Atmospheric processes important for ENSO forecasting include components of variability on seasonal timescales [1]. Anomalies of the Pacific atmospheric circulation in the winter–spring period can have a long-term impact on the tropical climate in the following summer through signal transmission via sea surface temperature (SST). In [14], an index based on sea level pressure (SLP) data from February to March in the region surrounding Hawaii is proposed, and it is shown that this region is the most informative part of a larger SLP structure that has a delayed impact on the ENSO through signal transmission via SST.

The Pacific meridional mode (PMM) is the leading mode of ocean–atmosphere interaction in the subtropical northeastern Pacific Ocean and plays an important role in the occurrence of the ENSO. One study [15] examined the performance of the Coupled Model Intercomparison Project Phase 5 and 6 (CMIP5/CMIP6) models in reproducing the PMM and its relationship with the ENSO. It was found that in the CMIP5/CMIP6, there is a large inter-model spread in the influence of the PMM spring state on the ENSO phase in the following winter. The PMM also plays an important role in modulating the development of the ENSO [16]. Some CMIP6 models can correctly reproduce the PMM characteristics, but show inconsistent PMM interactions with the ENSO due to differences in the feedbacks between wind, evaporation, and SST.

Knowledge of the ENSO drivers and underlying mechanisms is crucial to improve ENSO forecasting, which is still a research topic. The existing teleconnections between the tropical Atlantic and Pacific Oceans via atmosphere–ocean interactions during the first and last decades of the twentieth century indicate the possibility of ENSO forecasting [17]. Understanding this interdecadal modulation of the Atlantic–Pacific teleconnections may help improve forecasts of ENSO and its associated impacts.

Using a coupled ocean–atmosphere general circulation model, Luo et al. [18] showed that some past ENSO events could be predicted retrospectively with a lead time of more than 1 year. For example, the El Niño state in the winter of 1997/98 can be predicted to some extent with a lead time of about 1.5 years, but with a weak intensity and a large phase lag in predicting the onset of this exceptionally strong event. This is due to the influence of active and intense stochastic westerly wind bursts (WWB) [19] during the period from late 1996 to mid-1997, which are usually unpredictable on seasonal time scales.

Since 2002, the International Research Institute for Climate and Society (IRI), later in partnership with the Climate Prediction Center (CPC), has produced a multi-model ENSO prediction product, informally called the IRI ENSO Predictions Plume [20]. Barnston et al. [21] explore measures to improve the robustness and usability of this product, including bias and amplitude corrections, the multi-model ensemble method, the probability distribution formulation, and the output format. Barnston et al. [22] evaluate the ENSO hindcasts and real-time forecasts from the North American Multi-Model Ensemble for 1982–2015. Testing the eight individual models in this ensemble shows that some of them consistently produce better forecasts than others. In particular, the amplitudes of some ENSO model forecasts are too large when crossing the SPB.

Kirtman and Min [23] report the results of a large sample of coupled ENSO hindcasts for 1982–1998, using a forecasting system based on the National Center for Atmospheric Research (NCAR) CCSM3.0 model and an ocean data assimilation system provided by the National Oceanic and Atmospheric Administration (NOAA) Geophysical Fluid Dynamics Laboratory (GFDL). The robustness of the ENSO forecasting system is examined from both deterministic and probabilistic perspectives and compared with the current NOAA Climate Forecast System (CFS). It is then combined with the CFS to create a multi-model forecasting system. Using deterministic and probabilistic performance metrics, it is shown that the multi-model combination is generally as effective as the model showing the best ENSO forecasting ability.

Currently, ENSO can be predicted fairly well up to about 6 months in advance, but large errors and uncertainties remain in real-time forecasts [24]. At the same time, various approaches are used to improve the understanding of ENSO processes, and various models for ENSO forecasting have been developed, including linear statistical models based on the analysis of the principal oscillation patterns (POP) and convolutional neural networks (CNN). Using the artificial neural network method, Lubkov et al. [25] showed the possibility of forecasting El Niño and La Niña events, taking into account their division into two types, Central and East Pacific, with a lead time of 3 to 9 months.

In [26], a method based on a recurrent neural network called an echo state network (ESN) is developed that can be trained efficiently to predict different ENSO indices despite their relatively high noise levels. To do this, the ESN model is trained on the low-frequency variability of ENSO indices, and then potential future high-frequency variability is estimated from specific patterns in its history. This method reveals the importance of interactions across different time scales in the mechanisms underlying ENSO, and that the low-frequency component of ENSO carries substantial predictive power, which can be exploited.

In [27], a data-driven model for analyzing and forecasting spatially distributed time series is proposed. The model is based on the linear dynamical mode of decomposition of the observed data, which is derived from a recently developed nonlinear dimensionality reduction approach. The key point of this approach is its ability to account for simple dynamic properties of the observed system by identifying the dominant time scales of the system. This method is applied to the tropical SST anomaly field, where the main variability signal is the ENSO.

The ENSO response to global warming remains uncertain, challenging ENSO projections under warming conditions. Zheng et al. [28] investigate changes in ENSO characteristics and predictability in idealized simulations with fourfold greenhouse gas forcing based on seven general circulation models. ENSO variability is shown to weaken in warmer climates compared to control experiments, with the neutral state lasting longer, while ENSO active states are shorter and shift toward La Niña. However, under warming conditions, the six-month ENSO predictability decreases slightly in five models and increases in two models.

Some climate models exhibit systematic errors in the representation of the ENSO, which are significant in tropical seasonal forecasting even at short lead times [29]. ENSO-related forecast winter/spring SST anomalies are significantly broadened or westward-shifted, and are too persistent during the exit phase of El Niño and La Niña events.

Results from testing a global anomaly-linked ENSO forecasting system, including a general circulation model and an ocean data assimilation system, show that ensembled forecasts of El Niño and La Niña events are useful for up to 7–9 months lead time [30]. Thus, ENSO forecasting well in advance (more than a year) is still an important and challenging scientific task. Since a unified and complete ENSO theory has not yet been developed, related indicators such as the Niño3.4 index and the Southern Oscillation index (SOI) are used to forecast the ENSO using appropriate numerical models. However, since the ENSO is a very complex and dynamic phenomenon, and the Niño3.4 and SOI indices mix many low-frequency and high-frequency components, the prediction accuracy of modern numerical ENSO forecasting methods still remains quite low [31].

In early 2014, several forecast systems predicted a strong El Niño event similar to that of 1997/1998 for the following winter of 2014/2015 [32]. However, the result was small positive SST anomalies in the tropical Pacific. In contrast, the following winter of 2015/2016 saw one of the strongest El Niño events on record, which was not predicted in early 2015 [20]. Thus, despite the demonstrated hindcast predictability of the ENSO, real-time forecasts of El Niño and La Niña events and related phenomena outside the tropical Pacific are not always correct. There is therefore a need to examine in more detail the predictability of extratropical Pacific anomalous features known to be associated with the ENSO (e.g., the Indian monsoon and Sahel rainfall, southern African drought, and extratropical SST) [11]. In addition, the relative importance of different oceanic and atmospheric physical mechanisms for the ENSO needs to be established. The seasonal dependence of ENSO predictability is evident in many models, but the processes responsible are not fully understood, and its significance is still a matter of scientific debate. Moreover, there are marked interdecadal differences in the accuracy of ENSO forecasts, and the reasons for this are not yet fully understood.

In the works [33,34], based on observational data, reanalyses and CMIP5 models, it is shown that the ENSO is an element of a planetary phenomenon called the Global Atmospheric Oscillation (GAO). The GAO combines ENSO teleconnections with the tropics of both the Indian and Atlantic Oceans, as well as with temperate and high latitudes. Moreover, these teleconnections act in both directions: from the ENSO on the anomalies all over the Earth, and from the SAT and SLP anomalies outside the tropics of the Pacific Ocean on the ENSO. The structure of the GAO spreads from west to east; that is, it has the property of a planetary wave, making a complete revolution around the Earth in approximately 3.7 years. Decomposition of the observational data into empirical orthogonal functions (EOFs) showed that the GAO is the first principal component of the interannual variations in the global SAT and SLP anomalies [35]. In [36], based on the west–east propagation of the planetary spatial structure of the GAO, a mechanism for predicting El Niño and La Niña events is proposed. In [37], it is shown that some of the CMIP6 models reproduce the planetary structure of the GAO, and demonstrate strong ENSO teleconnections with SST and SLP anomalies outside the tropical Pacific Ocean. The aim of this work is to test whether the CMIP6 ensemble of models reproduces the west–east spread of the planetary structure of the GAO, and to assess the feasibility of successful ENSO forecasting by CMIP6 models based on the GAO.

2. Materials and Methods

The global monthly values of surface air temperature (SAT) and sea level pressure (SLP) obtained from the preindustrial control (piControl) experiment of the earth system climate models included in CMIP6 are analyzed [38]. In the piControl experiment, the forcing from any external forces is constant and is at the preindustrial period level. The only periodic external force affecting the climate system in this experiment is the annual cycle of solar heat input. The developer organizations, names, and main characteristics of the analyzed 50 CMIP6 models are presented in

Table 1. Main characteristics of the studied CMIP6 models and cross-correlation of the 1st empirical orthogonal functions of the interannual oscillations of the SAT and SLP anomalies with EONI and GAO1. All correlations presented in the table are 95% significant.

Model Developer	Model Name	Average Model Resolution (°lat × °lon)	Number of Years in the piControl Experiment	Max. Corr. EOF1 SAT and EONI	Shift at Max. Corr. EOF1 SAT and EONI (Months)	Max. Corr. EOF1 SLP and GAO1	Shift at Max. Corr. EOF1 SLP and GAO1 (Months)
AS-RCEC	TaiESM1	0.938 × 1.250	500	0.87	3	0.80	0
AWI	AWI-CM-1-1-MR	0.938 × 0.938	500	0.82	3	0.84	0
BCC	BCC-CSM2-MR	1.125 × 1.125	600	0.83	3	0.83	0
BCC	BCC-ESM1	2.812 × 2.812	451	0.81	1	0.69	0
CAMS	CAMS-CSM1-0	1.125 × 1.125	500	0.87	2	0.85	0
CAS	CAS-ESM2-0	1.406 × 1.406	550	0.93	2	0.77	0
CAS	FGOALS-f3-L	1.000 × 1.250	561	0.92	3	0.82	0
CAS	FGOALS-g3	2.250 × 2.000	700	0.89	3	0.86	0
CCCma	CanESM5	2.812 × 2.812	1000	0.83	3	0.80	0
CCCma	CanESM5-CanOE	2.812 × 2.812	501	0.84	3	0.79	0
CMCC	CMCC-CM2-SR5	0.938 × 1.250	500	0.87	3	0.83	0
CMCC	CMCC-ESM2	0.938 × 1.250	500	0.90	3	0.88	0
CNRM-CERFACS	CNRM-CM6-1	1.406 × 1.406	500	0.86	2	0.79	0
CNRM-CERFACS	CNRM-CM6-1-HR	0.500 × 0.500	300	0.83	2	0.60	0
CNRM-CERFACS	CNRM-ESM2-1	1.406 × 1.406	500	0.85	2	0.81	0
CSIRO-ARCCSS	ACCESS-CM2	1.250 × 1.875	500	0.87	1	0.78	0
CSIRO	ACCESS-ESM1-5	1.241 × 1.875	900	0.81	3	0.69	0
E3SM-Project	E3SM-1-0	1.000 × 1.000	500	0.84	2	0.82	0
EC-Earth-Consortium	EC-Earth3	0.703 × 0.703	501	0.85	2	0.78	0
EC-Earth-Consortium	EC-Earth3-AerChem	0.703 × 0.703	311	0.75	2	0.68	0
EC-Earth-Consortium	EC-Earth3-Veg	0.703 × 0.703	500	0.85	2	0.77	0
EC-Earth-Consortium	EC-Earth3-Veg-LR	1.125 × 1.125	501	0.80	2	0.70	0
FIO-QLNM	FIO-ESM-2-0	0.938 × 1.250	575	0.87	2	0.84	0
HAMMOZ-Consortium	MPI-ESM-1-2-HAM	1.875 × 1.875	1000	0.84	2	0.74	0
INM	INM-CM4-8	1.500 × 2.000	531	0.80	1	0.44	0
INM	INM-CM5-0	1.500 × 2.000	1201	0.80	1	0.60	0
IPSL	IPSL-CM6A-LR	1.259 × 2.500	800	0.87	3	0.81	0
KIOST	KIOST-ESM	1.875 × 1.875	500	0.89	3	0.77	0
MIROC	MIROC-ES2L	2.812 × 2.812	500	0.92	3	0.86	0
MIROC	MIROC6	1.406 × 1.406	800	0.92	4	0.88	0
MOHC	HadGEM3-GC31-LL	1.250 × 1.875	500	0.84	2	0.81	0
MOHC	HadGEM3-GC31-MM	0.556 × 0.833	500	0.89	2	0.84	0
MPI-M	MPI-ESM1-2-HR	0.938 × 0.938	500	0.89	1	0.80	0
MPI-M	MPI-ESM1-2-LR	1.875 × 1.875	1000	0.89	2	0.75	0
MRI	MRI-ESM2-0	1.125 × 1.125	701	0.87	2	0.78	0
NASA-GISS	GISS-E2-1-G	2.000 × 2.500	851	0.90	3	0.84	0
NASA-GISS	GISS-E2-1-H	2.000 × 2.500	401	0.89	2	0.77	0
NCAR	CESM2	0.938 × 1.250	1200	0.88	3	0.86	0
NCAR	CESM2-FV2	1.875 × 2.500	500	0.88	3	0.85	0
NCAR	CESM2-WACCM	0.938 × 1.250	499	0.87	3	0.83	0
NCAR	CESM2-WACCM-FV2	1.875 × 2.500	500	0.90	3	0.88	0
NCC	NorCPM1	1.875 × 2.500	500	0.85	3	0.83	0
NCC	NorESM2-LM	1.875 × 2.500	501	0.85	3	0.86	0
NCC	NorESM2-MM	0.938 × 1.250	500	0.85	3	0.86	0
NOAA-GFDL	GFDL-CM4	1.000 × 1.250	500	0.86	2	0.84	0
NOAA-GFDL	GFDL-ESM4	1.000 × 1.250	500	0.88	2	0.86	0
NUIST	NESM3	1.875 × 1.875	500	0.84	1	0.77	0
SNU	SAM0-UNICON	0.938 × 1.250	700	0.82	3	0.80	0
THU	CIESM	0.938 × 1.250	500	0.88	3	0.85	0
UA	MCM-UA-1-0	2.250 × 3.750	500	0.89	1	0.83	0

(columns 1–4). The duration of the piControl experiment ranged from 300 to 1200 years, depending on the model, with the most common value being 500 years (Table 1, column 4). The 50 models under study were selected based on the results of comparison of extratropical ENSO teleconnections in CMIP6 models [37]. For validation, the main results were additionally calculated separately for the first 150 years and the following 150 years of the piControl experiment (Supplementary Materials, Figures S8–S11). These results turned out to be close to the results obtained by the whole piControl experiment.

Additionally, to analyze the stability of the obtained results in the conditions of the modern changing climate, the Historical experiment data 58 CMIP6 models for 1850–2014 were studied, which includes radiative forcing from changes in the concentration of greenhouse gases in the atmosphere. The results obtained from the Historical experiment (Supplementary Materials, Tables S1–S3, Figures S1–S6) were close to the results obtained from the piControl experiment.

To identify El Niño and La Niña events, an index called the Extended Oceanic Niño Index (EONI) was used [39]. EONI represents the average values of SAT anomalies in the central-eastern equatorial Pacific Ocean (5° N–5° S, 170–80° W) (Figure 1a). The average annual cycle of SAT and SLP for the entire piControl experiment period was calculated for each of the analyzed models. Then, at each node of the model grid, SAT and SLP anomalies were calculated relative to this average annual cycle. The SAT anomalies were then averaged over the region (5° N–5° S, 170–80° W), and the resulting time series was smoothed with a 3-month moving average to obtain each model’s own EONI time series. El Niño and La Niña events were then determined separately for each model based on the resulting EONI time series. For El Niño (La Niña), EONI values should be continuously above +0.5 °C (be less than −0.5 °C) for 5 months or more. We also used the Oceanic Niño Index (ONI), calculated as average SST anomalies in the Niño3.4 region (5° N–5° S, 170–120° W), but it covers a smaller area than EONI and is therefore less representative in our opinion. However, the results obtained for these two indices turned out to be quite close.

In order to identify the current phase of the GAO, the GAO1 index was used, which is calculated as an algebraic combination of normalized values of the SLP anomalies in 10 areas that coincide with the maxima and minima in the spatial structure of the SLP anomalies of the GAO field, coinciding in phase with the ENSO (Figure 2a) [33]. GAO1 is calculated using the following formula:

GAO1 = P(5° S–5° N, 145–155° E) + P(5° S–5° N, 55–65° E) + P(5° S–5° N,
35–25° W) + P(55–65° N, 95–85° W) + P(65–55° S, 95–85° W) − P(5° S–5° N,
95–85° W) − P(45–55° N, 175–165° W) − P(45–55° N, 15–5° W) − P(55–45° S,
15–5° W) − P(55–45° S, 175–165° W),

where the Ps are the average SLP anomalies in areas with given coordinates. In essence, the GAO1 index is an extension to the entire tropical belt, as well as to the temperate and high latitudes, of the equatorial Southern Oscillation index (equatorial SOI), which is calculated as the normalized difference between the standardized SLP anomalies averaged in regions (5° N–5° S, 80–130° W) and (5° N–5° S, 90–140° E).

To predict El Niño and La Niña events, a predictor index based on the west–east distribution of the spatial structure of the GAO (Predictor GAO–PGAO) is used. Previously [36], based on cross-correlations, it was found that the PGAO is ahead of EONI and GAO1 by approximately 1 year. The PGAO is calculated as an algebraic combination of the normalized values of the average SAT (T) and SLP (P) anomalies in 15 regions of the GAO field, which is ahead of the El Niño phase:

PGAO = T(20–50° N, 160–130° W) + T(35–25° S, 160–80° W) + T(65–45° S,
150° E–160° W) + T(65–45° S, 60° W–0°) + T(40–70° N, 90° E–180°) − T(30° S–
30° N, 60° W–180°) − T(60–31° S, 90–120° E) − T(31–60° N, 100–40° W) + P(50–
70° N, 170° E–120° W) + P(70–50° S, 170° E–120° W) + P(60° S–20° N,
40–80° E) + P(30° S–30° N, 70–10° W) − P(0–40° N, 120° E–120° W) −
P(45–25° S, 120° E–70° W) − P(50–70° N, 50–90° E).

Methods for identifying the main modes of climate variability are the subject of works by various authors. In [40], by 6-yr highpass- and lowpass- filtering, the time variability of the leading EOF of the global SST field is separated into two components: one identified with the “ENSO cycle-related” variability on the interannual timescale, and the other a linearly independent “residual” comprising all the interdecadal variability in the record. In [41], it was found that the ENSO is the second most important dominant nonlinear mode for global SST anomalies after the annual cycle. We applied geographically weighted Principal Component Analysis (PCA), called the method of decomposition into empirical orthogonal functions (EOFs), to identify the main modes of variability of SAT and SLP anomalies on interannual periods of oscillations characteristic of the GAO and ENSO. For this purpose, preliminary bandpass filtering of SAT and SLP anomaly series from 2 to 7 years was performed using the Butterworth filter at each grid node. Bandpass filtering was performed for periods from 2 to 7 years, since it was found in [37] that it is in this range that the highest oscillation energies of the EONI index of the 50 CMIP6 models under study are observed. In addition, to suppress strong variability of the SAT and SLP anomalies at high latitudes compared to the tropics, the time series at each grid node were multiplied by the root of the cosine of the latitude of this node. Multiplying by the root cosine of the latitude applies geographic weighting to the data, so that grid cells with smaller areas will play a smaller role in finding the principal components. After that, PCA was applied to the obtained filtered and normalized SAT and SLP anomalies.

3. Results and Discussion

Based on the results of the piControl experiment of 50 CMIP6 models, indicated in Table 1, the ensembled field of the first principal component (PC1—spatial structure of EOF1) of the interannual variability of SAT anomalies (PC1 SAT) was constructed (Figure 1a). When constructing it for each model separately, the first 10 principal components were calculated for the SAT anomalies, normalized and filtered for oscillation periods from 2 to 7 years at each node of the model grid. The fields of the first principal component found for each model were interpolated onto a single 1° × 1° grid. Then, the interpolated fields of PC1 SAT were averaged over the 50 CMIP6 models under study to obtain the inter-model mean. The field of the corresponding inter-model spread (standard deviation) is given in the Supplementary Materials (Figure S12).

The structure of the inter-model mean PC1 SAT (Figure 1a) in its shape in many details repeats the global structure of the ensembled field of the difference in SAT anomalies between opposite ENSO phases for 50 CMIP6 models under study [37]—the planetary structure of the GAO SAT. Based on this, it can be assumed that the GAO and its ENSO element correspond to the main mode of global climate variability of SAT anomalies on interannual time scales for most of the 50 CMIP6 models under study. The relative degree of this correspondence for each model separately can be estimated from the maximum values of cross-correlations between EOF1 SAT and the EONI (Table 1, column 5), and the shifts at which these maximum cross-correlations are achieved (Table 1, column 6).

The assumption that the GAO is the main mode of global climate variability of the studied CMIP6 models on interannual time scales is confirmed by the ensembled field of the first principal component of interannual variability of the SLP anomalies (PC1 SLP), calculated in a way similar to PC1 SAT based on the results of the piControl experiment of 50 CMIP6 models (Figure 2a). The structure of the inter-model mean PC1 SLP also almost completely corresponds to the planetary structure of the ensembled field of the difference in SLP anomalies between opposite ENSO phases for 50 studied CMIP6 models [37]—the planetary structure of the GAO SLP. The relative degree of this correspondence for each of the 50 models studied separately can be estimated from the maximum values of cross-correlations between EOF1 SLP and GAO1 (Table 1, column 7), and the shifts at which these maximum cross-correlations are achieved (Table 1, column 8).

Due to preliminary bandpass filtering from 2 to 7 years, almost all fluctuations of SAT and SLP anomalies not related to the interannual variability were suppressed. However, the use of bandpass filtering significantly distorted the original SAT and SLP anomalies, and, accordingly, the PCA results. The obtained principal components relate only to the interannual variability of the SAT and SLP anomalies, and not to the SAT and SLP data as a whole. Nevertheless, without applying the said bandpass filtering, it would not have been possible to ensure that for each of the 50 studied models, the first principal component of the SAT and SLP anomalies corresponded to the strongest signal of interannual climatic variability, which is the ENSO. Therefore, it is partly due to the applied filtering that the first principal components of the SAT and SLP anomalies of 50 CMIP6 models under study correspond to the ENSO and GAO. And when they are averaged, the inter-model mean fields PC1 SAT and PC1 SLP are obtained (Figure 1a and Figure 2a), repeating in many details the planetary structures of the GAO obtained in [37] using the compositional method. Also, due to the features of the algorithm for decomposition into principal components, they were taken with the opposite sign if necessary. This was required so that when averaging, the PCs were of the same phase.

Table 1 contains the maximum values of cross-correlations of EOF1 of the global interannual variability of SAT anomalies (EOF1 SAT) and the EONI index (column 5) for shifts from −60 to +60 months, as well as the shifts to which these values correspond (column 6). All of these correlations turned out to be statistically significant at the 95% level. High values of cross-correlations of EOF1 SAT and the EONI (the average correlation is 0.86) as well as the El Niño’s SST anomalies spatial structure in the tropical Pacific Ocean (Figure 1a) confirm that the first principal components of the global interannual variability of SAT anomalies of the studied CMIP6 models relate precisely to the ENSO. At the same time, small lags (1–3 months) of EOF1 SAT with respect to the EONI are observed, which indicates that the planetary spatial structure of PC1 SAT (Figure 1a) is associated with the response to El Niño and La Niña events caused by the global ENSO teleconnections combined in the GAO. When analyzing the cross-correlation values in Table 1, it is necessary to take into account the bandpass filtering from 2 to 7 years, previously applied to the SAT anomalies, from which EOF1 was then calculated. Due to the application of this filtering, the correlation values turned out to be somewhat overestimated, and rather allow for a relative assessment for each model of the degree of connection between SAT EOF1 and the EONI in comparison with other CMIP6 models under study, than to determine the contribution of ENSO to SAT EOF1.

Table 1 shows the maximum cross-correlations of EOF1 of the global interannual variability of the SLP anomalies (EOF1 SLP) and the GAO1 index (column 7) for shifts from −60 to +60 months, as well as the shifts to which these values correspond (column 8). Due to the long series from which they were calculated, all of these correlations turned out to be statistically significant at the 95% level. The inter-model scatter of the cross-correlations of EOF1 SLP and GAO1 is higher than for EOF1 SAT and the EONI, indicating that CMIP6 model ensemble reproduces the planetary structure of the GAO SLP anomalies worse than it does the El Niño’s SAT anomalies in the equatorial region of the central-eastern Pacific Ocean. However, the shifts at which the maximum correlations between EOF1 SLP and GAO1 are observed are equal to 0 for all 50 CMIP6 models studied, which indicates the synchronicity of EOF1 SLP and GAO1.

Thus, the obtained ensembled planetary spatial structures of the first principal components of the global interannual variability of the SAT and SLP anomalies correspond with a small lag (1–3 months) to the maximum phase of the El Niño events and are quasi-synchronous with the positive phase of the GAO. The correspondence of the ensembled PC1 SAT and PC1 SLP to the positive phases of the ENSO and GAO is explained by the fact that, if necessary, the principal components were taken with the opposite sign for some models. Since bandpass filtering was applied only to the SAT and SLP anomalies before calculating the principal components, and bandpass filtering was not applied to the EONI and GAO1 indices, the high values of the maximum cross-correlations at small shifts (Table 1) confirm that the ENSO and GAO are the first principal components of the global interannual variability of the SAT and SLP anomalies.

To assess the ability to forecast El Niño and La Niña events by the CMIP6 models under study, the ensembled second principal components (PC2—spatial structures of EOF2) of the global interannual variability of the SAT anomalies (PC2 SAT) (Figure 1b) and the SLP (PC2 SLP) (Figure 2b) are of great importance. They were constructed similarly to the 1st components, but the second principal components of 28 out of 50 CMIP6 models under study were averaged. These 28 models were a subset of the complete set of 50 model members, and they were selected as follows. The cross-correlation values between the PGAO (the ENSO predictor index) and EOF2 were calculated for the SAT and SLP anomalies (

Table 2. Cross-correlations of the 1st and 2nd empirical orthogonal functions of the interannual oscillations of the SAT and SLP anomalies with each other, as well as with PGAO.

Model Name	Max. Corr. EOF1 SAT and EOF1 SLP	Shift of EOF1 SAT and EOF1 SLP (Mon.)	Max. Corr. EOF2 SAT and EOF2 SLP	Shift of EOF2 SAT and EOF2 SLP (Mon.)	Max. Corr. EOF1 SAT and EOF2 SLP	Shift of EOF1 SAT and EOF2 SLP (Mon.)	Max. Corr. EOF1 SLP and EOF2 SLP	Shift EOF1 SLP and EOF2 SLP (Mon.)	Corr. EOF2 SAT and PGAO	Corr. EOF2 SLP and PGAO	Is the Model Included in the PC2 Calculation (95% Significance)?
TaiESM1	0.88	1	0.67	1	0.60	10	0.43	7	0.75	0.75	yes
AWI-CM-1-1-MR	0.81	1	0.07	0	0.58	5	0.16	1	0.33	0.04	no
BCC-CSM2-MR	0.82	1	0.32	0	0.42	1	0.35	0	0.03	0.08	no
BCC-ESM1	0.75	1	0.29	0	0.19	0	0.63	0	0.02	0.08	no
CAMS-CSM1-0	0.83	2	0.11	0	0.65	5	0.22	0	0.35	0.01	no
CAS-ESM2-0	0.83	1	0.14	4	0.58	5	0.45	0	0.60	0.08	yes
FGOALS-f3-L	0.85	1	0.22	0	0.74	6	0.40	0	0.34	0.10	yes
FGOALS-g3	0.87	1	0.16	0	0.71	6	0.45	0	0.36	0.05	no
CanESM5	0.75	2	0.09	1	0.44	5	0.21	0	0.19	0.05	no
CanESM5-CanOE	0.77	2	0.08	−5	0.51	6	0.17	0	0.30	0.06	no
CMCC-CM2-SR5	0.86	1	0.42	1	0.51	11	0.25	−1	0.76	0.62	yes
CMCC-ESM2	0.89	1	0.65	1	0.68	11	0.47	15	0.57	0.86	yes
CNRM-CM6-1	0.77	1	0.42	1	0.53	5	0.15	0	0.31	0.14	yes
CNRM-CM6-1-HR	0.49	2	0.08	52	0.52	6	0.58	0	0.32	0.06	no
CNRM-ESM2-1	0.79	1	0.22	1	0.43	5	0.09	5	0.46	0.11	yes
ACCESS-CM2	0.84	1	0.55	1	0.57	6	0.43	−1	0.47	0.51	yes
ACCESS-ESM1-5	0.79	1	0.15	1	0.42	−11	0.59	0	0.11	0.03	no
E3SM-1-0	0.85	1	0.09	0	0.54	6	0.33	0	0.63	0.05	no
EC-Earth3	0.76	1	0.11	0	0.39	−10	0.06	2	0.30	0.01	no
EC-Earth3-AerChem	0.75	1	0.35	0	0.35	8	0.21	0	0.51	0.15	yes
EC-Earth3-Veg	0.82	1	0.14	0	0.41	6	0.16	0	0.45	0.05	no
EC-Earth3-Veg-LR	0.80	1	0.15	1	0.37	7	0.22	0	0.63	0.02	no
FIO-ESM-2-0	0.86	1	0.18	0	0.54	9	0.27	0	0.78	0.19	yes
MPI-ESM-1-2-HAM	0.76	1	0.06	17	0.47	6	0.55	0	0.26	0.00	no
INM-CM4-8	0.53	1	0.08	6	0.22	−6	0.66	0	0.24	0.09	yes
INM-CM5-0	0.58	1	0.04	2	0.44	0	0.38	0	0.32	0.03	no
IPSL-CM6A-LR	0.77	2	0.15	1	0.67	5	0.09	−3	0.38	0.04	no
KIOST-ESM	0.81	1	0.26	0	0.68	7	0.45	0	0.18	0.02	no
MIROC-ES2L	0.88	3	0.09	20	0.74	10	0.11	0	0.35	0.04	no
MIROC6	0.87	4	0.07	13	0.67	9	0.06	2	0.51	0.05	no
HadGEM3-GC31-LL	0.85	1	0.26	1	0.52	7	0.46	0	0.62	0.18	yes
HadGEM3-GC31-MM	0.84	1	0.16	6	0.57	7	0.28	0	0.51	0.20	yes
MPI-ESM1-2-HR	0.80	1	0.18	0	0.49	8	0.34	0	0.47	0.14	yes
MPI-ESM1-2-LR	0.77	1	0.35	0	0.53	8	0.56	0	0.43	0.17	yes
MRI-ESM2-0	0.86	1	0.09	1	0.58	7	0.18	0	0.64	0.04	no
GISS-E2-1-G	0.86	1	0.64	1	0.74	8	0.42	8	0.10	0.63	yes
GISS-E2-1-H	0.81	1	0.27	0	0.61	7	0.47	0	0.27	0.05	no
CESM2	0.90	1	0.56	1	0.63	7	0.37	9	0.67	0.75	yes
CESM2-FV2	0.91	1	0.72	1	0.70	7	0.54	8	0.59	0.85	yes
CESM2-WACCM	0.89	1	0.31	1	0.61	7	0.26	0	0.60	0.52	yes
CESM2-WACCM-FV2	0.92	1	0.79	1	0.71	7	0.56	8	0.57	0.84	yes
NorCPM1	0.86	1	0.44	0	0.58	−11	0.33	0	0.46	0.23	yes
NorESM2-LM	0.88	1	0.36	1	0.65	6	0.31	0	0.23	0.43	yes
NorESM2-MM	0.87	1	0.09	−2	0.63	6	0.27	0	0.29	0.14	yes
GFDL-CM4	0.82	1	0.13	0	0.53	6	0.14	0	0.62	0.10	yes
GFDL-ESM4	0.86	1	0.12	0	0.53	7	0.15	1	0.68	0.15	yes
NESM3	0.79	1	0.24	0	0.28	7	0.10	1	0.53	0.10	yes
SAM0-UNICON	0.89	1	0.41	1	0.61	9	0.34	−1	0.74	0.64	yes
CIESM	0.89	1	0.19	0	0.60	8	0.19	0	0.54	0.00	no
MCM-UA-1-0	0.79	1	0.56	0	0.25	14	0.14	1	0.31	0.27	yes

, columns 10 and 11). The fact that the PGAO leads the El Niño index by about 1 year was shown in [36]. If the values of these correlations of the analyzed model had a statistical significance of more than 95%, then this model was included in the ensemble for averaging the two principal components (Table 2, column 12). The reason for choosing such a criterion is that not all of the 50 studied CMIP6 models have a second principal component of global interannual oscillations of SAT and SLP anomalies associated with the GAO and leading the phase of El Niño events. Due to the criterion of the presence of a significant connection of EOF2 with the PGAO and the PGAO leading the EONI and GAO1 indices, those 28 models were selected whose second principal component corresponds to the GAO phase preceding the El Niño events. At the same time, in order for the PGAO index to lead the positive phase of the ENSO, if necessary, the principal components in some models were taken with the opposite sign.

The average-model PC2 SAT (Figure 1b) and PC2 SLP (Figure 2b) demonstrate spatial structures that correspond to the initial phase of El Niño events. The main difference between PC2 SAT and PC1 SAT is the presence of negative values in the tropics of the Atlantic and Indian Oceans. It is above these negative SST anomalies that, according to [33], an extensive area of positive SLP anomalies, which is part of the GAO, is formed at the beginning of El Niño events. At the same time, in the central and eastern parts of the equatorial region of the Pacific Ocean, a structure (“tongue”) of positive SAT anomalies characteristic of El Niño is already beginning to form (Figure 1b). In the high latitudes of the Pacific Ocean, positive SAT anomalies are observed in the PC2 SAT field, as well as in the PC1 SAT field, but in the PC2 SAT field, they are located farther to the west than in the PC1 SAT field, which indicates a west–east spread of the spatial structure of the GAO. The PC2 SLP structure is also located to the west of the PC1 SLP structure, which confirms the west–east dynamics of the GAO.

Table 2 presents the maximum absolute values of mutual correlations of the 1st and 2nd empirical orthogonal functions (EOF1 and EOF2) of the interannual oscillations of the global SAT and SLP anomalies. The shifts corresponding to these maximum values are also given. The correlations between EOF1 SAT and EOF1 SLP, with some exceptions, turned out to be quite high (the average correlation is 0.81) with small phase shifts (1–2 months) (Table 2, columns 2 and 3). This indicates strong relationships between the first principal components of the interannual oscillations of the global SAT and SLP anomalies, and a slight lead of EOF1 SLP relative to EOF1 SAT. Since SAT over the water surface strongly depends on the SST, which has a fairly high inertia, and the SLP is an atmospheric characteristic, then the lead of EOF1 SLP relative to EOF1 SAT indicates the leading role of the atmosphere relative to the ocean on interannual time scales. At the same time, there is reason to believe that the higher the correlation between EOF1 SAT and EOF1 SLP, the stronger and more stable the global teleconnections between the atmosphere and the ocean in the CMIP6 model under consideration on interannual periods characteristic of the ENSO and GAO.

The cross-correlation values between EOF2 SAT and EOF2 SLP (Table 2, columns 4 and 5) demonstrate a much larger inter-model spread than between EOF1s. Some models show a fairly strong relationship between EOF2 SAT and EOF2 SLP, with correlation values between them greater than 0.5 (highlighted in bold). However, at the same time, there are also models in which the correlations between EOF2 SAT and EOF2 SLP are close to zero, which indicates the absence of significant relationships between them in these models. Thus, the first principal components of the interannual oscillations of the global anomalies of SAT and SLP turned out to be quite strongly related to each other in almost all of the 50 CMIP6 models under study, and the second ones in a minority of the 50 CMIP6 models.

Since the ENSO is an atmosphere–ocean interaction process, and the PGAO predictor index depends on both SAT and SLP anomalies, the CMIP6 models that have strong global teleconnections of SAT EOF1 and SLP EOF1 (Table 2, column 2), as well as SAT EOF2 with SLP EOF2 (highlighted in bold in Table 2, column 4), are better suited for predicting El Niño and La Niña events based on the GAO. Moreover, in such models, the first principal components of the global interannual oscillations of the SAT and SLP anomalies should have strong connections with the EONI and GAO1 (Table 1, columns 5–8), and the second principal components of these oscillations should have strong connections with the PGAO (highlighted in bold in Table 2, columns 10 and 11).

Moreover, for the long-term ENSO forecast, the models should have such reliable predictors that it is possible to forecast strong El Niño and La Niña events with a lead time of approximately 1 year, i.e., allowing one to overcome the SPB. Consideration of the maximum cross-correlations between SAT EOF1 and SAT EOF2 (Table 2, column 6), as well as the shifts at which these correlations are observed (Table 2, column 7), allows us to conclude that SAT EOF2 of some CMIP6 models (highlighted in bold in Table 2, columns 6 and 7) is ahead of SAT EOF1 by 7–11 months, with correlation values greater than 0.5. For these models, the second principal component of the interannual variability of SAT anomalies can be used to forecast the first principal component, i.e., to forecast the ENSO. The same applies to the cross-correlations between EOF1 SLP and EOF2 SLP (Table 2, columns 8 and 9), but the relationships between the first principal components of the SLP were weaker than for the SAT.

Thus, based on the analysis of the cross-correlations indicated in Table 1 and Table 2, from the 50 CMIP6 models under study, it is possible to select those that have good GAO-based ENSO predictive potential. To such, a subset of CMIP6 models can be preliminarily attributed: AS-RCEC TaiESM1, CMCC CMCC-ESM2, NASA-GISS GISS-E2-1-G, NCAR CESM2-FV2 and NCAR CESM2-WACCM-FV2 (highlighted in bold in Table 2, column 1).

To determine the spatial structures of the SAT and SLP anomalies that occur before El Niño events and that can be used to forecast the ENSO with different lead times, average compositional fields of the SAT (Figure 3) and SLP (Figure 4) anomalies with shifts of −14, −12, −10, −8, −6, −4, −2, and 0 months from the moment of maximum development of El Niño events were calculated for 50 studied CMIP6 models. These fields demonstrate the west–east distribution of planetary structures of the SAT and SLP anomalies, which is typical for the GAO [33,35]. Due to this west–east distribution, the PGAO predictor index, which describes the patterns of SAT and SLP anomalies occurring approximately 12 months before the El Niño maximum (Figure 3b and Figure 4b), can predict the ENSO with approximately a year’s lead time [36].

Figure 5 shows the average composite fields of SAT and SLP anomalies based on the PGAO index for 50 CMIP6 models. These inter-model mean spatial structures are characteristic of the GAO phase that precedes the El Niño maximum by approximately 12 months. Green and purple colors in Figure 5 highlight the areas in which the average anomalies in SAT and SLP participate in the PGAO calculation with a positive and negative sign, respectively. It is evident that the areas marked in Figure 5 cover a significant part of the Earth’s surface. It is also possible to note the high similarity of the PC2 SAT field (Figure 1b) with the composite field of SAT anomalies that precede the El Niño maximum by 6 months (Figure 3d), and the PC2 SLP field (Figure 2b) with the composite field of SLP anomalies based on the PGAO index (Figure 5b). Thus, inter-model mean PC2 SLP is ahead of the maxima of El Niño events significantly more than PC2 SAT, which indicates the leading role of the atmosphere in this process.

Since the PGAO includes both SAT and SLP anomalies, in order to estimate the degree of relationships between the second principal components and the PGAO for each of the 50 CMIP6 models under study separately, the maximum values of cross-correlations between the sum of the absolute values of SAT EOF2 and SLP EOF2 (combined EOF2) with PGAO (

Table 3. Cross-correlations with shifts of the EONI, PGAO indices and their EOF1 and EOF2 combined for SAT and SLP.

Model Name	Max. Corr. of EOF2 and PGAO	EOF2 and PGAO Shift (mon.)	Max. Corr. of EOF2 and EOF1	EOF2 and EOF1 Shift (mon.)	Max. Corr. EOF2 and EONI	EOF2 and EONI Shift (mon.)	Max. Corr. PGAO and EONI Filter from 2 to 7 Years	PGAO and EONI Shift (mon.)	Coherence of PGAO and EONI Periods from 2 to 7 Years	Phase Relationships of PGAO and EONI (°)
TaiESM1	0.82	0	0.60	−9	0.64	−5	0.81	−11	0.88	−109.5
AWI-CM-1-1-MR	0.19	0	0.36	−6	0.40	−2	0.52	−9	0.73	−40.8
BCC-CSM2-MR	0.16	3	0.49	−1	0.35	0	0.36	−8	0.67	−24.6
BCC-ESM1	0.13	−4	0.18	0	−0.07	0	0.36	−6	0.63	−16.8
CAMS-CSM1-0	0.39	4	0.55	−4	0.54	−2	0.72	−11	0.83	−123.9
CAS-ESM2-0	0.38	1	0.40	−5	0.41	−4	−0.67	10	0.80	−126.5
FGOALS-f3-L	0.39	6	0.59	−6	0.62	−4	0.81	−12	0.76	−129.3
FGOALS-g3	0.30	3	0.59	−4	0.64	−2	0.61	−9	0.66	−143.9
CanESM5	0.20	1	0.52	−1	0.50	−2	0.41	−8	0.69	−34.5
CanESM5-CanOE	0.21	2	0.33	−4	0.39	−2	0.50	−10	0.68	−35.1
CMCC-CM2-SR5	0.81	0	0.44	−12	0.47	−7	0.66	−11	0.79	−123.2
CMCC-ESM2	0.79	0	0.64	−13	0.61	−9	0.78	−14	0.87	−121.7
CNRM-CM6-1	0.30	1	0.44	−5	0.47	−3	0.60	−11	0.73	−143.2
CNRM-CM6-1-HR	0.24	2	0.42	0	0.39	−3	−0.61	7	0.77	−137.5
CNRM-ESM2-1	0.41	1	0.40	−5	0.43	−3	−0.59	8	0.76	−141.9
ACCESS-CM2	0.59	1	0.49	−6	0.51	−5	0.65	−10	0.81	−136.9
ACCESS-ESM1-5	0.07	25	0.27	12	−0.35	−1	0.43	−13	0.71	−145.0
E3SM-1-0	0.41	1	0.35	−7	0.37	−4	0.68	−9	0.82	−58.0
EC-Earth3	0.24	1	−0.30	12	0.30	−3	0.58	−10	0.75	−47.7
EC-Earth3-AerChem	0.41	0	0.26	−9	0.29	−6	0.56	−9	0.71	−33.4
EC-Earth3-Veg	0.33	2	0.40	−5	0.45	−3	0.58	−11	0.73	−46.3
EC-Earth3-Veg-LR	0.41	0	0.28	−6	0.31	−5	0.50	−9	0.70	−43.5
FIO-ESM-2-0	0.63	0	0.51	−7	0.56	−5	0.75	−10	0.82	−62.0
MPI-ESM-1-2-HAM	0.23	1	0.56	−1	0.52	−2	0.46	−13	0.70	−143.5
INM-CM4-8	0.12	1	0.36	1	−0.15	−7	−0.17	9	0.58	−22.5
INM-CM5-0	0.24	0	0.30	0	0.23	−1	0.26	−11	0.64	−21.4
IPSL-CM6A-LR	0.23	0	0.43	−5	0.48	−2	0.40	−9	0.69	−31.8
KIOST-ESM	0.32	8	0.54	−5	0.56	−4	−0.80	6	0.84	−134.8
MIROC-ES2L	0.37	8	0.54	−10	0.50	−9	0.75	−17	0.82	−137.4
MIROC6	0.36	1	0.56	−7	0.57	−3	0.69	−14	0.85	−128.7
HadGEM3-GC31-LL	0.51	0	0.35	−7	0.39	−5	0.63	−10	0.81	−127.4
HadGEM3-GC31-MM	0.49	1	0.41	−7	0.44	−6	−0.70	7	0.81	−133.6
MPI-ESM1-2-HR	0.41	1	0.45	−6	0.46	−5	−0.56	7	0.76	−138.1
MPI-ESM1-2-LR	0.37	1	0.47	−1	0.45	−4	0.54	−13	0.74	−142.2
MRI-ESM2-0	0.49	1	0.41	−7	0.42	−5	0.72	−11	0.82	−127.0
GISS-E2-1-G	0.45	1	0.66	−8	0.66	−8	0.81	−18	0.87	−146.5
GISS-E2-1-H	0.31	11	0.52	−1	0.49	−4	−0.79	6	0.88	−134.1
CESM2	0.81	0	0.58	−8	0.56	−6	0.74	−10	0.83	−57.2
CESM2-FV2	0.78	0	0.69	−8	0.65	−5	0.79	−10	0.80	−46.8
CESM2-WACCM	0.69	0	0.52	−8	0.50	−5	0.70	−10	0.84	−58.1
CESM2-WACCM-FV2	0.76	1	0.69	−8	0.65	−6	0.81	−11	0.86	−125.2
NorCPM1	0.46	1	0.44	−7	0.51	−4	−0.71	9	0.86	−119.1
NorESM2-LM	0.41	1	0.50	−6	0.55	−4	0.75	−13	0.83	−122.2
NorESM2-MM	0.31	7	0.58	−4	0.66	−2	0.68	−11	0.79	−130.0
GFDL-CM4	0.50	1	0.49	−4	0.51	−3	0.58	−10	0.76	−47.4
GFDL-ESM4	0.56	0	0.50	−5	0.52	−4	0.64	−11	0.74	−132.6
NESM3	0.28	0	0.16	−6	0.20	−4	−0.53	7	0.71	−146.4
SAM0-UNICON	0.84	0	0.54	−10	0.54	−6	0.74	−9	0.86	−57.6
CIESM	0.41	2	0.52	−7	0.55	−4	0.75	−13	0.84	−121.8
MCM-UA-1-0	0.33	0	0.21	−14	0.22	−12	−0.45	8	0.70	−143.4

, column 2), and the shifts to which they correspond (Table 3, column 3) were calculated. It is evident that mainly the TaiESM1, CMCC-ESM2, GISS-E2-1-G, CESM2-FV2, and CESM2-WACCM-FV2 models marked in bold in Table 2, column 1 are characterized by fairly high values of these correlations, with a shift close to zero. In addition to the above subset of CMIP6 models, we can add some others that also have high and, at the same, time quasi-synchronous connections between the combined EOF2 and PGAO (highlighted in bold in Table 3, columns 2 and 3): CMCC-CM2-SR5, ACCESS-CM2, FIO-ESM-2-0, CESM2, CESM2-WACCM, GFDL-ESM4, SAM0-UNICON. Thus, there is a gradual expansion of the subset of CMIP6 models that have the potential to predict El Niño and La Niña events in advance. It is important that the models falling into this subset based on the correlations in Table 3 also have fairly high correlations in Table 1 and Table 2.

Some CMIP6 models have strong links between combined EOF1 and EOF2 (Table 3, column 4), with EOF2 leading by more than half a year (Table 3, column 5) marked in bold in these columns. CMIP6 models are also found to have strong links between combined EOF2 and EONI (Table 3, column 6), with EOF2 leading by half a year (Table 3, column 7) marked in bold in these columns. Also of major importance for the GAO-based ENSO lead time forecast is the presence of a strong link between the PGAO and EONI (Table 3, column 8) with the PGAO leading by approximately 1 year (Table 3, column 9). CMIP6 models with such couplings were further identified by the coherence between the PGAO and EONI (Table 3, column 10) and the corresponding phase relationships (Table 3, column 11) for oscillation periods typical for the ENSO (2–7 years), and are marked in bold in Table 3, columns 1 and 8–11. For the models with the highest correlation values between the PGAO and EONI (highlighted in bold), the shift at which this value is reached varies from −17 to −9 months. Thus, different CMIP6 models have different ENSO forecast lead times based on the west–east distribution of the GAO planetary structure. This is important for overcoming the SPB of ENSO events. The above parameters can be used for assessing the ENSO’s prediction skills of the individual CMIP6 models.

As shown in [37], some CMIP6 models reproduce the global structure of ENSO teleconnections better than other models. Due to this, as well as the reproduction of the west–east propagation of this global structure, some of the CMIP6 models have the ability to predict the ENSO with greater lead times than other ones. All the above-listed relative estimates of the relationships allowed us to select the following 21 CMIP6 models that can successfully forecast ENSO from the GAO with a lead time of about 1 year: TaiESM1, CAMS-CSM1-0, FGOALS-f3-L, CMCC-CM2-SR5, CMCC-ESM2, ACCESS-CM2, FIO-ESM-2-0, MIROC-ES2L, MIROC6, HadGEM3-GC31-LL, MRI-ESM2-0, GISS-E2-1-G, CESM2, CESM2-FV2, CESM2-WACCM, CESM2-WACCM-FV2, NorESM2-LM, NorESM2-MM, GFDL-ESM4, SAM0-UNICON, and CIESM. The graphs of the cross-correlation functions with shifts from −60 to +60 months of the EONI and PGAO indices for most of these models, and the average for 50 CMIP6 models are shown in Figure 6. These graphs allow us to state that the PGAO leads the EONI by approximately 1 year for most CMIP6 models, but there is a noticeable difference in the ability to forecast the ENSO using the GAO between the models. Thus, the cross-correlations with a shift of approximately −1 year for the 21 models listed above are higher than the inter-model mean of all 50 CMIP6 models under study. We did not display the cross-correlation graphs with high values for models from the same developers, for example from NCAR, in order not to clutter Figure 6. It should also be kept in mind that some models not included in this subset of 21 models may also have fairly good predictive ability of ENSO, but are inferior in the lead time.

The oscillatory shape of the cross-correlation functions presented in Figure 6 is noteworthy, with an average oscillation period of approximately 45 months (3.7 years). This corresponds to the maximum of the average model EONI spectrum found in [37] at a period of approximately 4 years. Thus, based on the obtained results, it can be concluded that the EONI, GAO1 and PGAO are indices of the same planetary process, called the GAO in [33], and including the ENSO as one of its elements. Moreover, an essential part of the 50 CMIP6 models considered reproduces both the planetary structure of the GAO itself and its west–east distribution, due to which it is possible to predict the ENSO based on these models, overcoming the SPB. The fields of cross-correlations between the EONI and SAT and SLP anomalies with shifts from −21 to +21 months (43 months at all) demonstrate planetary spatiotemporal dynamics of the GAO (Figure 7 and Figure 8). A close correspondence of the spatial structure of correlations at the −12 month shift with the SAT and SLP anomaly fields preceding the 12 months before the El Niño culmination is evident (Figure 3b and Figure 4b). The same correspondence can be traced in subsequent fields up to the shift of 0 months. It is also evident that the spatial structure of the GAO as a whole extends from west to east.

The west–east propagation of the spatial structure of the GAO in the CMIP6 ensemble of models is consistent with previously obtained results from observational data [33]. This west–east propagation is most clearly seen in a sequence of grid-point cross-correlations between the EONI and the SLP anomalies (Figure 8). First (at the month “−21 mon” that corresponds to La Niñas), the X-shaped structure on the Pacific is formed by positive correlations, and the tropical elliptic-shaped structure is formed by negative correlations. Then, both structures move eastward to the American continent and West Pacific, respectively. At about the month “−12 mon”, these structures start to transform themselves into their opposites, i.e., the X-shaped structure becomes more similar to the elliptic-shaped structure, and the elliptic-shaped structure becomes similar to the X-shaped one. These transformations end at the month “0 mon”, which corresponds to the maximum of the El Niño event. Analogous eastward movements and transformations take place during the subsequent months up to the month “+21 mon”, which corresponds to the maximum of La Niña event. Thus, the spatial dynamics of the GAO represent a transition from La Niño to El Niño and back, and demonstrate the wave nature of this process on a planetary scale. The use of the planetary nature of this wave process makes it possible to overcome the SPB in ENSO forecasts.

Thus, the CMIP6 ensemble of models as a whole reproduces those features of the west–east dynamics of planetary structures of SAT and SLP anomalies that were discovered earlier from observational data [33,35,36]. And ranking the models using correlations between different ENSO and GAO indices, as well as the EOFs of interannual variations in SAT and SLP anomalies, allows us to find the subset of models that best reproduce these features. The graphs in Figure 6 show that the values of the cross-correlation functions of the founded models exceed the CMIP6 ensemble mean at a shift of approximately −12 months. This suggests that using this subset of models may provide better results in predicting El Niño and La Niña events before SPB than using the average of the entire CMIP6 ensemble of models. It should be noted that among all other models, MIROC-ES2L and GISS-E2-1-G stand out, having maximum correlations at shifts of −17 and −18 months, respectively (Figure 6, Table 3).

The described structure of the GAO and its west–east dynamics is global. But one can find many of its regional manifestations. Thus, approximately the year after El Niño, the GAO SLP anomaly structure (Figure 8, “+18 mon”) resembling the positive phase of the North Atlantic Oscillation is formed in the North Atlantic [42]. Of interest for the Southern Ocean may be the patterns of positive correlations with SST. One of these patterns extends from west to east from the region south of Australia and New Zealand (Figure 7, “−21 mon”) to the Drake Passage (Figure 7, “+21 mon”) across the entire Pacific sector of the Southern Ocean. The EONI cross-correlation fields with SLP anomalies also show a west–east spread of a structure reminiscent of the Antarctic Oscillation pattern [43]. And many similar local manifestations of ENSO teleconnections can be found in many different regions of the earth [35]. Moreover, in [35], a review of the reverse influence of anomalies in regions outside the tropical Pacific on the ENSO is also made.

As a validation of the ENSO forecasts, the EONI and PGAO indices for the abstract 115 years of the CMIP6 models TaiESM1, CMCC-ESM2, FGOALS-f3-L, MIROC-ES2L, GISS-E2-1-G, and NorESM2-LM are shown as an example (Figure 9). Also, the Oceanic Niño Index (ONI) and PGAO graphs for 1900–2014 of the Historical experiment of the TaiESM1, CMCC-ESM2, FGOALS-f3-L, MIROC6, CESM2-FV2, and NorESM2-LM models are shown in Supplementary Materials (Figure S7). These graphs show that almost everywhere, the PGAO maxima lead the EONI maxima (El Niño events), and the PGAO minima lead the EONI minima (La Niña events). This lead is especially characteristic of strong events. This allows using PGAO as an ENSO predictor in the selected CMIP6 models. However, there are periods when El Niño and La Niña events weaken, and the predictive ability of the PGAO decreases. But this happens mostly when there are weak El Niño and La Niña events. Strong El Niño and La Niña events are generally preceded by PGAO maxima and minima, which confirms the possibility of forecasting these events by the selected CMIP6 models with a long lead time.

In [39], based on observational data, it was shown that the ENSO is affected by the following external quasi-periodic forces: the annual variation in heat input from the sun, the 14-month Chandler wobble in the Earth’s pole motion, the 11-year variation in solar activity, and the 18.6-year luni-solar nutation of the Earth’s rotation axis. Due to the incommensurability of the periods of action of these external forces, there is reason to believe that the ENSO dynamics are not chaotic, but can be described by the mathematical model of the strange nonchaotic attractor (SNA). Thus, the ENSO dynamics can be predicted by statistical methods with a sufficiently long lead time.

4. Conclusions

The interannual variability of global SAT and SLP anomalies of an ensemble of CMIP6 models is investigated for its relationship with ENSO. The results of the piControl experiment for a long model period of about 500 years without changes in external forcing, and the Historical experiment with changing radiative forcing from the anthropogenic greenhouse effect for the period 1850–2014, are used.

Based on three different methods—principal component analysis, compositional analysis and correlation analysis—it is shown that, in general, for the considered ensemble of 50 models, the first principal component of interannual variations in SAT and SLP anomalies is associated with the phase of maximum development of El Niño and La Niña events. It also turned out that in 28 out of 50 models considered, the second principal component corresponds to a phase that is ahead of the maximum of El Niño events.

It is shown that the ensemble of 50 CMIP6 models quite accurately reproduces the west–east distribution of the planetary spatial structure of the amplitude of oscillations of the SAT and SLP anomalies characteristic of the ENSO and called the GAO, as well as its periodicity (~3.7 years). Thus, it can be concluded that the most of considered CMIP6 models contain a physical mechanism for the interaction of the ENSO with the Indian and Atlantic oceans, as well as with temperate and high latitudes.

The degrees of EOF1 and EOF2 connections of interannual variability of SAT and SLP anomalies with the ENSO and GAO indices are estimated. The CMIP6 models with strong EOF1 links with the ENSO index EONI and the GAO index GAO1 are determined. The CMIP6 models with combined EOF2 having strong links with the PGAO predictor index are also found, which can be used to predict El Niño and La Niña events in advance based on the west–east spread of the GAO spatial structure. Based on the analysis of the identified relationships, the following 21 CMIP6 models were selected that allow ENSO forecasting with a lead time of approximately 1 year: AS-RCEC TaiESM1, CAMS CAMS-CSM1-0, CAS FGOALS-f3-L, CMCC CMCC-CM2-SR5, CMCC CMCC-ESM2, CSIRO-ARCCSS ACCESS-CM2, FIO-QLNM FIO-ESM-2-0, MIROC MIROC-ES2L, MIROC MIROC6, MOHC HadGEM3-GC31-LL, MRI MRI-ESM2-0, NASA GISS-E2-1-G, NCAR CESM2, NCAR CESM2-FV2, NCAR CESM2-WACCM, NCAR CESM2-WACCM-FV2, NCC NorESM2-LM, NCC NorESM2-MM, NOAA GFDL-ESM4, SNU SAM0-UNICON and THU CIESM.

Thus, based on the ensemble of CMIP6 models, it is shown that the ENSO is an element of a planetary phenomenon— the GAO. The GAO combines ENSO teleconnections with the SAT and SLP anomalies outside the tropics of the Pacific Ocean. The GAO structure spreads from west to east; that is, it has the property of a planetary wave, making a full revolution around the Earth in about 3.7 years. From the entire CMIP6 ensemble, those models were selected that reproduce the above-mentioned features better than others. That makes it possible to improve the advance forecast of the ENSO base on this subset of CMIP6 models.