Simulated ENSO asymmetry and nonlinear dynamical heating
The positive skewness of the eastern-Pacific temperature anomalies that characterizes ENSO asymmetry is poorly reproduced among 25 CMIP5 and 26 CMIP6 historical simulations (Supplementary Table 1), despite the fact that simulated ENSO SST amplitude is in a reasonable range compared to observations (Fig. 1e, f). Figure 2a shows that the standard deviation of detrended Niño-3 SST anomalies (σENSO; 150°–90°W, 5°S–5°N) is very close to the observations on average, however, the multi-model mean of the skewness (γENSO) cannot be statistically distinguished from zero. Thus, CMIP climate models fail badly in simulating ENSO skewness that is about 1 in observations, indicating no improvement in CMIP6 compared to earlier CMIP phases16,17,18,31.
Here we show that ENSO SST skewness is highly constrained by subsurface NDH variability in climate models, using multiple reanalysis datasets and the subset of 25 CMIP5 and 18 CMIP6 historical simulations available for evaluating NDH (not all models provide all the necessary ocean data fields that are required to calculate NDH—see “Methods” section and Supplementary Table 1). The subsurface NDH is derived from the nonlinear temperature advective terms in the equatorial eastern-Pacific box (100°W–180°, 1°S–1°N, 50–150 m)37, where an intense NDH gives rise to rectified warming along the mean thermocline and EUC in reanalysis (Fig. 1c, d). In the CMIP climate models, both the mean and variability of the subsurface NDH are too weak (Fig. 1g, h), indicating serious deficiencies in simulating subsurface ocean nonlinear dynamics.
The level of ENSO asymmetry (γENSO) linearly increases with respect to the relative strength of subsurface NDH variability to ENSO amplitude (σNDHsub/σENSO; referred to as NDH efficiency) among the CMIP models with a correlation coefficient of 0.78 (p < 0.00001; Fig. 2b). This explains ~60% of the model-to-model γENSO variance. However, only a few models have an NDH efficiency that is comparable to the observations at 0.30 month−1 on average (Fig. 2b). In contrast, the inter-model spread in surface NDH above 50 m depth does not explain the spread in γENSO (Supplementary Fig. 1). This is supported by an analytical nonlinear ENSO model42, which indicates that surface NDH does not guarantee to enhance the positive ENSO SST skewness. Therefore, we hypothesize that the inability to simulate subsurface NDH is a key dynamical factor responsible for the lack of ENSO asymmetry in climate models.
In addition to the abovementioned deficiencies in the subsurface NDH efficiency, weak ENSO asymmetry in climate models can be attributable to atmospheric nonlinearities in the dynamic and thermodynamic feedback processes of ENSO13,14,25,33,35,43,44. A previous study14 indicates that the El Niño-La Niña asymmetry in shortwave (SW) surface heat flux feedback, characterized by an east-west contrast pattern near the dateline in association with the location of the main convection region, is underestimated in most climate models. Thus, we examine the zonal contrast of equatorial Pacific SW anomalies (140°–170°E minus 140°–170°W, 5°S–5°N) in terms of its regression coefficient difference for positive and negative Niño-3 SST anomalies (∆xSW feedback asymmetry) in Fig. 2c. Indeed, the simulated ENSO SST skewness tends to increase with respect to ∆xSW feedback asymmetry with a moderate correlation coefficient of 0.37 (p = 0.008). However, the simulated ENSO skewness is generally too low even in the climate models that have higher feedback asymmetry. We also confirm that the dynamic feedback asymmetry derived from the zonal wind stress anomalies in the central-Pacific region (CP; 150°E–120°W, 5°S–5°N) is nearly zero in reanalysis datasets despite that most climate models simulate positive asymmetry in the dynamic feedback (i.e., stronger wind feedback for El Niño than La Niña; Fig. 2d). Thus, it does not appear that enhancing either the SW- or wind feedback nonlinearity would help to considerably improve simulated ENSO SST skewness, calling a need for other sources of the nonlinearity. Therefore, we hereafter focus on the subsurface NDH as a dynamic nonlinear source for simulating ENSO asymmetry.
Model biases in the dynamics of ENSO asymmetry
To detect the dynamics responsible for models’ common biases that prevent simulating realistic ENSO asymmetry, we classify CMIP into subgroups (Table 1) using the inter-model fidelity in the subsurface NDH efficiency (Fig. 2b). The group H is composed of 10 CMIP5 and 4 CMIP6 models that have higher levels of NDH efficiency than the multi-model mean (0.16 month−1) and are thus closer to the observations (0.30 month−1). Seven group H models with an NDH efficiency even greater than 0.20 month−1 are further classified as a subset, group HH. The other 29 models are categorized as group L (low NDH efficiency). Figure 3 shows the simulated nonlinearity in each model group in terms of the SST skewness and mean equatorial NDH (see also Supplementary Fig. 2). In group H, the skewness is broadly positive in the eastern- and negative in the western Pacific similar to the reanalysis (Fig. 1a, b), except that simulated positive skewness is still too weak. The group HH models perform even better (Fig. 3c). Furthermore, as seen in the reanalysis (Fig. 1c, d), the positively skewed intense NDH variability near the thermocline results in positive long-term mean residuals (Fig. 3e, f)—a rectification onto the climate mean state9,37. In contrast, the group L models fail to reproduce these key nonlinear properties (Fig. 3a, d).
The group L models, consisting of about 70% of the available CMIP ensemble but severely suffering from poor ENSO nonlinearity, show discernible differences from the reanalysis datasets (Fig. 4). In the composited mean states of group L, the eastern-Pacific cold tongue extends far to the west13,14,15,18,34,45,46 (Supplementary Fig. 3), accompanied by a too intense westward ocean surface current. Nevertheless, the simulated amplitudes of mean EUC and easterly trade winds are comparable with the reanalysis39 (Supplementary Fig. 4). In contrast, the anomalous CP zonal wind stress response to ENSO SST anomalies is too weak (Fig. 5a) even though the SST anomaly amplitude and precipitation response are close to observations (Supplementary Figs. 2 and 4). In addition to this SST-wind coupling bias, most group L models fail to reproduce the anomalous westward EUC response to westerly wind anomalies over the central Pacific39 (Fig. 5b). Indeed, the zonal current covarying with ENSO is less intense in the subsurface but too strong in the western-Pacific surface layer (Supplementary Fig. 4). These two biased linear dynamical coupling processes—from SST to winds and from winds to the EUC—can prevent the generation of wind stress induced subsurface NDH that would normally enhance positive skewness of subsurface temperature and SST anomalies in the eastern equatorial Pacific37 (Fig. 4c, d; see also Supplementary Fig. 5).
In the group H models that better simulate the subsurface NDH efficiency and ENSO SST skewness, these two linear coupling processes (both from SST to winds and from winds to EUC) are improved (Fig. 5) and also the excessive mean cold-tongue bias tends to be reduced (Supplementary Fig. 3). The linear wind-EUC coupling is further improved in group HH, which is a subset of the group H models that have higher NDH efficiency. A common bias still exists in groups H and HH with negative mean NDH in the western-Pacific surface layer (Fig. 3d–f), potentially due to too strong anomalous eastward ocean currents associated with El Niño and a too strong mean-state westward surface current (Fig. 4 and Supplementary Fig. 4). Nevertheless, the group H models perform noticeably better in simulating oceanic nonlinear dynamics as well as linear ENSO feedback processes.
The majority of the climate models fail to simulate the intensity of the linear coupled dynamic feedbacks from SST to winds and from winds to the EUC (Fig. 5). As previous studies suggested12,13,14,15, the errors in the SST-winds coupling tend to be compensated by errors in thermodynamic radiative feedbacks, resulting in a seemingly realistic growth rate and thus ENSO amplitude (Supplementary Fig. 6). However, this error compensation falls apart for ENSO asymmetry because the asymmetry is largely affected by errors in the dynamical coupling alone through the subsurface NDH. Thus, simulating both linear dynamical coupling and thermodynamic feedback correctly (i.e., reducing error compensations that have occurred in climate models) shall improve ENSO asymmetry. It will be important for modeling ENSO nonlinearity to elucidate why the wind stress response to ENSO SST anomalies is too weak and what prevents a realistic EUC response to a given wind forcing.
Constrained tropical response to ENSO amplitude change
The ENSO nonlinearity is a potential source that modulates the current climate mean state4,5,6,9,16,33 and future projections7,23,25,26. However, we have shown in the previous section that the majority of the CMIP climate models fail to reproduce the dynamics of ENSO asymmetry. Thus, we next address whether the future tropical climate response is related to ENSO amplitude change if we constrain the CMIP climate model ensemble to the L, H, and HH subsets. Hence, we use a dynamics-oriented criterium for ENSO asymmetry rather than one based only on SST statistics as is done in previous studies25,26 to hopefully provide a clearer illustration of the consequences of ENSO changes through nonlinear dynamics.
As the group H models are able to simulate the essential properties of subsurface NDH, we expect that in these models an increase of future ENSO amplitude should enhance the mean-state rectification effect that warms the equatorial eastern-Pacific subsurface and deepens the thermocline. In contrast, a future reduction of ENSO amplitude in these models should lead to less rectified warming in the equatorial eastern-Pacific subsurface. Importantly, we do not expect this effect for models that do not realistically simulate NDH (group L). Comparing the historical period to the future scenario simulations for 2051–2100 (RCP8.5 for CMIP5 and SSP5-8.5 for CMIP6; see “Methods” section), we find that the future change in the long-term mean of the subsurface NDH linearly follows ENSO amplitude change in group H (Fig. 6a). However, the projected ENSO amplitude change shows still a large spread even within group H: 7 of these models project a strengthening of ENSO (group H-sEN), while the other 6 models a weakening (group H-wEN). Quantifying ENSO changes in a future warm climate remains a challenge26,27,28 as it requires climate models to simulate ENSO processes realistically in the current climate20. As for the group H-wEN models at least (Supplementary Fig. 7), we confirm that the weakening of ENSO is attributable to the declining oceanic wave response to equatorial-Pacific zonal winds due to an intensified thermal stratification24. Nevertheless, the spread of ENSO amplitude change in group H does give us an opportunity to address the question of how ENSO asymmetry and amplitude together can affect subsurface NDH and thus ENSO’s nonlinear rectification onto the climate mean state.
What kind of surface warming pattern will emerge in the equatorial Pacific in response to greenhouse gas forcing is still a subject of active debate20,23,24,25,27,39,47,48,49,50,51,52. Here we show that the projected warming pattern is strongly related to ENSO amplitude change once climate models are conditioned by their fidelity in dynamics responsible for ENSO asymmetry (as reflected in the simulated NDH efficiency). A measure of El Niño-likeness of tropical Pacific warming, defined as the spatial correlation coefficient of the SST trend pattern in the scenario simulations with the ENSO-regressed SST anomaly pattern in the historical simulations over 90°E–60°W and 20°S–20˚N (see “Methods” section), is positive (i.e., more warming in the east than the west) in most CMIP models20,25,26,48 but highly uncertain ranging from −0.11 to 0.83 (Fig. 6b). In group H, we find a statistically significant increase of the El Niño-likeness as ENSO amplitude increases (r = 0.57, p = 0.042). This relationship is strengthened to a correlation coefficient of 0.89 (p = 0.0073) within the group HH models that have greater NDH efficiency. The El Niño-likeness in group L is positive, ranging from 0.33 to 0.73, and there exists no correlation with ENSO amplitude change (r = 0.00) as expected because a simulated ENSO without realistic NDH is not able to yield a noticeable rectification effect.
How the projected ENSO change affects future trends is detectable in the difference between the strengthening and weakening ENSO (sEN and wEN) simulations once the CMIP ensemble is conditioned by the fidelity to simulated ENSO nonlinear dynamics (Fig. 7). In a CMIP5 greenhouse warming scenario (RCP8.5), the tropical-Pacific surface warming trends of groups H-sEN and H-wEN are significantly different while those of groups L-sEN and L-wEN cannot be distinguished statistically (Fig. 7c, i). For instance, the surface warming in the eastern equatorial Pacific is enhanced in H-sEN whereas it is reduced in H-wEN (Fig. 7a, b), yielding more eastern-Pacific warming in their difference as the response to increasing ENSO amplitude (Fig. 7c). A difference between these two groups is also discernible in the Atlantic and Indian Oceans, as well as many land regions in East Asia, Africa, Australia, South America, and Europe. Correspondingly, precipitation increases more in the eastern equatorial Pacific for H-sEN but in the western Pacific for H-wEN (Supplementary Fig. 8a–c). In contrast, the group L models project El Niño-like tropical warming only regardless of ENSO amplitude change (Figs. 6b and 7g–i). The mean NDH change following the ENSO amplitude change in group H (Fig. 6a) leads to an east-west contrast of subsurface warming difference between H-sEN and H-wEN (Fig. 7d–f) while the mean NDH change is small in group L so that ENSO-related trends are barely discernible (Fig. 7j–l).
The different responses to ENSO amplitude change between groups H and L reveal the importance of constraining the CMIP ensemble by their fidelity of realistically simulating both linear and nonlinear ENSO dynamics for detecting the ENSO mean-state rectification effect. Only if nonlinear ENSO dynamics are captured realistically in climate models, enhancing of the subsurface rectification effect (that warms the eastern-Pacific subsurface temperature, deepens the mean-state thermocline, and thereby increases eastern-Pacific SST) can lead to a more pronounced El Niño-like global-warming pattern, whereas suppressing the rectification effect leads to less eastern Pacific warming. Without simulating ENSO asymmetry and relevant dynamics realistically, however, climate models project an El Niño-like warming pattern only50. This can be attributed to ENSO-independent processes such as reduced equatorial-Pacific mean upwelling due to a weakening of the Walker and Hadley circulations48,52,53, of which the latter is more prominent in group L (Supplementary Fig. 9). This mean state change might in turn alter ENSO properties19,20,21,22,26,54,55,56, but ENSO amplitude change shows no correlation to the warming pattern in group L (Fig. 6b).