The future of the El Niño–Southern Oscillation: using large ensembles to illuminate time-varying responses and inter-model differences

Future changes in the El Niño–Southern Oscillation (ENSO) are uncertain, both because future projections differ between climate models and because the large internal variability of ENSO clouds the diagnosis of forced changes in observations and individual climate model simulations. By leveraging 14 single model initial-condition large ensembles (SMILEs), we robustly isolate the time-evolving response of ENSO sea surface temperature (SST) variability to anthropogenic forcing from internal variability in each SMILE. We find nonlinear changes in time in many models and considerable inter-model differences in projected changes in ENSO and the mean-state tropical Pacific zonal SST gradient. We demonstrate a linear relationship between the change in ENSO SST variability and the tropical Pacific zonal SST gradient, although forced changes in the tropical Pacific SST gradient often occur later in the 21st century than changes in ENSO SST variability, which can lead to departures from the linear relationship. Single-forcing SMILEs show a potential contribution of anthropogenic forcing (aerosols and greenhouse gases) to historical changes in ENSO SST variability, while the observed historical strengthening of the tropical Pacific SST gradient sits on the edge of the model spread for those models for which single-forcing SMILEs are available. Our results highlight the value of SMILEs for investigating time-dependent forced responses and inter-model differences in ENSO projections. The nonlinear changes in ENSO SST variability found in many models demonstrate the importance of characterizing this time-dependent behavior, as it implies that ENSO impacts may vary dramatically throughout the 21st century.

1 Introduction

Understanding how the El Niño–Southern Oscillation (ENSO) will change under increasing greenhouse gas emissions is critical due to ENSO's widespread impacts, which include changes to floods, droughts, and fishery production (e.g., McPhaden et al., 2006; Cai et al., 2012, 2014, 2021; Taschetto et al., 2020). While previous work demonstrates model agreement on future intensification of ENSO's atmospheric impacts (e.g., precipitation; Yun et al., 2021; Power et al., 2013; Fasullo et al., 2018), debate remains as to whether ENSO-driven sea surface temperature (SST) variability will increase or decrease in the future as disparities between projections from different climate models persist (e.g., Cai et al., 2022; Beobide-Arsuaga et al., 2021; Stevenson et al., 2021). In this study we present ENSO SST projections for the first time from 14 individual single model initial-condition large ensembles (SMILEs) with the time-evolving forced response cleanly separated from internal variability.

Previous work finds a diverse range of projections of ENSO variability in multi-model ensembles (CMIP; coupled model intercomparison projects). In CMIP3 a range of ENSO variability projections were found, leading to a review stating that it was not yet possible to determine whether ENSO would change in the future (Collins et al., 2010). In CMIP5 the models again demonstrate a variety of responses with the multi-model mean change indistinguishable from zero (Stevenson, 2012; Beobide-Arsuaga et al., 2021). Other studies have further investigated ENSO projections by excluding models that have strong ENSO biases in the historical period. When selecting for models that correctly represent ENSO skewness, Cai et al. (2018) find an increase in eastern Pacific ENSO SST variability. This result is supported by the recent release of CMIP6 wherein Cai et al. (2022) demonstrate that 34 of 43 models show an increase in ENSO SST variability in SSP585 (27 of 39 in SSP126), in agreement with another early CMIP6 assessment (Fredriksen et al., 2020).

There is also some evidence for an increase in ENSO variability in the observational record, with paleoclimate data suggesting that ENSO variability strengthened post-1950 (Cai et al., 2021). Coral records indicate that ENSO variability has been 25 % larger over the last 5 decades compared to the pre-industrial era (Grothe et al., 2020; Cobb et al., 2013). This agrees with tree ring records, which show a strengthening of ENSO variance in the late 20th century (Li et al., 2013). A study using linear inverse modeling based on instrumental data agrees that ENSO SST variance has increased since 1976 (Capotondi and Sardeshmukh, 2017). The authors, however, note that such differences could be due to multi-decadal variability or the impact of climate change, and the cause of this observed change remains unclear. A recent study also noted the strong contribution of internal variability in the observed record (Gan et al., 2023). Gan et al. (2023) conclude that 65 % of the observed increase in extreme El Niño events is attributable to the internal variability (largely multi-decadal variability from the Atlantic Multidecadal Oscillation), with the rest of the increase attributed to a changing climate.

While consensus between multi-model ensembles and paleoclimate records suggests that ENSO variability may increase with increasing greenhouse gas emissions, other studies that look at the equilibration to strong greenhouse gas forcing complicate this result. A recent study using a very-high-resolution climate model (0.1^∘ ocean; 0.25^∘ atmosphere) finds that ENSO SST variability decreases substantially under strong (4 × CO₂) forcing (Wengel et al., 2021). This study further argues that CMIP-class models are too low in horizontal resolution to correctly capture some important aspects of ENSO ocean dynamics, which could lead to an artificial increase in SST variance in these models. On longer timescales, ENSO amplitude is also found to decrease in CMIP-class standard-resolution (1 to 2^∘) models (Callahan et al., 2021). Using output from the Long Run Model Intercomparison Project, Callahan et al. (2021) find diverging transient ENSO variability projections but a consistent decrease in SST variability across models when the response has stabilized. Indeed, Kim et al. (2014) suggest that ENSO variability in CMIP5 models is not linear in time, implying an important role for analyzing multiple time periods. Collectively, these studies demonstrate the importance of further research to reconcile these discrepancies.

Differences between the projections of ENSO variability in individual models have been linked to the pattern of mean-state warming across the tropical Pacific (e.g., Jin and Neelin, 1993; Fedorov and Philander, 2000; Jin et al., 2006; DiNezio et al., 2012). For instance, in CMIP5 and CMIP6 models there is a weak negative correlation between the magnitude of the tropical Pacific SST gradient and ENSO SST variability (Beobide-Arsuaga et al., 2021). In this case, an increase in ENSO SST variability is linked to a weakening of the tropical Pacific SST gradient. This relationship varies between models (Fredriksen et al., 2020) and is found to be more robust in models that have more realistic subsurface nonlinear dynamical heating and ENSO asymmetry when compared to observations (Hayashi et al., 2020). Periods of weaker tropical Pacific SST gradient are also linked to higher ENSO variability, with periods with a stronger gradient linked to periods of smaller ENSO SST variability (Capotondi and Sardeshmukh, 2017; Rodgers et al., 2004; Ogata et al., 2013; Choi et al., 2013; McPhaden et al., 2011). Models that show reduced ENSO variability in response to anthropogenic forcing also show a strengthening of the tropical Pacific SST gradient (Kohyama and Hartmann, 2017; Kohyama et al., 2017). This relationship between ENSO amplitude and the mean-state tropical Pacific SST gradient is identified in Paleoclimate Modelling Intercomparison Project (PMIP) models, although it is not constant through time (Wyman et al., 2020). However, paleoclimatic proxies suggest that there is no relationship between mean zonal SST gradient and ENSO variability (Cai et al., 2021), demonstrating that this link is either not found in our observations or, as in the PMIP models, is not constant in time. We emphasize that the east–west Pacific SST gradient is only one of many mean-state metrics that are potentially important to influence (and vice versa be influenced by) ENSO variability. In fact, different aspects of the climate mean state affect ENSO feedbacks in multiple ways, leading to a potentially high sensitivity of the linear ENSO growth rate (and hence SST variability) to these factors (Jin et al., 2006).

How the tropical Pacific SST gradient itself will change in the future is also a point of contention. Most CMIP-class models agree on a projected weakening of the SST gradient (El Niño-like warming; Kociuba and Power, 2015; Meehl and Washington, 1996; Fredriksen et al., 2020; Cai et al., 2021). This, however, must be reconciled with the recent observed increase in the SST gradient (La Niña-like warming; Kosaka and Xie, 2013; England et al., 2014; McGregor et al., 2014). Some studies suggest that the inconsistencies between models and observations can be explained by internal variability and observational uncertainty (Watanabe et al., 2021; Chung et al., 2019; Coats and Karnauskas, 2017), while others argue that the observations are truly outside the model range (Kociuba and Power, 2015; Seager et al., 2022). One hypothesis is that under transient forcing the ocean thermostat mechanism, whereby upwelling waters delay warming in the eastern equatorial Pacific, overwhelms the long-term weakening of the Walker circulation and SST gradient (Clement et al., 1996; Heede and Fedorov, 2021), leading to temporary La Niña-like warming. Another study concludes that a more accurate representation of ENSO nonlinearity leads to La Niña-like warming in at least one realistic climate model (Kohyama and Hartmann, 2017; Kohyama et al., 2017). As such, while most models project an El Niño-like warming, this is not the only possible future response. Alternatively, realistic ENSO nonlinearity does not necessarily lead to a La Niña-like warming pattern. While realistic representations of ENSO nonlinearities are a necessary condition for the simulation of rectified mean-state changes via nonlinear dynamical heating, diverging projected changes in future ENSO SST variability among models with largely realistic ENSO nonlinear dynamic heating can lead to either ENSO-induced damped or amplified eastern tropical Pacific mean-state warming (Hayashi et al., 2020). Additionally, it has been hypothesized that SST mean-state biases in the Atlantic Ocean are responsible for low-biased interbasin connectivity (and hence also low-biased SST variability) between the Atlantic and Pacific (McGregor et al., 2018). This suite of studies highlights the need for further research and continued observations to reduce uncertainties in how the tropical Pacific SST gradient will change in the future.

A primary reason for uncertainty in future changes in ENSO variability and the mean-state tropical Pacific SST gradient is that ENSO itself experiences large multi-decadal variability (Wittenberg, 2009), meaning that long averaging periods or large ensembles are needed to identify robust forced changes (e.g., Stevenson et al., 2010; Milinski et al., 2020; Maher et al., 2018; Deser et al., 2020; Ng et al., 2021). Projections using multi-model ensembles can be difficult to interpret because the spread due to internal variability is comparable to the spread due to inter-model differences (Maher et al., 2018; Ng et al., 2021). With this in mind, large ensembles are needed to quantify the time-varying ENSO response (Stevenson, 2012; Milinski et al., 2020). Without large ensembles, previous studies were forced to rely on long time averages to robustly evaluate changes in ENSO variability (e.g., 2000–2099 compared to 1900–1999).

Figure 1 Time series of Niño3 (blue), Niño3.4 (black), and Niño4 (red) variability for the December, January, and February average (DJF) in each model. Variability is computed as the running 30-year standard deviation for each ensemble member after detrending each member by removing the ensemble mean. The DJF mean is taken prior to computing the variability for this and all subsequent figures. The thick line shows the ensemble mean of the running standard deviation, and the thin lines show the 5 %–95 % range across the ensemble. The bottom right panel shows the running standard deviation for ERSSTv5 observations (linearly detrended). The diamond symbols in each panel show the observed variability for the periods 1951–1980 and 1990–2019. Model name, ensemble size, and scenario used are shown in the title of each panel for all figures.

In this study, we present ENSO projections from 14 SMILEs, for which the time-dependent response of ENSO to external forcing can be isolated from internal variability through ensemble averaging. The unprecedented number of SMILEs used allows a detailed examination of the model dependence of future ENSO projections. In particular, this paper aims to

• isolate the forced response of ENSO in individual models,

• evaluate the time evolution of projected changes in ENSO statistics,

• investigate the time-dependent mean-state response of the tropical Pacific (tropical Pacific SST gradient) to greenhouse gas forcing,

• determine whether there is a relationship between projected changes in ENSO variability and the tropical Pacific SST gradient in climate models, and

• investigate the role of greenhouse gases, aerosol, and natural forcing in historical ENSO variability changes using single-forcing SMILEs.

Ziehn et al. (2020)Kirchmeier-Young et al. (2017)Swart et al. (2019)Kay et al. (2015)Rodgers et al. (2021)Jeffrey et al. (2012)Döscher et al. (2022)Wyser et al. (2021)Sun et al. (2018)Burger et al. (2022)Delworth et al. (2020)Boucher et al. (2020)Tatebe et al. (2019)Hajima et al. (2020)Maher et al. (2019)

Table 1 SMILEs used in this study. CMIP6 models are highlighted in bold font.

2 Datasets

The 14 SMILEs used in this study include both CMIP5- and CMIP6-class models and use one of three external forcing scenarios (Table 1). Models used in this study have between 10 and 100 ensemble members for the historical and future scenarios. For those with single forcing we allow as few as three members due to the limited availability of models with any more than a single ensemble member. Previous work has investigated how many ensemble members are needed to investigate ENSO characteristics and find that this depends on the metric used, the length of time averaging, the level of acceptable error, and the model itself (Maher et al., 2018; Milinski et al., 2020; Lee et al., 2021). Lee et al. (2021) find, for example, that for ENSO amplitude 10–47 members are sufficient depending on which model is used to make the estimate, while this range widens to 1–45 members for ENSO seasonality. Maher et al. (2018) determine that to look at ENSO variance over 30-year periods 10 members are sufficient for a 20 % level of acceptable error, while 30–40 are needed to reduce this to 10 %, with 12–20 sufficient for a 15 % error (Milinski et al., 2020). Fewer members are needed to estimate mean-state ENSO metrics than its variability (Milinski et al., 2020; Lee et al., 2021). In this study we accept models with a minimum of 10 ensemble members to include as many models as possible; however, we note that there are larger errors associated with estimates of ENSO variability for the smaller ensemble sizes. We additionally show the time series of ENSO indices in Fig. 1 as is acts to demonstrate the inherent variability due to different ensemble sizes.

While the purpose of this study is to provide a detailed overview of ENSO projections in all available SMILEs, it is useful to provide some observational context. A detailed comparison of ENSO characteristics in both CMIP models and SMILEs can be found in the following studies for a large range of ENSO metrics, and we do not repeat this here, although we note that in general CMIP6 models outperform CMIP5 in 8 out of 24 metrics and are only degraded in 1 (Planton et al., 2021; Lee et al., 2021) and that ENSO-related SST biases and its teleconnections are also improved (Fasullo et al., 2020). We do, however, compare model-simulated and observed El Niño and La Niña composites (Figs. S1 and S2 in the Supplement). We find diverse errors in ENSO patterns, amplitudes, and transitions. All models overestimate the westward extension of ENSO SST anomalies and either overestimate or underestimate the ENSO amplitude, and some models overestimate the duration of ENSO events (e.g., MIROC6 and MIROC-ES2L).

When considering ENSO SST variability we find that only some models capture the observed variability within the ensemble spread (Figs. 1 and S3), but we note that this is not a reliable predictor for future change. For example, both MPI-GE and EC-EARTH3 compare favorably with the observed magnitude of ENSO SST variability but have different future responses (Fig. 1). MPI-GE shows no change in ENSO SST variability, while EC-EARTH3 has strongly increasing ENSO SST variability. We additionally note that there is large multi-decadal variability in the amplitude of ENSO variability in individual ensemble members. This indicates that we need a SMILE to robustly detect a forced change, but also that an actual forced change in the real world might be masked by this large multi-decadal variability.

Figure 2 Spatial pattern of relative change in DJF SST variance in the (2–7 years) ENSO band between the late 21st century (2070–2099) simulations and the historical (1951–1980) period. The differences in variance have been normalized by the variance of the historical (1951–1980) period (contours). The multi-ensemble mean of the individual model ensemble means can be found in the last panel. Red indicates heightened variability, and cool colors indicate dampened variability. For this and following figures using DJF, 30 years are used in the calculations; e.g., for 1951–1980 we include 30 DJFs starting with December 1950 and DJF variance is computed as the variance of the DJF means. We note that a maximum ensemble size of the first 20 members is used in this analysis, with all members of models with 20 or fewer members used.

3 ENSO projections

In this section we evaluate the response of ENSO variability in each individual SMILE in one of three future scenarios dependent on availability (RCP8.5, SSP370, SSP585). We note that there may be differences between SSP370 and RCP85/SSP585, particularly regarding early 21st century sulfate aerosols, which have potential consequences for the evolution of ENSO (O’Neill et al., 2014; Fasullo and Richter, 2022).

3.1 Amplitude

Projections of ENSO amplitude are nonlinear in time and differ between models (Fig. 1). A large proportion of the models demonstrate an increase in ENSO amplitude as characterized by December–January–February (DJF) SST variability of three ENSO indices (Niño3, Niño3.4 and Niño4; ACCESS-ESM1.5, CESM1-LE, CSIRO-Mk36, EC-EARTH, IPSL-CM6A-LR, MIROC6, MIROC-ES2L). This increase is often nonlinear in time, and there is a large spread across models in the magnitude of the increase. Other models show a decrease (CanESM2) or limited change (MPI-GE) in ENSO amplitude. Finally, five models display nonmonotonic behavior (CESM2-LE, GFDL-CM3, GFDL-ESM2M, CanESM5, GFDL-SPEAR-MED). In these models, ENSO amplitude first increases, then plateaus, and lastly decreases. These results are consistent for all three ENSO indices, with subtle differences for some models (e.g., CSIRO-Mk36 and IPSL-CM6A-LR). Given the time-dependent response of ENSO amplitude exhibited in these time series plots, we consider projections for two future time periods (2021–2050 and 2070–2099) in the following sections.

3.2 Pattern

Forced changes in the spatial pattern of ENSO-related SST variability differ between models (Figs. 2 and S4). While the multi-ensemble mean (MEM; used for ease of comparison to CMIP multi-model mean results) shows an increase in variability in the central Pacific, each model has its own unique pattern of change (Fig. 2). Some models demonstrate a general decrease in SST variability along the Equator (CanESM2, CESM2-LE, MPI-GE), while others show a general increase in this region (CESM1-LE, EC-EARTH, MIROC6, MIROC-ESL2L, GFDL-SPEAR-MED). There are also differences in the zonal pattern of variability changes. A group of models has large increases in variability in the central Pacific (ACCESS-ESM1.5, CanESM5, IPSL-CM6A-LR), while another group shows a clear increase across the Pacific (CESM1-LE, EC-EARTH, GFDL-SPEAR-MED, MIROC-ES2L). These diverse patterns are generally similar between the two time periods considered but smaller in magnitude for the earlier period (Fig. S4; 2021–2050) compared to the later period (Fig. 2; 2071–2100), with the MEM pattern consistent between time periods. This result is, however, not true for all models. CESM2-LE has a central equatorial Pacific response of opposite sign between the two periods considered, while GFDL-CM3 and CSIRO-Mk36 have very different patterns of spatial change in the two periods considered. These figures give the pattern of overall change in ENSO SST variability; however, El Niño and La Niña are not symmetric and may evolve in different ways. Due to this we next consider changes in El Niño and La Niña events individually.

Figure 3 Longitude–time sections of the difference of equatorial Pacific (5^∘ S–5^∘ N) SST anomalies (^∘C; shading) composited for El Niño events during 2070–2099 compared to 1951–1980 in 14 SMILEs and the multi-ensemble mean. The SST anomalies for El Niño events during 1951–1980 are overlaid (^∘C; contours at intervals of 0.4 ^∘C; zero contours thickened and negative contours dashed). El Niño events are defined when the Niño3.4 index exceeds 0.75 standard deviations in Dec⁰ but not in Dec⁻¹. Monthly SST anomalies are calculated by removing the ensemble mean from individual ensemble members.

3.3 ENSO event evolution

Composites of El Niño events in each period show increased El Niño SST amplitude relative to the baseline period 1951–1980 in most models and the MEM (Fig. 3; 2070–2099, Fig. S5; 2021–2050). In the baseline period, most models have maximum SST anomalies in the eastern equatorial Pacific sometime between November and February, followed by a westward propagation of warm anomalies and then a switch to cold anomalies in the following year. In the MEM, the longitudinal location of the maximum SST anomaly shifts westward by approximately 15^∘ in 2070–2099 compared to that in 1951–1980, indicating that El Niño SST variability gets stronger over the central Pacific in the future period, consistent with results shown in Fig. 2. In models wherein the amplitude of El Niño increases, the amplitude of the subsequent La Niña also tends to increase (Figs. S6 and S7), presumably due to a stronger heat content discharge, and vice versa for models wherein the amplitude of El Niño decreases. When considering the asymmetry between El Niño and La Niña event composites, the MEM shows that the amplitude of La Niña intensifies more than El Niño in the western-central equatorial Pacific, and the duration of La Niña becomes longer following the amplitude intensification, highlighting the nonlinearity of their responses (Fig. S6). The individual model responses, however, demonstrate strong inter-model differences, which might be related to the diverse model biases in simulating ENSO (Figs. S1 and S2). Overall, ENSO SST changes in the period 2021–2050 show patterns similar, but with weaker amplitude, to those in the periods in 2070–2099 (Figs. S5 and S7). CESM2-LE is an exception that has opposite changes in El Niño SST amplitude and La Niña duration between the two periods.

3.4 Seasonal synchronization

We next investigate the seasonal synchronization of projected changes in ENSO SST variability (Figs. 4 and 5). We aim to determine during which seasons ENSO SST variability changes in each SMILE and determine whether the ENSO seasonality is likely to increase or decrease. We find that models demonstrate different magnitudes of change for different months (Fig. 4). Most models show the largest change (generally an increase) in boreal winter with limited changes in boreal summer. Individual models have opposite-signed changes between the two seasons. Overall, there is an increase in ENSO seasonality (Fig. 5); i.e., increases in ENSO variability are concentrated in the season when ENSO SST variability is climatologically largest. These changes are largely consistent across models, with an overall increase in ENSO SST variability found in boreal fall and winter and a decrease in the variability in spring and summer.

Figure 4 Change in ENSO SST variability, taken as the ensemble standard deviation in the Niño3.4 region for each SMILE and the multi-ensemble mean. Blue lines show the climatology in 1951–1980, shading shows the 5-year running-mean anomalies with respect to 1951–1980, and red lines show the climatology in 2070–2099. Black lines show the observational climatology based on ERSST5 in the period 1951–1980. We note that the climatology is computed as 1970–1999 for EC-Earth in blue and will be updated upon revision.

4 Mean-state projections

The warming rate over the equatorial Pacific varies seasonally and as a function of longitude (Fig. 6). While some models have larger warming trends in the eastern Pacific (El Niño-like warming) for all months (MIROC6, CESM1-LE, EC-EARTH, CESM2-LE), others demonstrate a more complicated seasonal cycle of warming rate, with the location of the largest warming dependent on which month is considered. In the MEM, the eastern Pacific warms more in all seasons; however, from December to June the strong warming trend extends into the central Pacific, which warms with almost equal magnitude to the east. The MEM warming rate is larger over the eastern equatorial Pacific during late boreal spring (the climatological warm season of the equatorial Pacific) compared to boreal autumn (the cold season), but the season of strongest warming is model-dependent. For example, GFDL-CM3 has the largest eastern equatorial Pacific warming in boreal summer, while CESM1-LE has the largest warming in boreal fall and winter.

The time evolution of tropical Pacific SST gradient changes is investigated by considering the SST difference between the eastern and western equatorial Pacific (Fig. 7). A decrease in this gradient indicates stronger warming in the eastern Pacific, i.e., El Niño-like warming (red in Fig. 7), while a strengthening of the gradient corresponds to La Niña-like warming (blue in Fig. 7). We find El Niño-like warming in most models, which is evident in the MEM. This is, however, not the case for CSIRO-Mk36 and is seasonally dependent for GFDL-CM3, GFDL-ESM2M, IPSL-CM5A-LR, and MPI-GE. With the exception of CSIRO-Mk36, models with a stronger climatological gradient tend to have stronger El Niño-like warming (Fig. S3b), although the correlation between the gradient climatology and the gradient change is not significant. The season during which projected changes are strongest is again model-dependent, but SST gradient changes are often concentrated in the season during which the climatological SST gradient is strongest in the specific model considered. The following section relates the projected changes in ENSO variability to the changes in the tropical Pacific SST gradient to investigate whether these changes are linked.

5 Relationship between ENSO variability and mean-state changes

The projected change in ENSO SST variability is linearly related to the projected change in the tropical Pacific SST gradient across ensemble members from all models in the period 2021–2050 (Fig. 8). This relationship breaks down in the later period (2070–2099) for models that have large SST gradient changes (e.g., CESM2-LE, CSIRO-Mk3.6, CanESM2). In general, an ensemble member that experiences El Niño-like warming is more likely to see an increase in ENSO variability, while a La Niña-like warming is associated with the opposite response. While this does not occur for all ensemble members, the larger the response in either the tropical Pacific SST gradient or ENSO SST variability the more likely that a concurrent change in the other variable will be seen.

Figure 5 Change in ENSO SST variability taken as the ensemble standard deviation in the Niño3.4 region for the period 2070–2099 as compared to 1951–1980. The grey lines are individual large ensembles, while the black line is the multi-ensemble mean and the black dashed line is the multi-ensemble median. Red shading indicates an increase in this metric, while blue indicates a decrease.

We further investigate the relationship between changes in ENSO SST variability and the change in the tropical Pacific SST gradient by plotting the evolving relationship between the two variables over time (Fig. 9). Many models have a hook-shaped ensemble mean (forced) response that corresponds to the following time-dependent response. First, there is an increase in ENSO SST variability concurrent with uniform warming across the Pacific (i.e., no change in tropical Pacific SST gradient) or weak El Niño-like warming. Then there is El Niño-like warming concurrent with a plateau in ENSO SST variability. In some models, this is followed by a decrease in ENSO SST variability as the mean state continues to warm in an El Niño-like fashion. This behavior is seen in seven models (CESM1-LE, GFDL-SPEAR-MED, MIROC6, MIROC-ES2L, CanESM5, CESM2-LE, IPSL-CM6A-LR), with an additional four that demonstrate a portion of the hook-shaped response (EC-EARTH3, ACCESS-ESM1-6, CanESM2, MPI-GE). The three outliers are GFDL-ESM2M, which exhibits an increase then a decrease in both ENSO SST variability and the mean-state tropical Pacific SST gradient, GFDL-CM3, which shows an increase then a decrease in ENSO SST variability without much change in the mean-state tropical Pacific SST gradient, and CSIRO-Mk36, which has strong La Niña-like warming and an increase in ENSO SST variability.

Figure 6 Longitude–time sections of the ensemble mean linear trend of SST over the equatorial Pacific (5^∘ S–5^∘ N) during 2015–2099 normalized by the global mean SST trend. The bottom right panel shows the multi-ensemble mean climatological seasonal cycle of SST in the equatorial Pacific over 2015–2044. Vertical black lines delineate the averaging regions used in Figs. 7–10, with the solid lines indicating the western Pacific region and the dashed lines indicating the eastern Pacific region.

The time evolution of ENSO SST variability and the tropical Pacific SST gradient (Fig. 9) provide context for the results in Fig. 8. The internal variability relationship that is the primary contribution to Fig. 8a clearly shows the role of rectification in the mean state (e.g., Hayashi et al., 2020), but for the forced changes this is only one of several mechanisms at play, so the forced changes can depart from this linear relationship (Figs. 8b and 9). The observational changes between 1951–1980 and 1990–2019 also show a notable departure from this linear relationship, with a strong La Niña-like mean-state trend together with a weak increase in Niño3.4 variance (Fig. 9). The multi-decadal La Niña-like mean-state change in observations is near the edge of what could result from the internal variability found in models (z score = −1.88), suggesting that this mean state is likely partially forced (see Seager et al., 2022; Wills et al., 2022). Our results highlight the value of using SMILEs to analyze the time-dependent response of the tropical Pacific in climate models. We additionally note that the individual ensemble member spread is much larger for ENSO SST variability projections than the tropical Pacific SST gradient (Fig. 9). This means that more ensemble members are needed to retrieve the forced change in ENSO SST variability under warming (as well as higher-order moments of other Earth system variables; Rodgers et al., 2021), but somewhat fewer are needed to isolate the forced change in the tropical Pacific SST gradient (e.g., Wills et al., 2020; Rodgers et al., 2021).

Figure 7 Tropical Pacific SST gradient climatology and mean-state changes, taken as the ensemble mean difference between the eastern equatorial Pacific (90–150^∘ W, 5^∘ S–5^∘ N) and the western equatorial Pacific (120–180^∘ E, 5^∘ S–5^∘ N) for each large ensemble and the multi-ensemble mean (see Fig. 6 for region outlines). Black lines show the climatology in 1951–1980, shading shows the 5-year running-mean anomalies with respect to 1951–1980, and red lines show the climatology in 2070–2099. Black dashed lines show the observational climatology based on ERSST5 in the period 1951–1980. Red anomalies show El Niño-like changes, and blue anomalies show La Niña-like changes.

6 A single-forcing perspective

Using two SMILEs (CanESM5 and MIROC6) that have single-forcing simulations available, we separate the contribution of aerosols, greenhouse gases, and natural forcing to historical changes in both ENSO SST variability and the tropical Pacific SST gradient (Fig. 10). Both models have forced increases in ENSO variability over this time period (1985–2014 compared to 1951–1980), although this is not evident in all ensemble members. This increase is consistent with observations and is driven by anthropogenic forcing (greenhouse gases and aerosols), not natural forcing. The larger ensemble (MIROC6) has strong variability across ensemble members in the natural forcing ensemble, with large increases or decreases observed in a single realization (or observations) due to this variability alone. This indicates that ENSO SST variability changes have not yet emerged from internal variability, consistent with previous work (e.g., Ying et al., 2022).

Over the historical period, the tropical Pacific SST gradient weakens in MIROC6 and does not change in CanESM5 in response to all forcings. This is a result of the combination of greenhouse gas and aerosol forcings, which generally oppose each other (Fig. 10). Specifically, greenhouse gases act to weaken the gradient in both models, while aerosols strengthen the gradient in CanESM2 and have a minimal effect in MIROC6. The observed change in the tropical Pacific SST gradient sits right at the edge of the ensemble spread in the all-forcing historical scenario. This suggests that the observed gradient could be due to internal variability; however, it would be an extreme event. Differences between these models and observations could also be due to an incorrect model response to external forcing.

Figure 8 Time differences of the Niño3.4 standard deviation (y axis) and the difference in SST between the eastern and western equatorial Pacific (defined as in Fig. 7; x axis) for each ensemble member of the 14 SMILEs, indicated by colored symbols. (a) Difference between 2021–2050 and 1951–1980; (b) difference between 2071–2099 and 1951–1980. Niño3.4 time series are detrended by subtraction of the time-varying ensemble mean prior to computation of the standard deviation.

7 Discussion

While the MEM change in ENSO SST variability in this study compares well with previous literature, the SMILEs add information by allowing for the isolation of the time-dependent ENSO response in individual models. In most models, although not all, ENSO amplitude is projected to increase, in agreement with Cai et al. (2018, 2022) and Fredriksen et al. (2020). However, while previous work used the differences between long time periods (i.e., 2000–2099 compared to 1900–1999) to isolate the projected change, we leveraged SMILEs to separate the time-dependent response of ENSO amplitude to external forcing from internal variability. Individual SMILEs have a range of projected evolutions of ENSO amplitude that are not linear in time. Additionally, in CESM2-LE, GFDL-CM3, and GFDL-ESM2M ENSO amplitude first increases, then decreases. In these models, the response captured using long time averaging is dependent on the specific period used. Here, SMILEs add value as an important new tool to fully capture the time-dependent ENSO variability projections.

The spatial pattern of projected changes in ENSO variability also differs between models. In the MEM, we find a general shift toward more variability in the central Pacific. However, the individual model patterns that make up the MEM are diverse, and more work is needed to understand the model differences that we find. The MEM increase in the central Pacific variability is also clear in the MEM composite of El Niño event evolution. Additionally, La Niña events appear to strengthen in the MEM more than El Niño events, highlighting nonlinear behavior. In general, for both pattern and event evolution changes, the individual model responses are similar between the two time periods considered (2021–2050 and 2070–2099), albeit smaller in magnitude for the earlier period. However, this is not true for all models, with CESM2-LE a clear outlier showing opposite changes in the earlier and later periods. Model differences could be related to different climatological biases, different patterns of transient mean-state warming, and different ENSO feedbacks and dynamics (e.g., Bellenger et al., 2014; Capotondi et al., 2015; Planton et al., 2021; Wills et al., 2022), with additional work needed to diagnose why individual models behave the way they do.

Figure 9 Time evolution of the 30-year running-mean DJF tropical Pacific SST gradient (defined as in Fig. 7; x axis) and 30-year running standard deviation of DJF Niño3.4 variance (y axis) in each SMILE. Symbols with black edges show the ensemble mean for overlapping 30-year periods every 15 years between 1850 and 2100, with the center year of the averaging period computed as indicated by the fill color. The 30-year running-mean DJF tropical Pacific SST gradient and 30-year running standard deviation of DJF Niño3.4 in each ensemble member is shown for the first and last averaging periods (symbols without black edges). The observational values for ERSST5 during the period 1951–1980 and 1990–2019 are shown with light blue and light green stars, respectively. The internal variability in each quantity expected in a single 30-year period, estimated from the multi-model mean of the ensemble mean variance in each quantity, is depicted with a 95 % confidence interval behind the light blue star in the bottom right panel.

We find that the seasonal synchronization of ENSO SST variability increases in most models investigated in this study. This manifests in an increase in ENSO SST variability in boreal fall and winter and a decrease in spring and early summer. These changes have implications for understanding ENSO's intrinsic dynamics (e.g., Stuecker et al., 2013; Stein et al., 2014; Chen and Jin, 2022), changes to remote ENSO impacts (e.g., potentially amplifying the ENSO impact in one season for a given region and reducing it in a different season), and potential changes to ENSO predictability (e.g., Balmaseda et al., 1995; Timmermann et al., 2018). Thus, changes in variability should be assessed seasonally and not by using annual averages. We note that determining the dynamical cause for the increased ENSO seasonal synchronization in most of the models will require a detailed ENSO feedback analysis (e.g., Chen and Jin, 2022), including assessing potential future changes in the southward wind shift mechanism (e.g., McGregor et al., 2012; Stuecker et al., 2013).

The tropical Pacific SST gradient response also agrees with previous work showing an El Niño-like warming in most models and the MEM (e.g., Kociuba and Power, 2015; Fredriksen et al., 2020; Lian et al., 2018; Cai et al., 2021). However, again the individual model responses that go into the MEM are diverse, with the warming varying both seasonally and longitudinally. In general, the season of strongest SST gradient weakens the most; however, this is again not consistent across all models. In most models, weakening of the SST gradient intensifies in the second half of the 20th century, consistent with a weakening of the Walker circulation gradually overwhelming the ocean thermostat mechanism (Heede et al., 2020; Heede and Fedorov, 2021).

Similar to the previous hypotheses, we find a link between the projected tropical Pacific SST gradient and the change in ENSO SST variability (DiNezio et al., 2012; Beobide-Arsuaga et al., 2021; Fredriksen et al., 2020; Hayashi et al., 2020; Wyman et al., 2020; Capotondi and Sardeshmukh, 2017; Rodgers et al., 2004; Ogata et al., 2013; Choi et al., 2013). However, the SMILEs provide a much clearer picture of the time evolution of anomalies. We find a large subset of models that have an increase in ENSO SST variability that precedes the mean-state change, with ENSO SST variability plateauing, and in some models decreasing as the mean state warms in an El Niño-like fashion. This result begins to conceptually link together previous work that found an increase in ENSO SST variability over the coming century (e.g., Cai et al., 2022) and a decrease in the longer-term equilibrated state (e.g., Callahan et al., 2021). The individual model responses may also put into context a recent high-resolution study (Wengel et al., 2021) that finds a decrease in ENSO SST variability and argues that CMIP-class models cannot capture all aspects of ENSO ocean dynamics correctly, leading to a potentially incorrect projected increase in ENSO SST variability. Based on our study, their model could be capturing a longer-term decrease in ENSO SST variability due to the strong forcing used (4 × CO₂). Alternatively, we do find CMIP-class SMILEs that show a decrease in ENSO SST variability, highlighting the need to use multiple models for robustness. Based on these results, high-resolution modeling studies that use multiple models are needed to reconcile projections and determine the robustness of the higher-resolution result compared to standard-resolution CMIP-class models.

Figure 10 Change in the Niño 3.4 index standard deviation and tropical Pacific SST gradient in DJF in the two available single-forcing ensembles. (a) Change in the Niño3.4 index standard deviation over the second half of the 20th century (1985–2014 minus 1951–1980). Individual ensemble members from CanESM5 (squares), MIROC6 (triangles), and the combination of the two ensembles (circles) are shown with open symbols. The average change is shown with a filled marker. Observed changes (ERSSTv5) are shown by the dashed black line. (b) Change in the tropical Pacific SST gradient (the difference in SST between the eastern and western equatorial Pacific; defined as in Fig. 7; x axis) over the same time period. We use the same visual conventions. Note that the sum of individual forcings may not add up to the “all forcings” scenario, particularly in relatively small ensembles.

In the two models that have single-forcing SMILE experiments, greenhouse gas and aerosol forcings contribute to a historical increase in ENSO SST variability. While one might expect aerosol forcing to have the opposite effect of greenhouse gas increases, aerosols are hemispherically asymmetric unlike greenhouse gas forcing, leading to shifts in the Intertropical Convergence Zone that can also influence ENSO (Kang et al., 2020; Pausata et al., 2020). In contrast, greenhouse gas and aerosol forcings have competing influences on historical trends in the tropical Pacific SST gradient. The observed change in the tropical Pacific SST gradient sits right at the edge of all ensemble members of the all-forcing historical scenario (consistent with McGregor et al., 2014; Seager et al., 2019, 2022; Wills et al., 2022). This suggests that the modeled SST gradient response to greenhouse gas forcing could be too strong or incorrect, the modeled SST gradient response to aerosol forcing could be too weak, the observed change could be an extreme event caused by internal variability, the modeled internal variability could be too small, or some combination of the above. Further investigation is needed as more single-forcing ensembles are released to determine the relative roles of these possibilities.

8 Summary and conclusions

In this study we highlight the value of SMILEs for investigating ENSO projections. SMILEs allow the isolation of the time-dependent response of ENSO SST variability to increasing greenhouse gases. By using SMILEs we can better understand projected changes in ENSO variability and the tropical Pacific SST gradient and isolate inter-model differences in their response to anthropogenic forcing.

Our results show the following.

1. ENSO SST projections

   • a.

     Projections of ENSO amplitude are not linear in time.

   • b.

     Models differ in their projections of the pattern and time evolution of ENSO SST variability and El Niño and La Niña event evolution.

   • c.

     The MEM projects an increase in ENSO SST variability, with more variability particularly in the central Pacific.

   • d.

     The seasonality of ENSO SST variability increases. ENSO SST variability increases in boreal fall and winter with a decrease in spring and early summer. This is qualitatively consistent across all models.

2. Mean-state projections

   • a.

     Most models project El Niño-like warming, although some models project the opposite.

   • b.

     Individual models are different in their longitude and season of maximum warming in the tropical Pacific.

3. Relationship between ENSO SST variability and the tropical Pacific SST gradient

   • a.

     When considering individual ensemble members, tropical Pacific SST gradient changes are linearly related to changes in ENSO SST variability in most models.

   • b.

     This response is time-dependent. In many, but not all, models ENSO variability first increases, and then the tropical Pacific SST gradient weakens as ENSO variability plateaus or decreases.

4. A single-forcing perspective

   • a.

     Increases in ENSO SST variability in the two single-forced SMILEs result from anthropogenic (aerosol and greenhouse gas) forcing.

   • b.

     More single-forced SMILEs are needed to understand the tropical Pacific SST gradient change.

These results present an extensive picture of future changes in ENSO in 14 SMILEs. The use of SMILEs means that the diverse responses across models can be truly attributed to model differences, rather than including contributions from internal variability. These diverse responses demonstrate a need for further investigation into the processes causing model differences in ENSO projections and provide a baseline for future research. Additionally, our results highlight time-dependent behavior including nonlinear changes in ENSO SST variability and changes in the tropical Pacific SST gradient that intensify near the end of the 21st century. This has important implications for ENSO's teleconnections and impacts, as a nonlinear change in ENSO SST variability likely has nonlinear time-dependent changes in its impacts as well. Our results may also help to reconcile previous work that suggests a transient increase in ENSO SST variability, but a long-term equilibrated decrease, by isolating the time-dependent behavior of ENSO SST variability projections. While many models agree on the trajectory of projected changes, not all models behave the same way. This highlights the need for further research on the mechanisms of inter-model differences in ENSO projections. There is a rich diversity of future ENSO changes projected by climate models, and more work is needed to understand which aspects of these projections are robust.