Robust Decadal Variations in ENSO Diversity, and its Impact on Future Scenarios

El Niño-Southern Oscillation (ENSO) shows a large diversity of events, whose modulation by 27 climate variability and change, and their representation in climate models, limit our ability to 28 predict their impact on ecosystems and human livelihood. Here, we introduce a new framework 29 to analysze probabilistic changes in event-location and -intensity, which overcomes existing 30 limitations in studying ENSO diversity. We find robust decadal variations in event intensities 31 and locations in century-long observational datasets, which are associated with perturbations 32 in equatorial wind-stress and thermocline depth, as well as extra-tropical anomalies in the 33 North and South Pacific. A large fraction of CMIP5 and CMIP6 models appear capable of 34 simulating such decadal variability in ENSO diversity, and the associated large-scale patterns. 35 Projections of ENSO diversity in future climate change scenarios strongly depend on the 36 magnitude of decadal variations, and the ability of climate models to reproduce them 37 realistically over the 21 st century.


40
El Niño-Southern Oscillation (ENSO) is the leading mode of tropical climate variability, with 41 impacts on ecosystems, agriculture, freshwater supplies and hydropower production spanning 42 much of the globe 1-3 . The majority of impact studies, including seasonal to multi-year 43 predictions, have developed from a "canonical" representation of ENSO, as characterised by 44 sea-surface temperature anomalies (SSTa) in the central-eastern Pacific 4-6 . However, ENSO 45 shows large differences from one event to another in terms of its intensity, spatial patterns and 46 temporal evolutions 7-9 . For instance, while the 1997/98 El Niño displayed extreme SSTa in the 47 eastern equatorial Pacific (EP-ENSO), the largest SSTa in the winter of 2002/03 were weaker 48 and primarily confined to the central equatorial Pacific (CP-ENSO). Differences in the 49 longitudinal location and intensity of ENSO events are sensitively associated with different 50 impacts on regional climate throughout the world 10,11 . Such differences in ENSO patterns, 51 referred to as "ENSO diversity" 7 , and their representation in climate models thus strongly 52 influence the skill of impact prediction systems 12 , and underscore the need for an appropriate 53 characterization and further mechanistic understanding of ENSO diversity. 54

55
The post-1990s increases in the frequency of CP El-Niño events 13-15 has led some researchers 56 to associate such changes in ENSO patterns to the influence of global warming 16  considered more conducive to either CP-or EP-ENSO. However, to date, no robust decadal 62 variation in event locations has been detected in long-term observational records 24,25 . Similarly, 63 we know very little about the ability of climate models to realistically reproduce ENSO 64 diversity and its low-frequency variability 7,8 . CMIP3 and CMIP5 models showed large biases 65 in the background mean state of equatorial Pacific SST, leading to an excessive westward 66 extension of the ENSO patterns 26-28 , and limiting the models' ability to simulate realistic ranges 67 of ENSO diversity 9,29 . Further research is thus needed to determine how these decadal 68 variations may affect the reported historical and future strengthening of ENSO 30-32 , as well as 69 the increasing occurrence of CP-versus EP-events 13-16 . 70

71
Studying ENSO diversity has however been largely limited by various technical shortcomings. 72 Such studies are traditionally based on several indices of tropical Pacific SSTa, designed to 73 monitor the variability of either a "canonical" ENSO or two extreme ENSO flavours 33-37 . Yet, 74 these indices emphasize El Niño events, while La Niña events, which are associated with 75 weaker anomalies than El Niño events 38 , are typically located further West 7 and tend to show 76 a more limited range of pattern diversity 39 , have received less attention. Hence, traditional 77 indices tend to neglect the existing asymmetry between warm-and cold-phases, and between 78 CP-and EP-events. Besides, most of these indices are significantly inter-correlated (Fig.1a), 79 suggesting that they provide redundant information on ENSO. More importantly, the 80 probability distribution of the peak location of SSTa over the tropical Pacific (when each index 81 is exceeding ±0.5°C standard deviation; Fig.1b), indicates that most of these indices are neither 82 representing solely CP-and EP-events, but rather different combinations of both ENSO 83 flavours. 84 #Figure 1# 85 It is therefore essential to develop a new framework allowing for a more precise assessment of 86 changes in the location and intensity of warm-and cold ENSO events. In this regard, the 87 "Center of Heat Index" (CHI) provided the basis for a more flexible framework, allowing the 88 longitudinal centre of SSTa to vary as a function of time instead of being geographically 89 fixed 24,25 . Yet, the CHI approach defines El Niño (La Niña) events based on anomalously warm 90 (cold) areas over very large longitudinal extents (> ~5550km), hereby yielding a very smooth 91 representation of ENSO diversity, and limiting the detectable ranges of spatial patterns. Here, 92 we introduce a novel "non-parametric" framework to allow for a more precise assessment of 93 probabilistic changes in ENSO diversity, notably in event locations and intensity (cf. Methods). 94 This new framework is first applied to a suite of observational datasets, including century-long 95 reconstructions, reanalysis and high-resolution satellite-derived estimates (cf.

113
Observed changes in the likelihood of ENSO location and intensity 114 The average distribution of ENSO location and intensity in 20-year overlapping periods 115 between 1850 and 2017, is examined using five observational datasets (Fig. 2). Hence, observed ENSO diversity is much broader than previously suggested, and exceeds a 149 bimodal view, consistently with results from neural-network based clustering of equatorial 150 Pacific SSTa 44 . This is particularly true for La Niñas, whose diversity was strongly questioned 151 in previous studies 39 . More importantly, ENSO-event locations do not follow a normal 152 distribution describing a randomly distributed continuum of events 24 , converging toward a 153 "canonical" location. Instead, ENSO diversity is clearly linked to low-frequency variations, 154 with multiple preferential locations, which may modulate potential trends in ENSO behaviour.

156
Evaluation of ENSO diversity in climate models within our new framework 157 We compare the simulated probability distributions of ENSO location and intensity from 26 158 CMIP5 and 28 CMIP6 models to our five observational datasets in Fig. 3. (leftmost columns of the CNRM group of models), but only for El Niño events (Fig. 3, top). 167 Other models show highly asymmetrical probability distributions, with clear tendencies to 168 favour either EP/Canonical-or CP-events ( Fig. 3a-b, top), including some extreme cases (e.g. 169 CSIRO-Mk3-6-0) with events always centred further west than the observational range. The In summary, most models simulate some range of pattern diversity for El Niño and La Niña 186 events. In particular, few models present a range of event locations in relatively good 187 agreement with observations, and minimal biases in their intensity. The most common biases 188 concern either the tendency of models to favour one type of events or the events' intensity. 189 Notably, larger biases are found for the models that produce erroneous westward extensions of In observations, SST regressions on El Niño longitude yield an EP-type event, with largest 249 anomalies extending westward along the Equator from the coast of South America (Fig. 6a,   250 top-left). The associated strong westerly wind anomalies extend to the eastern Pacific, where 251 the thermocline is significantly deeper, whilst slight, but significant, easterly wind anomalies 252 and shallower thermocline are found in the western Pacific (Fig. 6a, top-and bottom-left). 253 These patterns indicate that El Niño events tend to be located further east when trade-winds  (Fig. 6c). The North Pacific anomalies associated with ENSO locations are 285 also significant in IPSL-CM6A-LR (Fig. 6c), like in many other models ( Supplementary Fig.   286 3). These relationships between ENSO and Pacific extra-tropical variability however show 287 large ensemble spread in IPSL-CM6A-LR (Fig. 6c), highlighting that these relationships are and La Niña (Fig. 6a-b, top-right). In addition, we found that observed El Niño (La Niña) is 293 more intense when the mean thermocline is deeper (shallower) and the trade-winds are 294 consistently weaker (stronger) over the equatorial Pacific ( Fig. 6a-b, top-and bottom-right). 295 Compared to the large-scale patterns associated with ENSO locations, changes in ENSO 296 intensity are associated with larger wind-stress and thermocline depth anomalies over the 297 central-eastern equatorial Pacific (Fig. 6a-b). ENSO intensity also appears associated with 298 extra-tropical SST and wind anomalies that are more symmetric about the Equator compared 299 to those associated with the location (Fig. 6a- thermocline responses strongly differ from one simulation to another (Fig 6c). According to IPSL-CM6-LR and UKESM-0-LL, most ensemble members converge to more 320 CP-ENSO over the second half of the 21 st century (Fig. 7). This shift to more westward events  (Fig. 7). Such decadal variations remain stronger than potential trends throughout 324 the 21 st century in CNRM-CM5 (Fig. 7). By contrast, the last generation of the same model 325 (i.e. CNRM-CM6-1), which underestimates decadal variability (Fig. 5), shows a clear shift 326 toward more CP-ENSO in the second half of the 21 st century (Supplementary Fig. 6). diversity, according to our framework (Fig. 7). In most models, event-intensity and the 339 frequency of extreme events appear, at least, as variable in the 21 st century as during the 340 historical period (Fig. 7). However, some models, such as IPSL-CM6-LR (in the second half 341 of the 21 st century; Fig. 7, left) and MIROC6 (from the early-to mid-20th century; 342 Supplementary Fig. 6, middle), do show an intensification of ENSO events. In addition, as 343 highlighted in previous studies 30,31 , those models show an increase in the frequency of extreme 344 events ( Fig. 7; Supplementary Fig 6). Although the reliability of MIROC6 simulations is 345 questionable, considering their generally weaker performances in simulating event-location, 346 their results suggest a potential role of anthropogenic climate change in altering ENSO intensity 347 over the 21 st century. Thus, our results highlight that future changes in ENSO characteristics 348 are not necessarily monotonic, as usually assumed, but may undergo large-amplitude decadal 349 variations, leading to the suppression or enhancement of the impact of anthropogenic climate 350 change on ENSO diversity from one decade to another.

392
Observational reference datasets 393 We use five observational datasets, covering all the different ways to reconstruct long-term 394 variability for SST, as well as different resolutions (Supplementary Table 1      Variability in ENSO's, observed and simulated, most likely location and intensity. (a) Maximum power spectrums of the running 10-year El Niño most likely location and intensity (i.e. the mode), as determined using continuous wavelet analysis, and using four long-term observational reference datasets (left: ERSST.v5, COBESST.v2, HadSST1, SODA.si3), 95/250 CMIP5/6 historical runs (middle), as well as 26/28 CMIP5/6 piControl runs (right). (b) same as (a) but for La Niña. Signi cance of variability patches are tested at p=0.05 based on 1000 Monte-Carlo simulations of the red noise background spectrum. 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 spectrums are estimated using the full length of each time series.
The maximum power spectrums are weighted by the signi cance, and only signi cant variability is shown. Figure 5 CMIP5/6 bias in decadal variability of ENSO's most likely location and intensity. (a) Average ratio of standard deviation (rSD) between historical runs and observed decadal variance (>10 year) in the running 10-year most likely location and intensity of El Niño and La Niña events. (b) same as (a) but using pi-Control runs. Statistical signi cance is assessed by performing a two-sided Fisher's F-test at p=0.05 between every 100-yr segments through the course of climate simulations and every 100-yr segments in the four longer-term observational SST datasets (i.e. 27,740 ≤ n ≤ 209,000 replicates), to quantify a rate of success (i.e. the number of times observations and simulations showed equal variance). Black dots highlight simulations for which the rate of success is lower than 10%, showing signi cantly different variance at p=0.1.

Figure 6
Large-scale patterns driving long-term variability in ENSO location and intensity. (a) Observed regressed SST (blue to red shades), wind-stress (vectors) and Z20 anomalies (lines) associated with changes in El Niño location (right) and intensity (left) and using multiple observational data sets (SST: ERSST.v5, COBESST.v2, HadSST1 and SODA.si3; wind-stress: NOAA-20CR.v3; Z20: SODA.v2.2.4). (b) Same as (a) but for La Niña events. (c) same as (a-b) but using the IPSL-CM6A-LR large ensemble (32 members). While SST and wind-stress anomalies are displayed at the pan-Paci c scale based on the median changes in observations, simulated regressed anomalies are assessed through the ensemble median (top) and ensemble spread (standard deviation [SD]; middle]). Z20 anomalies are estimated through the median changes between 5°S and 5°N (bottom). Red and Blue shades on the Z20 anomalies indicate the spread between the four SST observational data sets and within the IPSL-CM6A-LR large ensemble (light to dark: maximum/minimum, 10/90th, 30/70th and 45/55th 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 in uence the ensemble spread. Statistical signi cance is assessed at p=0.05 using 1000 permutations, and displayed as black contour for SSTa, and blue/red crosses for Z20 anomalies. Only signi cant wind-stress anomalies at p=0.05 are displayed.  trend, to allow comparison with observations. The same baseline period was used to estimate the 90th and 10th percentiles.