How Does El Niño–Southern Oscillation Change Under Global Warming—A First Look at CMIP6

The latest generation of coupled models, the sixth Coupled Models Intercomparison Project (CMIP6), is used to study the changes in the El Niño–Southern Oscillation (ENSO) in a warming climate. For the four future scenarios studied, the sea surface temperature variability increases in most CMIP6 models, but to varying degrees. This increase is linked to a weakening of the east‐west temperature gradient in the tropical Pacific Ocean, which is evident across all models. Just as in previous generations of climate models, we find that many characteristics of future ENSO remain uncertain. This includes changes in dominant time scale, extratropical teleconnection patterns, and amplitude of El Niño and La Niña events. For models with the strongest increase in future variability, the majority of the increase happens in the Eastern Pacific, where the strongest El Niño events usually occur.


Introduction
El Niño-Southern Oscillation (ENSO) is characterized by irregular fluctuations between cold (La Niña) and warm (El Niño) conditions in the eastern and central equatorial Pacific on a time scale of 2-7 years. The warm phase is associated with a weakening of the trade winds and eastward shift of convection, which brings the warm waters of the west Pacific eastward. This decrease of the east-west gradient in sea surface temperature (SST) is concomitant with a deepening (shoaling) of the thermocline in the eastern (western) equatorial Pacific. Due to its global teleconnections, ENSO is not only the dominant mode of tropical interannual variability but also the leading source of forecast skill on seasonal to interannual time scales in many other parts of the world (Barnston, 2016;Jin et al., 2008). It has important impacts on fisheries, agriculture, hurricanes, droughts, floods, and other severe weather events.
A rich body of work has studied the response of ENSO to global warming in previous generations of climate models, but there has been no clear consensus on how ENSO will change under global warming (e.g., Berner et al., 2020;Cai et al., 2015a;Collins et al., 2010;Guilyardi et al., 2012;Stevenson, 2012;Taschetto et al., 2014;Yeh et al., 2012). The CMIP6 archive provides a new opportunity to study the ENSO response to prescribed radiative forcing ) across a number of state-of-the-art climate models. Here, we will focus on the following question: "To which degree do CMIP6 models agree on ENSO changes in different global warming scenarios?" Assessments of ENSO future changes must account for the diversity of ENSO spatial patterns (Ashok et al., 2007;. ENSO events display a broad spectrum of anomaly centers ranging from the dateline (CP events) to the far eastern equatorial Pacific (EP events, Capotondi et al., , 2020, and the exact location of the warming centers may be model dependent . This diversity can have very important consequences for atmospheric teleconnections and worldwide impacts (Ashok et al., 2007;Larkin & Harrison, 2005;Patricola et al., 2018) and needs to be considered when examining ENSO response to global warming.
Several studies have shown that ENSO characteristics, such as period and growth rate are highly dependent on the tropical Pacific mean state (Battisti & Hirst, 1989;Fedorov & Philander, 2001). In particular, the mean temperature in the eastern equatorial Pacific, which controls the temperature gradient between the West Pacific warm pool and Eastern Pacific cold tongue, as well as the zonal slope of the equatorial thermocline, are important factors controlling ENSO's stability characteristics and ENSO diversity (Capotondi & Sardeshmukh, 2015;Fedorov & Philander, 2000, 2001. Specifically, a deeper thermocline in the eastern equatorial Pacific, accompanied by reduced easterly winds and weaker zonal SST gradient, favors longer periods and larger SST anomalies in the eastern equatorial Pacific, as observed, for instance, in the 1980s and 1990s (Capotondi & Sardeshmukh, 2017;Fedorov & Philander, 2001) relative to previous decades. Another controlling factor is cross-equatorial winds, which can significantly influence ENSO properties but with larger uncertainties in future scenarios (Hu & Fedorov, 2018). The ENSO mean state relationship is further complicated by the presence of ENSO asymmetries, with warm events typically stronger than cold events in the eastern Pacific (positive skewness), and cold anomalies somewhat larger than warm anomalies in the central Pacific (negative skewness), an aspect of ENSO that may be indicative of system nonlinearities. Such nonlinearities may, in turn, lead to a "rectification" of ENSO variations into the mean state, resulting in a low-frequency modulation of equatorial SSTs that are El Niño-like. Indeed, Karamperidou et al. (2017) find a significant relationship between ENSO amplitude changes and the correlation between the patterns of ENSO and SST trends.
Analyses of previous generations of climate models reported a weakening of the zonal SST gradient and of the atmospheric Walker circulation across the majority of the models, a consensus that did not translate, however, in a consistent change in ENSO amplitude, as measured by commonly used ENSO indices (e.g., the Niño 3.4 index) that are averages of SST anomalies at a fixed location. Given the differences in ENSO spatial patterns across models, a better model agreement was found, albeit for a selected group of models, when indices that accounted for the ENSO patterns unique to each model were used Carréric et al., 2019). The selection criterion was based on a metric of model nonlinearity, as encapsulated by the coefficient ɑ of the nonlinear relationship between the two leading Principal Components (PCs) of SST in the equatorial Pacific (Karamperidou et al., 2017). Models in the CMIP5 archive with a parameter ɑ in the preindustrial control simulations close to the observed value (−0.29) appeared to have a balance of (linear) ENSO feedbacks in better agreement with observations and exhibited a warming trend in the eastern equatorial Pacific (Karamperidou et al., 2017). The parameter ɑ also appeared to be associated with values of SST skewness in the eastern and central equatorial Pacific similar to the observed values, as well as a "realistic" separation of EP and CP ENSO events . An increase in ENSO amplitude was found by Cai et al. (2018) in those models with values of ɑ relatively close to the observed in the historical model simulations. This increase in amplitude was attributed to both the enhanced mean warming and increased vertical stratification in the eastern equatorial Pacific. In this study, we revisit the relationship between changes in ENSO amplitude and mean state changes in the latest generation of climate models with a primary focus on the connection between changes in ENSO amplitude and changes in the mean zonal SST gradient.
A rather robust response detected in previous generations of climate models is the projected poleward shift of the jet stream (Yin, 2005), which changes the atmospheric meridional gradients, and thus affects tropical-extratropical teleconnection patterns (Stevenson, 2012;Stevenson et al., 2012). Hence, extratropical teleconnection patterns may change under climate change scenarios, even if ENSO itself does not change significantly. This aspect is also examined in our analysis of the CMIP6 models.
To assess the degree of agreement among the CMIP6 models on the ENSO change under different warming scenarios, we use a number of established diagnostics. First, we investigate changes in power spectral densities, total variance, and zonal SST gradient. Then we examine variance changes in the context of ENSO diversity, and at last, we investigate global sea-level pressure teleconnection patterns during future El Niño and La Niña events.

Data and Methods
Our analysis focuses on the projected change under four different Shared Socioeconomic Pathways (SSP) (O'Neill et al., 2016). The four SSP scenarios are expected to have an approximate forcing of 2.6, 4.5, 7.0, and 8.5 W/m 2 in the year 2100, as denoted by the last two digits of the names of the scenarios. Results obtained for these four future scenarios (covering the 86-yr period 2015-2100) are compared to those from the control simulation (piControl). Since different models have a different number of ensemble members, for a fair comparison we use only one member for each future scenario.
However, due to a large level of internal variability, one ensemble member may be insufficient to robustly detect interscenario differences in variance. Hence we focus only on the significant changes in each scenario relative to the preindustrial control. Establishing significant variance changes between scenarios, or between scenarios and the historical period requires several ensemble members, for example, like the 33-member ensemble of CESM1 used by Berner et al. (2020).
The analyses are performed for 11 models having a control simulation with at least 499 years. For all piControl simulations longer than 500 years, we only study the first 500 years. Changes from the piControl are analyzed to focus on the models' response to anthropogenic forcing rather than evaluating model skill over the historical period. The confidence intervals of the piControl simulations used for analyses in Figures 1, 2, and 3 are computed by first splitting the 500-year records into 86-year segments with 56-year overlap. For each segment, we compute the quantity of interest in the same way as we would for the future scenarios and determine the range of possible values across segments. Changes from the piControl are considered statistically significant if they are outside this range.
Details regarding the detrending method, spectral analysis, definition of ENSO diversity, and ENSO teleconnections diagnostics are provided in Text S1 in the supporting information.  Figure S1 in the supporting information) are shown with a thick line. The green numbers denote the maximum power for the model MIROC-ES2L.

Power Spectral Density of Niño 3.4 Index
The Niño 3.4 index, defined as the area average of monthly SST anomalies in the region 5°S to 5°N, 170-120°W, is a commonly used index to describe variability associated with ENSO (see Text S1 for details on method). The power spectra of the Niño 3.4 index obtained from detrended scenarios show a wide range of variability with spectral peaks in the 2-7 year range (Figure 1 and Figure S1 in the supporting information), demonstrating the CMIP6 models' ability to produce a quasi-oscillatory behavior that is reminiscent of ENSO in nature. Figure 1 shows the seven models with the most marked changes in the spectra of future scenarios. The most significant increases in variability are for the models MIROC6 and MIROC-ES2L. MIROC-ES2L is the model with the most regular periodic variations, centered at periods of about 5 years ( Figure S1). CESM2-WACCM and CESM2 also show significant increases in power for most future scenarios with CESM2-WACCM showing the largest increase around 3 years and CESM2 showing an increase at 1.5-3 year periodicities. CanESM5 shows a small increase in power, mainly at the seasonal cycle and at periods of 3-4 years. The limited duration of the scenario simulations (86 years) makes it difficult to estimate subtleties in the future  5°N, 170-120°W] for each model and scenario after detrending and removing the seasonal cycle (see Text S1), shown as differences from the piControl estimate using the first 500 years (black line). (b) The black lines show the piControl mean of the east-west temperature gradient, and the colored bars the mean change from piControl over the 86-year period 2015-2100 in the future scenarios. Temperatures in the west are averaged over the region 5°S to 5°N, 120-170°E, and in the east over the Niño3 region [5°S to 5°N, 150-90°W]. Positive values mean the west is warmer than the east. In both panels, the light gray shading denotes the spread of these quantities in 86-year overlapping segments from the first 500 years of piControl, ranging from the minimum to the maximum estimates. (c) Scatterplot and correlation coefficient of data in (a) and (b). (d) As in (c), but with E-index variance from the next section along the x-axis. change of spectra. In particular, it is hard to assess whether the increase in power is proportional to the radiative forcing.

Changes in Variance of Niño 3.4 Index and Zonal Temperature Gradient
Another important metric to assess ENSO changes is the variance of the Niño 3.4 index. All but one model show a significant increase in variance for most scenarios (Figure 2a), but we observe no correlation between the magnitude of forcing and variance change. As the change in mean SST in the Niño 3.4 region is proportional to the forcing ( Figure S2), we conclude that there is no obvious relation between changes in variance and mean SST either.
The strength of ENSO variability has been linked to the east-west SST gradient in the Tropical Pacific. As this gradient weakens, westerly wind anomalies can more readily extend eastward and initiate strong warm events (Xie et al., 2018). Strikingly, all models agree that the east-west SST gradient weakens in future scenarios. Furthermore, in most models this weakening is proportional to the magnitude of the radiative forcing. This suggests an anticorrelation between the change in gradient and change in SST variance. However, the nonmonotonic increase in variance with increasing forcing suggests that there may also be other controlling mechanisms to be further investigated, such as, changes in thermocline depth, model sensitivity to aerosol forcing, changes in zonal, and meridional wind patterns.
Ten out of the 11 CMIP6 models show a decrease in SST gradient concurrent with an increase in SST variance ( Figure 2c) with a statistically significant (p-value = 1.52 · 10 −4 ) correlation of −0.55. We expect that some of the scatter in Figure 2c is due to the large internal ENSO variability (Berner et al., 2020), and the uncertainty in the functional relationship displayed in Figure 2c might be reduced if models with several ensemble members were analyzed.

ENSO Diversity
A single index, like the Niño 3.4 index, is insufficient to capture the full range of ENSO expressions and the temporal evolution of ENSO events. In particular, more than one index is needed to describe differences in ENSO spatial patterns. Several indices have been proposed to describe this diversity in El Niño spatial patterns (Capotondi et al., 2020). Here we use the approach introduced by Takahashi et al. (2011) to

Geophysical Research Letters
construct the E-and C-indices, which describe events with enhanced variability in the eastern and central Pacific, respectively (see Text S1 and Figure S3 for examples of their associated patterns). These indices are computed as linear combinations of the two leading principal components (PCs) of SST anomalies in the equatorial Pacific, and thus describe the patterns of variability typical of each model. Figure 3a shows the standard deviation of the E-index in the four climate change scenario simulations relative to the control simulation, whose standard deviation was normalized to one. Instead of arbitrarily selecting the models by their ɑ value in the control simulation, we include all the models and order them by increasing the magnitude of ɑ, to highlight the impact of this parameter on the changes in ENSO variance. While the only model showing a significant decrease in ENSO variance is the model with the smallest absolute value of ɑ, no clear relationship can be seen in Figure 3a between the variance changes and the magnitude of ɑ. In addition, the parameter ɑ may change in the scenario simulations relative to the control simulations (see Figure S5, where the nonlinear fit of PC1 and PC2 is shown for all models and simulations) and is not an intrinsic property of each model, as implied in the studies of Karamperidou et al. (2017) and Cai et al. (2018).
When the change of the variance in the E-index is plotted against the change of the east-west SST gradient, we see an even stronger relationship between the two quantities, as quantified by a correlation coefficient of −0.65 with p-value = 1.97 · 10 −6 ( Figure 2d): generally, a weakening of the SST gradient will lead to an increase in the variance of Eastern Pacific (EP) events.
Previous studies (Bellenger et al., 2014;Cai et al., 2015b;Capotondi, 2015;Stevenson et al., 2012) have also suggested an increase in the frequency of extreme La Niña events with global warming due to a strengthened zonal temperature gradient between the Maritime Continent and the central Pacific, where La Niña events typically peak. This should be reflected in the standard deviation of the C-index (Figure 3b). Robust increases in the C-index standard deviation are seen in some models, but in some cases the standard deviations show a significant decrease (BCC-CSM2-MR, CESM2, MCM-US-1-0) or insignificant changes (CanESM5, CESM2-WACCM, UKESM1-0-LL), indicating a larger degree of uncertainty in the projected changes of La Niña's (as well as CP El Niño) amplitude relative to the EP El Niño events.

ENSO Teleconnections
ENSO, primarily a tropical ocean-atmosphere coupled process, has an influence globally via atmospheric and oceanic teleconnections (Alexander et al., 2002;Deser et al., 2012;Yeh et al., 2018). The teleconnection diagnostics used in this study follows that of Stevenson et al. (2012). El Niño and La Niña composites for the ensemble mean of CMIP6 are computed for sea-level pressure (SLP) anomalies. The mean SLP anomalies of the ensemble mean show the canonical features of the Aleutian low deepening during warm events and anomalous higher pressures during cold events with pressures of opposite sign in the Southern Hemisphere at the same longitude. Teleconnection changes are then evaluated by comparing the ssp585 and piControl. The changes in the future climate across the CMIP6 ensemble are shown in the black contours in Figure 4. Marked changes in the Aleutian island region and the Southern Ocean region are observed. The atmospheric teleconnections show a weakening signal in the ssp585 scenario compared to the piControl. This has been studied in previous versions of similar climate models and partly been attributed to the increase in atmospheric static stability in warmer climates (Ma et al., 2012;Stevenson et al., 2012).
The spatial patterns of the teleconnection of warm events in the future scenario shifts poleward and eastward in the CMIP6 ensemble mean over the Aleutian island regions and the Southern Ocean regions. The eastward shift can be seen as a weakening of the pressure anomalies in the west, and a strengthening in the east in Figure 4a. This has also been seen in the model versions from CMIP5 and CMIP3 (Meehl & Teng, 2007;Stevenson et al., 2012). The teleconnection patterns for La Niña events show a zonal elongation over the Aleutian region instead of a spatial shift, and a general weakening of the Southern Ocean anomalies as seen in Figure 4b. The standard deviation across the ensemble members is shown in Figures 4c and 4d, where the largest variance is observed over the Aleutian region, which is also the region of the strongest teleconnection from the Tropical Pacific. This ensemble spread indicates the uncertainty in the observed changes to the ENSO teleconnections due to internal variability and differences in model physics.

Summary and Conclusion
In this study, we have provided a first look at the projected change of ENSO in four CMIP6 future scenarios. Our analysis focused on understanding to which degree the various models agree about projected changes. As reported for the previous intercomparison projects CMIP3 and CMIP5, ENSO is characterized by a high degree of variability and diversity (e.g., Collins et al., 2010;Guilyardi et al., 2012;Taschetto et al., 2014;Yeh et al., 2012) across models and long data records are needed to establish statistically significant changes in its characteristics. While there continues to be no across-model consensus on the change in variance and spectra of ENSO, we see agreement on some emerging signals: 1. In all 11 models the east-west gradient of SST decreases in the future, with larger decreases in the scenarios with higher radiative forcing. A weaker gradient has been associated with increased likelihood of strong East Pacific warm events, which have large socioeconomic impacts. 2. Out of the 11 models 10 show a significant increase in variance of SST in the Niño 3.4 region for at least one SSP and four models for all future scenarios. This increase in variance is likely linked to the decrease in the zonal temperature gradient and increase of strong warm events. 3. While all CMIP6 models are able to produce quasi-oscillatory behavior reminiscent of ENSO, there is a wide range of variability with spectral peaks in the 2-7 year range. Seven out of the 11 CMIP6 models show a significant increase in power spectral density in the ENSO band with periods ranging from 3-7 years. 4. In eight out of the 11 models we see a significant increase in the standard deviation of the E-index for at least one SSP. Previous studies (Cai et al., 2015a;Cai et al., 2018) have linked the change in variability of the E-index to model's nonlinearities, a relationship that does not seem to be as robust for the models studied here. 5. In nine of the 11 models, the centers of the extratropical teleconnection pattern shift eastward and poleward for warm events. However, since the centers of the teleconnections coincide with the regions of largest internal variability, it is hard to establish significance for this shift.
The 11 CMIP6 models analyzed here appear to be in better agreement than the models contributing to the previous intercomparison projects CMIP3 and CMIP5. However, their projections still differ in many key 10.1029/2020GL090640

Geophysical Research Letters
aspects of ENSO, such as the spectra, the representation of ENSO diversity and the change in extratropical teleconnection patterns.
No attempt has been made here to evaluate the models' skill in representing observed ENSO variability. A careful assessment of the models' fidelity in representing ENSO during the historical period together with in-depth process-level analysis might enable to further constrain the projected change in ENSO in current and future CMIP simulations.