ENSO diversity shows robust decadal variations that must be captured for accurate future projections

El Niño-Southern Oscillation (ENSO) shows a large diversity of events that is modulated by climate variability and change. The representation of this diversity in climate models limits our ability to predict their impact on ecosystems and human livelihood. Here, we use multiple observational datasets to provide a probabilistic description of historical variations in event location and intensity, and to benchmark models, before examining future system trajectories. We find robust decadal variations in event intensities and locations in century-long observational datasets, which are associated with perturbations in equatorial wind-stress and thermocline depth, as well as extra-tropical anomalies in the North and South Pacific. Some climate models are capable of simulating such decadal variability in ENSO diversity, and the associated large-scale patterns. Projections of ENSO diversity in future climate change scenarios strongly depend on the magnitude of decadal variations, and the ability of climate models to reproduce them realistically over the 21st century. El Niño-Southern Oscillation event intensities and locations show pronounced decadal variations, which need to be represented in models for better projections of future ENSO diversity, suggest analyses of observations and climate model simulations.

E l Niño-Southern Oscillation (ENSO) is the leading mode of tropical climate variability, with impacts on ecosystems, agriculture, freshwater supplies and hydropower production spanning much of the globe [1][2][3] . The majority of impact studies, including seasonal to multi-year predictions, has developed from a canonical representation of ENSO, as characterised by seasurface-temperature anomalies (SSTa) in the central-eastern Pacific [4][5][6] . However, ENSO shows large differences from one event to another in terms of its intensity, spatial pattern, and temporal evolution [7][8][9] . For instance, while the 1997/98 El Niño displayed extreme SSTa in the eastern equatorial Pacific (EP-ENSO), the largest SSTa that coincided with the 2002/03 event was weaker and primarily confined to the central equatorial Pacific (CP-ENSO). Differences in the longitudinal location and intensity of ENSO events are sensitively associated with different impacts on regional climate throughout the world 10,11 . Such differences in ENSO patterns, referred to as "ENSO diversity" 7 , and their representation in climate models thus strongly influence the skill of impact-prediction systems 12 , and underscore the need for an appropriate characterisation and further mechanistic understanding of ENSO diversity, as well as its projected changes.
The post-1990s increase in the frequency of CP El-Niño events [13][14][15] has led some researchers to relate such changes in ENSO patterns to the influence of global warming 16 . However, decadal variations in ENSO characteristics have also emerged in unforced climate-model simulations 17,18 , as well as in empirical models 19 , suggesting that natural low-frequency variability may also play an important role in the modulation of ENSO diversity. In particular, the relative frequency of EP-or CP-ENSO has been associated with different phases of the Pacific Decadal Oscillation (PDO) 20,21 , which may itself result from the superposition of different processes 21 , and the Atlantic Multidecadal Oscillation 22 . CP events have also been linked to decadal variations of the North Pacific Gyre Oscillation 23,24 , a result consistent with the presence of spectral power in the decadal range for CP events 25,26 . This decadal ENSO modulation may arise from changes in the system dynamics 27 , from natural and/or anthropogenic external forcing 28 , or occur purely by chance 29,30 . Whatever their nature is, decadal variations in ENSO diversity can be expected to obscure the detection of the climate-change signal, and influence our assessment of 21 st -century projections of ENSO diversity, highlighting the importance of considering the system trajectory, and not only differences between specific future and historical periods, for those assessments.
The detection of low-frequency variations in ENSO diversity has however been largely limited by various technical shortcomings. Most observational studies use relatively short records (≤50 years), which are arguably of insufficient length to study ENSO diversity, and to benchmark models [7][8][9][31][32][33] . Long-term observational SST datasets are reconstructed based on the more densely sampled recent decades, and may be biased toward the patterns of recent events 34 . Century-long ocean reanalysis are typically forced by atmospheric reanalyses, which in turn are constrained by a limited set of observations that become increasingly sparse in earlier periods, thus introducing large uncertainties in the long-term characterisation of ENSO diversity. Climate models can provide long records of ENSO, but their skill in simulating ENSO diversity is generally limited 7,9,26,33,[35][36][37] .
Furthermore, aforementioned studies are traditionally based on several indices of tropical Pacific SSTa, designed to describe the variability of either a canonical ENSO or two extreme ENSO flavours 32,38-41 . These indices have been primarily used to classify El Niño events, since La Niña events, which are typically associated with weaker anomalies than El Niño events 42 , tend to be located further West 7 , and show a more limited range of pattern diversity 43 . Hence, some of these indices may neglect the existing asymmetry between warm-and cold phases, and between CPand EP events. Most of these indices are significantly intercorrelated (Fig. 1a), highlighting that they provide redundant information on ENSO. More importantly, the probability distribution of the peak location of SSTa over the tropical Pacific (when each index is exceeding ±0.5°C standard deviation; Fig. 1b), indicates that most of these indices are neither representing solely CP-and EP events, but rather different combinations of both ENSO flavours.
Hence, it is essential to adopt an approach that allows for a more flexible and precise assessment of changes in the location and intensity of warm and cold ENSO events, by considering the varying longitude of the largest SSTa at any given time, for both El Niño and La Niña events (cf. Methods). This approach is not Fig. 1 Relationships between ENSO indices, and their ability to disentangle CP and EP. a Pearson's correlations between 13 ENSO indices (Niño boxes; 38 PC-based EP-and CP-ENSO; TNI; 40  new. For example, the "Center of Heat Index" (CHI) provides the basis for a more flexible framework, allowing the longitudinal centre of SSTa to vary as a function of time instead of being geographically fixed 34,44 . A similar concept was recently used to provide a single index capturing the longitudinal shifts of ENSOrelated atmospheric convection 45 , and to identify the emergence of double-peaked El Niño events 46 .
In addition, the classification of EP-and CP events is dependent on the specific dataset used (Fig. 1b), suggesting that a probabilistic definition that accounts for the level of agreement among diverse datasets is needed to more robustly assess changes in ENSO diversity. Here, we revisit this classification approach for both El Niño and La Niña events across multiple observational datasets, including century-long reconstructions, ocean and coupled reanalyses, and high-resolution satellite-derived estimates (cf. Supplementary Table 1) to unravel observational uncertainties, while achieving a probabilistic description of the temporal evolution of event location and intensity. The same approach is then applied to the two latest multimodel ensembles, i.e., CMIP5 47 and CMIP6 48 (cf. Supplementary Table 1), to evaluate model fidelity in simulating ENSO diversity within the context of a more flexible framework, which (i) allows for a continuum of patterns, rather than relying on a more rigid bimodal view of ENSO; and (ii) examines models' performance against multiple observational references. Historical and multicentennial preindustrial control simulations (piControl: 400-1200 years, radiative forcing fixed at the 1850 conditions) are used to examine the evolution of location and intensity of ENSO events in the absence of anthropogenic influences (piControl) or under historical radiative forcing (historical), and assess model performance based on multiple metrics (cf. Methods). A subset, including highest-and lowest-performing models, is then selected to examine how events' locations and intensity evolve in the future, and whether significant changes in those quantities can be assessed in the presence of low-frequency modulations.
This study aims at addressing the following questions: (i) can we robustly detect decadal variations in ENSO diversity across a variety of observational/reanalysis products?; (ii) do models also exhibit such low-frequency variations, and how do they impact future scenarios?; do both El Niño and La Niña events undergo similar evolutions? The answers to these questions will provide a novel perspective to the assessment of ENSO 21 st century projections by accounting for system trajectories that do not necessarily evolve monotonically with climate change.

Results
Observed changes in the likelihood of ENSO location and intensity. The probability distribution of ENSO location and intensity is examined in six observational datasets, using 20-year running periods between 1850 and 2018 (Fig. 2). The results obtained using 10-or 30-year running windows are very similar, but the 20-year window is found to optimise the signal-to-noise ratio ( Supplementary Fig. 1).
On average, El Niño events occur over multiple longitudinal locations (Fig. 2a), but show stronger probability to occur either over the central Pacific (~162-171°W; CP-Niño) or the eastern Pacific (~89-132°W; EP-Niño). We note larger observational uncertainty in the probability of EP-Niño (Fig. 2a), which is substantially less likely in ERSST.v5, CERA20C, and especially in OISST.v2, which covers a period dominated by CP events. Comparing the statistical distribution of all observational datasets using 20-year running windows, larger uncertainties in the location of ENSO events are found in the early records (Fig. 2b, grey and black histograms). Nevertheless, looking at the agreement of high probability and the most likely location across all datasets, we identify a coherent and progressive shift to predominant CP-Niño events over the course of the 20 th century (Fig. 2b). This result confirms findings from previous studies, using different approaches, on the recent increase in CP-El Niño frequency [13][14][15] , and on the uniqueness of El Niño patterns 36 . La Niña events also appear characterised by multiple preferential locations, and show a larger degree of consistency in the probability distribution across the different datasets (Fig. 2c). CP-Niña events are systematically more likely than EP-Niña events, especially in COBESST.v2 and CERA20C (Fig. 2c). Looking at temporal changes in the statistical distribution across all datasets, coherent low-frequency variations emerge in the most likely locations of La Niña (Fig. 2d). CP-Niña events are more frequent in the 1930-40s and after the 1970s than during the 1950-60s (Fig. 2d).
All observational datasets show quasi-normal distributions for event intensity, with an average most probable intensity of +0.84°C and −0.92°C for El Niño and La Niña events, respectively (Fig. 2e, g). In all observational datasets, the probability distribution of La Niña intensity shows relatively small variations over time using a 20-year running window. These variations, however, exceed the range of variability that would be obtained in a random field (Fig. 2c, dashed lines and grey shading), and are substantially larger using a 10-year window ( Supplementary Fig. 1j). Decadal variations in intensity are stronger for El Niño events (Fig. 2e). The most likely El Niño intensities display relative minima in the late 19 th century, the 1920-50s, and 1980-90s, as compared with the early 20 th century, the 1960-70s and post-2000s (Fig. 2f). The reduced ENSO activity in the 1920-1950s coincides with a period of larger disagreement between observational datasets (Fig. 2f, grey and black histograms), which could explain reduced skill in the prediction of ENSO indices over the same period 49 . In addition, we note that the upper tail of the distribution, as characterised by the frequency of events exceeding the 90 th percentile (Fig. 2f, red histogram), is larger at the beginning of the century and after the 1970s, with a peak in the 1980-90s, which is likely associated with the extreme El Niño events of 1982/83 and 1997/98. This result is consistent with coral oxygen isotopic reconstructions 50 and simulated long-term future changes 51 . Similarly, an increase in the frequency of extreme La Niña events is identified after the 1960s (Fig. 2h, blue histogram), and especially in recent decades, consistently with model projections 52 . Importantly, however, unlike the recent increase in the frequency of extreme La Niña events, the intensification of El Niño events does not seem to exceed the historical range of variability (Fig. 2f, h, blue and red histograms).
Hence, consistent with a clustering analysis of equatorial Pacific SSTa 53 , observed diversity of ENSO patterns is large, even for La Niña events, whose diversity was questioned in previous studies 41 . More importantly, ENSO diversity is clearly linked to low-frequency variations, with progressive shifts toward more eastward or westward locations, and more or less intense events, which seem to substantially modulate recent trends in its behaviour. In addition, these low-frequency variations in ENSO location and intensity exceed the range of variations that would occur in a random field during some periods, notably the westward shift in the location of El Niño events (Fig. 2b, dashed lines and grey shading), suggesting that such variations are unlikely to occur purely by chance.
Evaluation of ENSO diversity in climate models. We compare the simulated probability distributions of ENSO location and intensity from 26 CMIP5 and 28 CMIP6 models to our six observational datasets in Fig. 3.
As seen in Fig. 2, El Niño and La Niña events show a primary location centred around 162-171°W, and a secondary location centred around 89-132°W (Fig. 3a, d). Yet, the range of simulated event locations varies significantly across models and model generations (Fig. 3b, c, e, f). Many models are capable of simulating multiple preferential locations (Fig. 3b, c, e, f, grey dots), and between 37 and 50% of simulations do not show significant biases in the mean location of El Niño and La Niña events in historical and piControl runs ( Supplementary Fig. 2). However, estimating biases in the mean location of ENSO events can conceal the presence of asymmetrical and multi-modal probability distribution. Indeed, only two CMIP6 models (IPSL-CM6A-LR and UKESM1-0-LL, one of the MOHC group of models) and one CMIP5 model (CNRM-CM5; leftmost columns of the CNRM group of models) display a broad range of event location that is somewhat comparable to observations (Fig. 3b, c, e, f). Most models show highly asymmetric probability distribution, with clear tendency to favour either EP-or CP events (Fig. 3b, c, e, f), including some extreme cases (e.g., CSIRO-Mk3-6-0) with events always centred further west than the observational range. Such excessive westward-shifted mean location of ENSO could be explained by anomalous westward extensions of the equatorial SSTa in CMIP5 54-56 and CMIP6 models.
Simulated-intensity distributions of ENSO are mostly consistent with observations, and tend to follow a quasi-normal distribution in most CMIP5 and CMIP6 models, but clear discrepancies emerge in the mean intensity and the probability of extreme events (Fig. 3g-l). In approximately 50% of the simulations, biases in the mean intensity of ENSO events are nonsignificant for warm-and cold events, using both historical and piControl simulations ( Supplementary Fig. 2). Large and significant overestimations and underestimations nevertheless persist in a little less than half of the models ( Supplementary  Fig. 2).
In summary, CMIP5 and CMIP6 models simulate some range of pattern diversity for El Niño and La Niña events. However, very few models present a range of event locations in relatively good agreement with observations, and minimal biases in their intensity. Most models hardly depart from an eastward or westward location, favouring either EP-or CP events, and presenting limited ENSO diversity. Notably, larger biases are found for the models that produce erroneous westward extensions of SSTa (cf. MIROC6, Supplementary Figs. 3-6).
Robust decadal variations in ENSO-preferred location and intensity. Here, we use spectral analysis to examine whether robust and significant low-frequency variations are found in the most likely location and intensity of ENSO, in both the five longterm observational datasets, CMIP5 and CMIP6 models. We first examine the timescales on which ENSO behaviour varies using maximum power spectrum (Fig. 4), a method that accounts for nonstationarity of ENSO spectral characteristics (cf. Methods), before comparing the observed and simulated magnitude of decadal variations in ENSO-event location and intensity (Fig. 5).
According to previous observational studies 57 , as well as fossil coral oxygen-isotope records 50 , the observed intensity of ENSO presents significant variability on interdecadal timescales at 16-yr and, especially, 32-yr periods (Fig. 4g, j). CMIP5 and CMIP6 also simulate significant interdecadal variability at these timescales ( Fig. 4h, i, k, l), in agreement with previous studies using climate models 17,18,30 . Many models (40-62% in historical and piControl simulations) display statistically equal decadal variance in the intensity of both El Niño and La Niña (Fig. 5), in agreement with observations. A substantial fraction of simulations (28-42% in historical and piControl runs) significantly overestimate the decadal variance in ENSO intensity (Fig. 5), as identified in the Niño3 + 4 index using CCSM4 26 and CESM2 35 .
Our results confirm the existence of a significant interdecadal modulation in ENSO intensity in accordance with several studies based on observations, proxy records, and climate models 17,18,30,34,50,57 . While previous studies reported an underestimation of decadal variability by climate models at both global 28 and Pacific Ocean 58,59 scales, this statement does not appear to be true for ENSO diversity in many CMIP5 and CMIP6 models, when considering their nonstationary behaviour (cf. Methods). Our results also reveal significant decadal modulations in the maximum likelihood of ENSO locations, which are robust and consistent in both observations and climate models. Nevertheless, many models show recurrent biases in simulating realistic magnitudes of decadal variance in ENSO diversity.
Large-scale patterns linked to decadal variability in ENSO location and intensity. To identify large-scale patterns of variability associated with spatiotemporal variations in ENSO, separately for El Niño and La Niña years, we compute linear regressions of pan-Pacific SSTa, wind stress, and equatorial 20°C isotherm depth (Z20) on the location and intensity of events, using five long-term observational datasets and 32 historical simulations from the IPSL-CM6A-LR large ensemble (Fig. 6). We focus on the Pacific region, as regressed SSTa is much lower, and often nonsignificant, in the other ocean basins. Similarly, we choose to focus on the IPSL-CM6A-LR model because it provides The results obtained using other CMIP5 and CMIP6 models are very similar, especially for patterns associated with ENSO intensity, and approximate the skill of the SODA.si3 reanalysis (Supplementary Figs. 3 and 4). In addition, regressed SSTa patterns associated with the changes in the location of ENSO events provide similar information than the difference between composite of SSTa for EP-and CP events ( Supplementary Fig. 5).
In observations, SST regressions on El Niño longitude yield an EP-type event with the largest anomalies extending westward along the Equator from the coast of South America (Fig. 6a). The associated strong westerly wind anomalies extend to the eastern Pacific, where the thermocline is significantly deeper, while slight, but significant, easterly wind anomalies and shallower thermocline are found in the western Pacific (Fig. 6a, b). These patterns indicate that, during El Niño years, EP-Niño events are more likely to occur when trade-winds weaken (strengthen), and the thermocline is significantly deeper (shallower), over the eastern (western) Pacific; meanwhile, opposite wind and thermocline anomalies seem to favour the occurrence of CP-Niño events. This is consistent with previous studies stressing the importance of the initial zonal thermocline slope as a discriminating factor for the selection of EP and CP events 60 . Compared with regression patterns associated with El Niño longitude, regressions on La Niña longitude show much stronger (weaker) signals in the western-central (eastern) Pacific (Fig. 6a, c). However, this might only denote that, in observation, La Niña events are more (less) likely to occur over the central (eastern) Pacific than El Niño events (Fig. 2a, c). As illustrated in Supplementary Fig. 5a, e, this is also because EP-Niño events tend to extend further west than EP-Niña events. Regressions on La Niña longitude result in a pattern that is reminiscent of a CP-Niño pattern 41 , with cold anomalies in the far eastern Pacific and warm anomalies in the central Pacific (Fig. 6c). In this case, strong westerly wind anomalies and deeper thermocline are found in the central Pacific, where they may contribute to the zonal advective feedback 61 , while weaker easterly anomalies and deeper thermocline are present in the western Pacific (Fig. 6c, d). Such patterns indicate that, during La Niña years, EP-Niña events are more likely to occur when trade-winds strengthen (weaken), and the thermocline is significantly shallower (deeper), over the eastern (western) Pacific; on the contrary, opposite wind and thermocline anomalies would favour the occurrence of CP-Niña.
These tropical signals are statistically significantly related to extra-tropical SSTa (Fig. 6a, c): colder (warmer) North Pacific SSTa is found, when El Niño events are located further east (west), and La Niña events are located further west (east). While these results corroborate previous study on changes in the frequency of CP and EP events during different phases of the PDO 21,62 , such North Pacific SSTa is also consistent with the Pacific Meridional Mode 63 , the Victoria mode [64][65][66] , and with changes in the intensity and location of the Aleutian Low and North Pacific High in response to EP-and CP-Niño 67 . Similar regression patterns are found in IPSL-CM6A-LR, and in other models ( Supplementary Figs. 3 and 4), which can produce . d-f Same as (a-c) but for La Niña. g-l Same as (a-f) but for El Niño and La Niña intensity. Only significant variability patches at p = 0.05, as determined using 1000 Monte-Carlo simulations of the red noise-background spectrum, are shown. Dashed red lines and grey shading indicate the area where variability can be underestimated because of edge effects, wraparound effects, and zero padding. As the continuous-wavelet analysis allows to account for temporal changes, the maximum power spectra are estimated using the full length of each time series.
realistic changes in zonal wind stress and thermocline depth, associated with shifts in ENSO locations (Fig. 6e-j). Patterns associated with changes in El Niño and La Niña locations are, however, much more symmetric in models than in observation ( Fig. 6e-j; Supplementary Figs. 3 and 4). Like other models, IPSL-CM6A-LR shows large internal variability in thermocline-depth anomalies (Fig. 6g, j), with a clear tendency to underestimate thermocline-depth anomalies during El Niño events (Supplementary Fig. 6), and this could explain larger ensemble spread in equatorial Pacific SSTa associated with shifting ENSO locations (Fig. 6f, i). The North Pacific anomalies associated with ENSO locations are also significant in IPSL-CM6A-LR (Fig. 6e, h), like in many other models (Supplementary Fig. 4). These relationships between ENSO and Pacific extratropical variability however show large ensemble spread in IPSL-CM6A-LR (Fig. 6f, i), highlighting that these relationships are highly sensitive to internal variability, as suggested by previous studies 21,68 .
Looking at regressed patterns associated with event intensity, patterns of SSTa are extending in the central-eastern Pacific, for both El Niño and La Niña (Fig. 6k, m). In addition, we found that observed El Niño (La Niña) events are more intense when the mean thermocline is deeper (shallower) and the trade winds are consistently weaker (stronger) over the equatorial Pacific ( Fig. 6k-n). Compared with the large-scale patterns associated with ENSO locations, changes in ENSO intensity are associated with larger wind stress and thermocline-depth anomalies over the central-eastern equatorial Pacific (Fig. 6k-n). ENSO intensity also appears associated with extra-tropical SST and wind anomalies that are more symmetric about the Equator compared with those associated with the location (Fig. 6k, m), and are somewhat reminiscent of the extra-tropical signature of the Interdecadal Pacific Oscillation 63 (IPO). Other studies discussed the separate importance of North and South Pacific climate variability on ENSO intensity at interannual to decadal timescales 68,69 . Although it systematically underestimates both zonal wind-stress and Z20 anomalies compared with observations ( Supplementary Fig. 7), IPSL-CM6A-LR exhibits large-scale anomalies associated with event intensity that are similar to observations (Fig. 6o-t). Other models also show similar results ( Supplementary Figs. 3 and 4). Most of them simulate coherent changes in wind stress and thermocline-depth anomalies over the equatorial Pacific and extratropical regions that are comparable to observations, during both El Niño and La Niña events. Interestingly, IPSL-CM6A-LR shows very little ensemble spread in equatorial Pacific SSTa, while the strength of extra-tropical anomalies and equatorial thermocline responses strongly differs from one simulation to another (Fig. 6p-t).
Impact of decadal variations on future scenarios for ENSO diversity. We next examine ENSO location and intensity in climate-change projections, using three models (IPSL-CM6-LR, UKESM-1-0-LL, and CNRM-CM5) that produce realistic ranges of event location and intensity (Fig. 3), decadal variability (Figs. 4 and  (Fig. 6), during the historical period (Fig. 7). For comparison, future scenarios of ENSO diversity are also examined in three other models that have a less realistic range of event location, and favour either EP-or CP-ENSO during historical and preindustrial periods (Supplementary Fig. 8).
According to IPSL-CM6-LR and UKESM-0-LL, most ensemble members converge to more CP-ENSO over the second half of the 21 st century (Fig. 7a, b, d, e). This shift to more westward events appears quite early in IPSL-CM6-LR, while it only emerges in the second half of the 21 st century in UKESM1-0-LL, as the ensemble mode is dominated by decadal variations over the first half of the 21 st century (Fig. 7a, b, d, e). Such decadal variations in the ensemble mode remain stronger than potential trends throughout the 21 st century in CNRM-CM5 (Fig. 7c, f). The impact of decadal variations in modulating potential trends in ENSO location is even larger when looking at individual simulations ( Fig. 7a-f). For instance, different simulations of CNRM-CM5 suggest antipodal increase in the probability of CP-or EP-Niño events in the near-and long-term future (2020-2040 and 2060-2090), depending on the phase of decadal variations (Fig. 7a-f). Similar diverse behaviour in simulations of future evolution of ENSO event location is found in CNRM-CM6-1, MIROC6 and CESM2 ( Supplementary Fig. 8a-f). In all cases, future trajectories for ENSO location are strongly influenced by the presence of decadal variations, which may themselves be model dependent and not necessarily consistent with observations. These two factors could explain why a recent study of ENSO uniqueness, using future CMIP5 model scenarios, and comparing the multimodel mean over specific decades, could not detect robust changes in ENSO diversity over the 21 st century 36 . Fig. 6 Large-scale patterns driving long-term variability in ENSO location and intensity. a Regressed SST (blue-to-red shades) and wind-stress (vectors) anomalies associated with changes in El Niño location using multiple observational datasets (SST: ERSST.v5, COBESST.v2, HadSST1, SODA.si3, CERA20C; wind stress: NOAA-20CR.v3; Z20: SODA.v2.2.4). b Same as (a) but for Z20 anomalies. c, d Same as (a, b) but for La Niña events. e-j same as (a-d) but using the IPSL-CM6A-LR large ensemble (32 members). k-t Same as (a-j) but for El Niño and La Niña intensity. While SST and wind-stress anomalies are displayed at the pan-Pacific scale based on the median changes in (a, c, k, m) observations, simulated regressed anomalies are assessed through (e, h, o, r) the ensemble median and (f, i, p, s) ensemble spread (standard deviation [SD]). Z20 anomalies are estimated through the median changes between 5°S and 5°N (b, d, l, n, g, j, q, t). Red and blue shades on the Z20 anomalies indicate the spread between the four SST observational datasets and within the IPSL-CM6A-LR large ensemble (light to dark: maximum/minimum, 10/90 th , 30/70 th , and 45/55 th percentiles), for El Niño and La Niña, respectively. Group 1 (black lines) and Group 2 (grey lines) illustrate how two opposed types of equatorial Z20 anomalies influence the ensemble spread. Statistical significance is assessed at p = 0.05 using 1000 permutations, and displayed as black contour for SSTa, and blue/red crosses for Z20 anomalies. Only significant wind-stress anomalies at p = 0.05 are displayed.
Nevertheless, two of the models with realistic representation of ENSO diversity during the historical period suggest a potential westward shift in the location of both El Niño and La Niña events over the 21 st century. Whether such changes could be attributed to anthropogenic climate change or to internal variability would however require using a larger ensemble.
In our subset of models, event intensity and the frequency of extreme events both show very large decadal variations over the historical period and in the 21 st century (Fig. 7g-l). These decadal variations are as large as any potential trends in CNRM-CM5, UKESM1-0-LL, CNRM-CM6-1, and CESM2 (Fig. 7h, i, k, l; Supplementary Fig. 8g, i, j, l), but these models strongly overestimate decadal variations in event intensity (Fig. 5). In IPSL-CM6-LR (in the second half of the 21 st century; Fig. 7g, j) and MIROC6 (from the early-to mid-20th century; Supplementary Fig. 8h, k), which show weaker biases in the magnitude of decadal variability in event intensity, a clear trend toward more intense events is found. Such an intensification of ENSO events is consistent with previous findings, using CMIP5 models, of an increased frequency of extreme El Niño 51 and La Niña 52 , although the criterion for model selection was different in those studies. In our approach, information on event intensity does not discriminate between EP-and CP events, but since EP events tend to be stronger than CP events 7 , the increasing ENSO intensity seen in some of the models used here is also consistent with previous results indicating enhanced EP-ENSO 70,71 . More importantly, such trends toward more intense ENSO events may remain largely modulated by significant decadal variations (Fig. 7g-l;  Supplementary Fig. 8g-l). Indeed, as already discussed over the historical period (Fig. 4), these decadal variations in ENSO location and intensity exceed the range of variations that would occur in a random field during some periods ( Fig. 7; Supplementary Fig. 8), suggesting that they cannot entirely occur by chance. This is particularly noticeable in individual runs, where these variations are not averaged or flattened ( Fig. 7; Supplementary Fig. 8). Thus, future changes in ENSO characteristics are not necessarily monotonic, as often assumed, but may undergo variations that need to be better understood and accounted for in assessments of future changes. Indeed, due to the system's nonlinearity, these decadal variations can modulate the ENSO response to anthropogenic forcing, as ENSO variations over a given period appear to be related to its own variability over the previous decades 72,73 .

Discussion
To overcome existing limitations in analysing ENSO diversity, this study uses multiple observational datasets to provide a . d-f Same as (a-c) but for La Niña events. g-l Same as (a-f) but for El Niño and La Niña intensity. Grey histograms on the bottom axis of the intensity panels indicate the average number of extreme events (as defined in Fig. 2) within the full-model ensemble. Dashed lines and grey shading indicate the range of variability that could be obtained in a random field at p = 0.05. SSTa is estimated by removing the 1850-2014 monthly climatology and trend, to allow comparison with observations. The same baseline period was used to estimate the 90 th and 10 th percentiles.
probabilistic description of historical variations in event location and intensity, and to benchmark models, before examining future system trajectories. Using multiple century-long observational datasets, CMIP5 and CMIP6 multi-model ensembles, we first identified robust long-term changes and variability in the likelihood of El Niño and La Niña location and intensity. Although the majority of models favour either EP-or CP-ENSO, we found that ENSO diversity is closely linked to significant decadal variations in both observations and climate models. These decadal variations do not only modulate event intensity, as already highlighted in previous studies 17,18,30,57 , but also affect event location.
We identified large-scale variability patterns associated with low-frequency changes in ENSO location and intensity that are similar in both observations and models, despite the large models' underestimation of equatorial zonal wind stress and thermoclinedepth anomalies. On the one hand, long-term changes in event location are associated with zonal perturbations in equatorial wind stress, which, according to previous studies 21, [62][63][64][65][66]68 , may be related to the North Pacific climate variability, and with significant modulations of the thermocline response over the central Pacific. On the other hand, low-frequency changes in event intensity are associated with strong equatorial wind stress and thermocline response, whose variability appears associated with the North and South Pacific climate variability. However, disentangling the complex mechanisms driving changes in the ENSO location and intensity remains a long-standing challenge. This is beyond the scope of this study and warrants further investigation. Although these decadal variations appear to be statistically significant during some periods in the observational record and in the models, they may also arise as residuals of ENSO events 74,75 , and/or be associated with influences from other ocean basins 22 , which could themselves be stochastic in nature. External forcing cannot be ruled out either, especially in the context of climate change 76 .
Our analysis of a small set of climate models, representing key features from the full range of uncertainties in terms of ENSO diversity and its variability, suggests a tendency toward more CP-El Niño and -La Niña over the 21 st century. This result is thus in line with a previous examination of projected changes in ENSO diversity, using CMIP3 models, which suggested an increased frequency of CP-El-Niño 16 . However, our results also highlight that 21 st -century projections of ENSO diversity, in terms of both location and intensity, may depend on the magnitude of the models' decadal variations, and on the models' skill at reproducing realistic ranges of location and intensity. This is consistent with a recent study of ENSO uniqueness in CMIP5 models, which was unable to detect changes in ENSO diversity over the 21 st century in the multimodel mean, and over specific decades 36 . Decadal variations, which are found to be more or less pronounced in different models and in different simulations of the same model, may also influence recently reported trends in the frequency of extreme El Niño 51 and La Niña 52 events over the course of the 21 st century.
This study also illustrates how the future evolution of ENSO diversity may vary significantly from model to model, and stresses the importance of considering the system's trajectory, rather than differences between fixed future and historical periods to assess changes. This is crucial because future ENSO responses to anthropogenic climate change may strongly depend on the magnitude of decadal variations over the historical period, via self-modulation processes 72,73 , which need to be adequately captured in climate models. In addition, the exact evolution of the system is very important in itself, as a monotonic change in ENSO characteristics may result in very different impacts than changes undergoing large decadal variation, as seen in some models. The nature of these decadal variations of ENSO diversity remains unclear. Large model ensembles will be key to separate the natural and forced components of these decadal variations, and will be considered in a future study.

Methods
Observational reference datasets. We use six observational datasets, covering all the different ways to reconstruct long-term variability for SST, as well as different resolutions (Supplementary Table 1 To examine the potential large-scale patterns associated with changes in the ENSO spatiotemporal variability, surface wind stress was derived from surface zonal and meridional winds for the period 1870-2015, using the NOAA-CIRES-DOE Twentieth Century Reanalysis version 3 84 (NOAA-20CR.v3). The NOAA-20CR.v3 uses SODA.si3 and HadSST1 as boundary forcing, and therefore provides consistent atmospheric circulations for those SST datasets. Because subsurface potential temperature data are not currently available in SODA.si3, we use SODA.v2.2.4, with NOAA-20CR.v2 as boundary forcing, to provide the most consistent estimate of thermocline depth, using the 20°C isotherm depth (Z20) as a proxy.
CMIP5/6 simulations. We use 95 ensemble members of historical simulations from 26 CMIP5 models 47 , and 250 members from 28 CMIP6 48 models, together with longer piControl runs (Supplementary Table 1), to evaluate how climate models perform in simulating ENSO diversity. Each individual member of historical simulations allows inferring climate variability from the mid-19 th to the early 21 st century, due to changes in anthropogenic and natural forcings, while accounting for uncertainties associated with internal variability 85 . Similarly, piControl simulations enable assessing the uncertainties associated with the limited length of reliable historical records. In addition, to discuss the implications of our results for future scenarios of ENSO diversity, we use the highest emission scenario or forcing level (8.5 W.m −2 ), i.e., the Representative Concentration Pathway RCP8.5 in CMIP5 models, and the Shared Socio-economical Pathway 5 that updates the highest forcing level, i.e., 8.5 W.m −2 (SSP5-85) in CMIP6 models. The number of available realisations is substantially lower in future scenarios than historical runs (Supplementary Table 1). Monthly fields of SST, zonal and meridional wind-stress and potential temperature (from which we estimated the thermocline depth from Z20) are used. To ensure consistency with the observational datasets, and to optimise the detection of changing locations and intensity in ENSO, model simulations were all interpolated onto a regular 1.25 × 1.25°grid in the ocean and the atmosphere.
Examining long-term variability and changes in ENSO location and intensity. To better account for multidimensionally varying properties of ENSO, building on the CHI concept 34,44 , and other recent studies 37,45,46 , we use a flexible framework, which allows for a continuum of patterns, rather relying on a more rigid bimodal view of ENSO, to estimate the location and intensity of El Niño and La Niña events. The location of ENSO events has been defined as the longitudinal location of the absolute maximum of SSTa, within a strip that spans the tropical Pacific from 150°E to 60°W (excluding the warm-pool region), and averaged between 5°S and 5°N over the boreal winter months (December-February). To harmonise the results over variable grid resolutions, and reduce the noise in the signal, the location of the absolute maximum of SSTa has been estimated using a 2°longitudinal smoothing (corresponding to the resolution of the coarser observational SST dataset; cf. Supplementary Table 1). The intensity of events is given by the value of the absolute maximum of SSTa at that location and during the same season. While in previous studies, a threshold for SSTa intensity was set to ±0.5°C 34,44 , 1 standard deviation 46 , or to a convective threshold 45 , to qualify as an El Niño or La Niña event, here, we classify all events greater (lower) than 0°C as El Niño (La Niña). This allows for a better consideration of neutral-or weak-ENSO events, and provides smoother interevent transitions, which are crucially important when analysing low-frequency variability. Similarly, El Niño and La Niña events, including subsequent events, are examined separately to better account for nonlinearity between warm-and cold phase of ENSO. SSTa is calculated by removing the mean and trend of each month, before detrending the data using a locally weighted regression (LOESS) 86 , which is less sensitive to outlier than linear regressions and better suited in the presence nonlinear trends. We adopted a probabilistic approach to quantify the likelihood of event location and intensity, and the degree of agreement among diverse observational datasets. This is achieved by estimating the probability-density function (PDF) over every 20-year segments of each observational dataset, and calculating the most likely values (i.e., the mode in statistical terms), as well as multidataset agreements of high probability (i.e., probability exceeding 0.01 and 0.4 for event location and intensity, respectively). Specifically, the agreement of high probability enables us to examine the temporal evolution of the heads of PDFs across multiple datasets, allowing for detection of multiple heads, while limiting the influence of outliers. Here, focus is given to the results of a 20-year running window, as it is found to optimise the signal-to-noise ratio, but very similar results are found using 10-or 30-yr running windows (Supplementary Fig. 1). We also quantify the percentage of disagreement between them, by comparing their respective PDF using a Kolmogorov-Smirnov (KS) test at p = 0.05 and 0.1, using the same 20-year running window. Significance of the variations in the most likely location and intensity of events are then tested against the variations that would be occurring in a random field at p = 0.05. This is achieved by replicating the above procedure using 1000 random fields, generated through autoregressive models of order 1 with the same mean and variance as the observational field. In addition, as the most likely location and intensity might not always be representative of changes affecting the tails of the distribution (or extremes), we also examine the changes in the frequency of extreme El Niño (La Niña) events, by quantifying the 20-year average number of events exceeding (lower than) the 90 th (10 th ) percentile across all datasets.
Then, we further explore the long-term variability using the 10-year most likely location and intensity of El Niño and La Niña events. Continuous-wavelet analyses are used to estimate the maximum power spectrum over the full length of observational and simulated records, while accounting for temporal changes 87 . Using continuous-wavelet analysis enables us to account for nonstationary significant patches of variability, which might not be significant over the full-length of the records, and would not be identified using fast Fourier transform.
Significance of variability patches is tested at p = 0.05, based on 1000 Monte-Carlo simulations of the red noise-background spectrum.
Testing robustness in climate models, identifying large-scale patterns, and implications for future scenarios. We first examine whether historical and piControl runs, from CMIP5 and CMIP6 models, are able to reproduce a realistic range of locations and intensities for both El Niño and La Niña events, by comparing the simulated PDF to multiple observational datasets. This visual comparison of the PDF skewness and kurtosis (or range) is combined with two statistical tests: (i) test for multimodality, i.e., the presence of multiple peaks on the PDF, based on kernel-density estimators and the quantification of excess mass; 88 (ii) test for difference in the mean using a two-sided student t-test (cf. Supplementary  Fig. 1). Statistical significance of these tests is calculated using 1000 permutations.
Second, we investigate whether significant decadal variability is detectable in climate models, by comparing the simulated maximum power spectra with observations. We then compare the simulated magnitude of decadal variability to the observed one using the centred ratio of standard deviation rSD ¼ 1 À sdðENSO int loc ½obs 10yr Þ sdðENSO int loc ½sim 10yr Þ 100 . Statistical significance is then assessed by performing a two-sided Fisher's F-test at p=0.05 between every 100-yr segment through the course of climate simulations, as compared with the period with significant decadal variance in the four longer-term observational SST datasets, from which the rate of success is quantified. Third, we compare the observed large-scale patterns associated with long-term variability in the location and intensity of El Niño and La Niña events to historical simulations. This consists of examining the differences in the patterns of pan-Pacific SSTa, wind stress, and thermocline depth at the Equator (5°S-5°N), which are computed using linear regression during composite El Niño and La Niña years, separately. Statistical significance of the regression patterns is calculated using 1000 permutations. As illustrated in Supplementary Fig. 5, regressed SSTa patterns associated with the changes in the location of ENSO events provide similar information than the difference between composite of SSTa for EP-and CP events. After comparing all the simulated patterns of pan-Pacific SSTa ( Supplementary  Fig. 3), we focus on a set of climate models that represents key features from the full range of uncertainties associated with CMIP5 and CMIP6 models in terms of ENSO diversity. In particular, we focus on the IPSL-CM6A-LR large ensemble model, which displayed similar mean and ranges of ENSO locations and intensity, as well as similar decadal variance to observation. Comparisons with other models presenting similar performance (UKESM1-0-LL), and models showing a tendency to favour EP-or CP events (CNRM-CM6-1, MIROC6), are provided in Supplementary Figs. 4-6.
Finally, using an extended set of climate models, including three highest-and lowest-performing models, we examine future trajectories for ENSO diversity (i.e., location and intensity) using RCP8.5 and SSP5-8.5 emission scenarios, using the same procedure as for the observational datasets. Specifically, we analyse how decadal variations can influence detection of future changes, and how the results differ in two contrasted simulations of the same models that only vary by their initial conditions, and depending on model performance. Once again, highestperforming models include models that represent similar mean and ranges of ENSO locations and intensity, as well as similar decadal variance and climatic background conditions, to observation during the historical period. Meanwhile, lowest performing models are models that present a tendency to favour either EPor CP events, as well as diverse biases in decadal variance.

Code availability
The code used in this study to produce the data analysed was developed in R programming, and can be provided upon reasonable request to BD.