ENSO atmospheric feedbacks under global warming and their relation to mean-state changes

Two atmospheric feedbacks play an important role in the dynamics of the El Niño/Southern Oscillation (ENSO), namely the amplifying zonal wind feedback and the damping heat flux feedback. Here we investigate how and why both feedbacks change under global warming in climate models participating in the 5th and 6th phase of the Coupled Model Intercomparison Project (CMIP5 and CMIP6) under the business-as-usual scenario (RCP8.5 and SSP5-8.5, respectively). The amplifying zonal wind feedback over the western equatorial Pacific (WEP) becomes significantly stronger in two third of the models, on average by 12 ± 7% in these models. The heat flux damping feedback over the eastern and central equatorial Pacific (EEP and CEP, respectively) increases as well in nearly all models, with the damping effect increasing on average by 18 ± 11%. The simultaneous strengthening of the two feedbacks can be explained by the stronger warming in the EEP relative to the WEP and the off-equatorial regions, which shifts the rising branch of the Pacific Walker Circulation to the east and increases the mean convection over the CEP. This in turn leads to a stronger vertical wind response during ENSO events over the CEP that strengthens both atmospheric feedbacks. We separate the climate models into sub-ensembles with STRONG and WEAK ENSO atmospheric feedbacks, as 2/3 of the models underestimate both feedbacks under present-day conditions by more than 40%, causing an error compensation in the ENSO dynamics. The biased mean state in WEAK in 20C constrains the ENSO atmospheric feedback projection in 21C, as the models of the WEAK sub-ensemble also have weaker ENSO atmospheric feedbacks in 21C. Further, due to the more realistic dynamics and teleconnections, we postulate that one should have more confidence in the ENSO predictions with models belonging to the STRONG sub-ensemble. Finally, we analyze the relation between ENSO amplitude change and ENSO atmospheric feedback change. We find that models simulating an eastward shift of the zonal wind feedback and increasing precipitation over the EEP during Eastern Pacific El Niño events tend to predict a larger ENSO amplitude in response to global warming.


Introduction
The El Niño/Southern Oscillation (ENSO), originating in the tropical Pacific, is the most dominant climate mode on interannual timescales, with major ecological and socio-economic impacts in the Pacific region and beyond (Philander 1990;McPhaden 1999). ENSO events are characterized by large variations in the sea surface temperature (SST) of the eastern and central tropical Pacific. The ENSO-warm phase is termed El Niño, its cold phase La Niña. Large progresses have been made in ENSO research during the last years. It is still unclear, however, how ENSO properties such as its amplitude, will respond to global warming (Van Oldenborgh et al. 2005;Meehl et al. 2007;Latif and Keenlyside 2009;DiNezio et al. 2012;Stocker et al. 2013;Kim et al. 2014a;Cai et al. 2015Cai et al. , 2020bChen et al. 2017;Beobide Arsuaga et al. 2021). The growth and decay of ENSO events are governed by an interplay of positive (amplifying) and negative (damping) feedbacks (Jin et al. 2006). The two most important atmospheric feedbacks are the positive wind-SST feedback and the negative heat flux-SST feedback (Lloyd et al. 2009). During El Niño anomalously warm SST in the eastern and central equatorial Pacific (EEP and CEP, respectively) shift the convection over the western equatorial Pacific (WEP) to the east, which in turn leads to positive zonal near surface wind (U) anomalies over this region, constituting the positive U feedback (Bjerknes 1969;Jansen et al. 2009). This in turn leads to a deepening of the thermocline in the EEP that reinforces the initial SST anomaly. Further, the eastward shift of the convection decreases the shortwave radiation (SW), resulting in a negative SW feedback that is strongest over the WEP (Lloyd et al. 2009). Due to the warmer SST the latent heat (LH) loss over the EEP increases, resulting in a negative LH feedback (Lloyd et al. 2011). Together, the SW and LH feedbacks form the negative net heat flux (Q net ) feedback over the entire equatorial Pacific (Lloyd et al. 2009).
El Niño and its cold phase, La Niña, in general involve the same positive and negative feedbacks, but show important asymmetries (Guan et al. 2019). El Niño events are on average stronger and the SST anomalies are further in the east than during La Niña events . Further, the center of El Niño events varies considerably between the EEP and CEP. The terms Eastern Pacific (EP) El Niño and Central Pacific (CP) El Niño were introduced to classify El Niño events and the EP El Niño events are on average stronger than the CP El Niño events (Kao and Yu 2009;Takahashi et al. 2011;Capotondi et al. 2014;Timmermann et al. 2018). Another important feature of ENSO is its phase locking to the annual cycle. ENSO events typically peak around the end of the calendar year (Tziperman et al. 1998).
Many state-of-the-art climate models from the 5th and 6th phase of the Coupled Model Intercomparison Project (CMIP5 and CMIP6) strongly underestimate both ENSO atmospheric feedbacks (U and Q net feedback, hereafter termed ENAF) Bellenger et al. 2014;Planton et al. 2020). This can be explained to a large extent by the equatorial cold SST bias, which results in a La Niñalike mean state with the rising branch of the Pacific Walker Circulation (PWC) too far in the west (Lloyd et al. 2009(Lloyd et al. , 2011(Lloyd et al. , 2012Bayr et al. 2018Bayr et al. , 2020Ferrett et al. 2018). As both ENAF strongly depend on the position of the rising branch of the PWC, an error compensation between the underestimated U and Q net feedbacks is observed in many climate models Kim et al. 2014b;Bayr et al. 2019bBayr et al. , 2020. This further results in biased ENSO dynamics; roughly half of the CMIP models simulate a hybrid of wind-driven and shortwave-driven ENSO dynamics and not purely wind-driven ENSO dynamics as observed (Dommenget 2010;Dommenget et al. 2014;Bayr et al. 2019bBayr et al. , 2021. Finally, the biased ENSO dynamics can also explain why many climate models have problems in simulating important ENSO properties like the asymmetry and phase locking of ENSO (Wengel et al. 2018;Bayr et al. 2019bBayr et al. , 2021, or fail to realistically simulate the ENSO teleconnections to e.g. the North Pacific (Bayr et al. 2019a).
There is an ongoing debate on how ENSO will change under global warming, with, for example, little agreement in ENSO amplitude change in the CMIP models (Cai et al. 2015(Cai et al. , 2020b. The lack of consistency among models can be partly explained by competing effects of amplifying and damping feedbacks in the ENSO dynamics (Dommenget and Vijayeta 2019;Cai et al. 2020bCai et al. , 2021. Early studies suggested an increase in the ENSO amplitude due to a larger sensitivity of the SST to ENSO-related wind anomalies due to increased ocean stratification (Timmermann et al. 1999;Collins 2000). In later versions of the climate models, it became evident that also the damping feedbacks strengthen under global warming. Therefore the overall amplitude change depends on the balance of the changes in the amplifying and damping feedbacks (Cai et al. 2020b). Little attention has been paid so far on how and why both ENAF will change under global warming. Therefore, our focus here is to investigate the U and Q net feedback changes under global warming and the implications for reducing ENSO-response uncertainty. Further, in recent studies mounting evidence was found that the climate models with the most realistic simulation of ENSO dynamics under present-day conditions may more reliably project ENSO changes under global warming (Cai et al. 2020a(Cai et al. , 2021Hayashi et al. 2020). To this end, we separate the climate models into two sub-ensembles, one comprising models with relatively realistic and the other models with less realistic ENAFs, and describe the response differences between the two sub-ensembles. The manuscript is organized as follows: in Sect. 2, we describe the data and methods used in this study and in Sect. 3, the change in ENAF under global warming. The change in the atmospheric mean state is analyzed in Sect. 4 and the relation between atmospheric mean state and ENAF changes in Sect. 5. This study is summarized in Sect. 6.

Data and methods
Observed SSTs are taken from HadISST for the period 1979-2018 (Rayner et al. 2003). Near-surface (10 m) zonal wind (U), tendency of air pressure at 500 hPa, here termed vertical wind (Ω), and heat fluxes for the period 1979-2018 are taken from ERA-Interim reanalysis (Simmons et al. 2007). Observed precipitation for the period 1979-2018 is taken from CMAP data set (Xie and Arkin 1997).
The data is interpolated onto a regular 2.5° × 2.5° grid and we use all models with the required data available (see the legend of Fig. 1 for the models).
The Niño1 + 2 region is defined as 80°W-90°W and 10°S-0°, the Niño3 region as 90°W-150°W and 5°S-5°N, the Niño3.4 region as 120°W-170°W and 5°S-5°N, the (e) (f) (g) Fig. 1 Atmospheric feedbacks in ERA Interim reanalysis, CMIP5 and CMIP6 models, in a 10 m zonal wind (U) feedback in the Nino4 region on the x-axis vs. net heat flux (Q net ) feedback in the combined Nino3 and Nino4 region on the y-axis, both for the historical scenario (1950-1999 = > 20C); in b same as a but here for the RCP8.5 (CMIP5) and SSP5-8.5 (CMIP6) scenarios (2050-2099 = > 21C); in c) U feedback, in 20C on the x-axis vs. 21C on the y-axis; in d same as c but here for the Q net feedback; in e the atmospheric feedback strength, i.e. the average of the U feedback and Q net feedback, both normalized by ERA Interim value, on the x-axis vs. the relative SST in the Nino3.4 region on the y-axis; in f same as e, but here for the 21C; in g same as e, but here the difference 21C − 20C; the red and blue color of the numbers indicate the models of the STRONG and WEAK sub-ensemble, respectively, and the triangles indicate the sub-ensemble mean; The solid lines in a, b, and e-g are the regression lines for the CMIP models and the correlation between the CMIP models is given in the upper right or left corner, with one, two or three stars indicating a significant correlation on a 90%, 95% or 99% confidence level, respectively; The dashed lines in c and d are the unity line, so that models on that line have similar strength of the corresponding feedback in both time periods Niño4 region as 160°E-150°W and 5°S-5°N and the Niño5 region as 130°E-160°E and 0°-10°N. Monthly-mean values are used in the analyses below. Anomalies are calculated with respect to the climatological seasonal cycle and the quadratic trend is subtracted. We define the U feedback as the SST anomalies (SSTA) averaged over the Niño3.4 region regressed on the U anomalies over the western equatorial Pacific (Niño4 region). We use here U instead of zonal wind stress for most models, as U is more comparable due to differences in wind-stress calculation between the models (see Bayr et al. 2019b for more details). For six models, which are indicated by a star in the legend of Fig. 1), U is not available and zonal wind stress is used instead. Wind stress is shown in units of 10 -2 N/ m 2 , which is roughly comparable to 1 m/s of U (Bayr et al. 2019b). The surface net heat flux (Q net ) feedback is defined as the SSTA of Niño3.4 regressed on Q net over the eastern equatorial (Niño3 region) and western equatorial Pacific (Niño4 region). The asymmetry of ENSO is shown here as the difference in skewness of the SSTA in Niño3 and Niño4 region; and phase locking is measured by the standard deviation of the SSTA in Niño3.4 in November to January (NDJ, the months with the highest variability in observations) divided by the SSTA in April to June (AMJ, the months with the lowest variability in observations).
The contribution of the thermodynamics and ocean dynamics to the SST change during ENSO growth in the Niño3.4 region is estimated by a method proposed by Bayr et al. (2019bBayr et al. ( , 2021: the net heat flux Q net is integrated over the 6 months preceding the maximum of each ENSO event: where c p = 4000 J kg −1 K −1 is the specific heat capacity at constant pressure of seawater, ρ = 1024 kg m −3 the average density of seawater, H = 50 m the mixed layer depth (MLD) in meters, and t the time in months. The MLD is set to a constant value of 50 m, as in the Niño3.4 region the differences in MLD between the models have a negligible influence on the results of this method . We normalize the SST change (dSST) and the SST change caused by the heat fluxes (dSST HF ) by dSST at the peak time (t = 0), yielding an SST change per Kelvin warming and enabling comparison of the individual simulations irrespective of their amplitude. The integration is also performed separately for the shortwave (SW), longwave (LW), sensible (SH) and latent heat (LH) flux components in order to estimate their contribution to the total SST tendency. Finally, as the observed SST change is the balance of the SST changes due to ocean dynamics (dSST oc ) and net heat flux (dSST HF ), we can estimate the contribution of the ocean dynamics from the results above: In this way, it is possible to quantify the contributions of thermodynamics and ocean dynamics during ENSO growth and to quantify if ENSO is purely ocean dynamics-driven, like in observations, or also partly shortwavedriven, as in many climate models with an overly strong equatorial cold tongue. We note that dSST HF is also influenced by the surface winds as the LH and SH components strongly depend on the surface wind speed. But as shown in Bayr et al. (2019b), this method is suited to quantify the relative contributions of thermodynamics and ocean dynamics to ENSO growth when compared with e.g. more sophisticated methods like the BJ index (Jin et al. 2006). To estimate the uncertainty in the SST change by the heat fluxes, we apply a bootstrapping approach. We randomly choose 1,000 times only two thirds of the ENSO events and calculate the SST change. The uncertainty is shown as error bars, indicating the 90% quantile of the estimated values.
Climate models tend to simulate different global mean temperatures, which shifts the threshold of convection . Therefore, we calculate the relative SST (rSST) by subtracting the area averaged SST over the tropical Pacific (120°E-70°W, 15°S-15°N). Further, due to moistadiabatic adjustment of the tropical atmosphere, the relative SST is more appropriate to describe the relation between SST and convection than the absolute SSTs (Johnson and Xie 2010;He et al. 2018;Izumo et al. 2019). We use the relative SST over the Niño3.4 region as a reference, as this region has the strongest influence on the ENAF strength in the coupled models (Bayr et al. 2018).
As a measure of the zonal circulation along the equator, we use the zonal stream function as defined in Yu and Zwiers (2010) and Bayr et al. (2014): Here u D is the divergent component of the zonal wind, , a the radius of the Earth, p the pressure and g the gravity. The zonal wind is averaged over the meridional band 5°N-5°S and integrated from the top of the atmosphere to the surface. In the figures shown below, only the levels below 100 hPa are depicted, as the stream function is nearly zero above that level.

Change in atmospheric feedbacks under global warming
The U feedback and Q net feedback vary considerably among the CMIP models in 20C (Fig. 1a), with some models having feedback strengths close to ERA-Interim and others strongly underestimating both the U feedback and Q net feedback. The strong linear relation between the two feedbacks reveals the error compensation between the too weak wind forcing and the too weak heat flux damping seen in many climate models Kim et al. 2014b;Bayr et al. 2019b). We define the ENAF strength of each model as the average strength of U and Q net feedback, after normalizing each by the ERA-Interim value. The ENAF strength of the CMIP models varies between 21 and 84% of the ERA-Interim strength. We separate the models into two sub-ensembles: models with ENAF strength of greater than 60% of the ERA-Interim are in the STRONG sub-ensemble and the others in the WEAK sub-ensemble, as indicated by the red and blue color in Fig. 1a), respectively. This allows us to highlight the differences in the global warming response between models with unrealistic (WEAK sub-ensemble) and more realistic ENAF and ENSO dynamics (STRONG sub-ensemble). As shown in Bayr et al. (2019b), models with strongly underestimated ENAF have, due to the error compensation, strongly biased ENSO dynamics, as ENSO is partly driven by a positive SW feedback, which compensates the too weak wind forcing. These models also tend to exhibit biases in other important ENSO properties such as the asymmetry between El Niño and La Niña (Fig. 2a) or the phase locking of ENSO to the seasonal cycle (Fig. 2c). Under global warming (21C), the linear relation between the U and Q net feedback becomes slightly weaker but remains highly significant (Fig. 1b). The sub-ensemble means (red and blue triangles in Fig. 1a,b) depict an increase in both feedbacks in 21C relative to 20C. More specifically, Fig. 1c) shows that the U feedback strength increases significantly-statistical significance is measured by the 95% confidence level-in 32 out of 47 models (68%), while 3 models simulate a significant decrease and 12 models no significant change. The sub-ensemble mean increases from 1.05 m/s/K in 20C to 1.17 m/s/K in STRONG and from 0.63 to 0.73 m/s/K in WEAK. The Q net feedback increases significantly in 44 out of 47 models, while only 3 models simulate no significant change (Fig. 1d). The sub-ensemble mean increases from − 11.5 W/m 2 /K in 20C to − 13.9 W/m 2 /K in 21C in STRONG and from − 5.1 to − 8.4 W/m 2 /K in WEAK (see also Table 1). The ENSO asymmetry in 21C also is considerably larger in the models with stronger ENAF than in models with weak ENAF (Fig. 2b), while the phase locking shows a weak but still significant correlation (Fig. 2d).
In 20C, the ENAF strength is strongly related to the relative SST (rSST) in the Central Pacific (Niño3.4 region, Fig. 1e). As shown in Bayr et al. (2018Bayr et al. ( , 2020, the meanstate rSST in this region determines the position of the rising branch of the PWC and therefore the atmospheric response to SSTA, which in turn determines the strength of both feedbacks. We hypothesize that the change in rSST in Niño3.4 also determines the change of ENAF in 21C. Indeed, a strong relation between the rSST and the ENAF strength can be observed in 21C (Fig. 1f) and in the difference between 21 and 20C too (Fig. 1g). Further, only two models simulate a slight decrease in the ENAF strength (Fig. 1g), while all other models show an increase, with quite a large spread regarding the level of increase amongst the models. Finally, both sub-ensembles exhibit an overall large spread in ENAF change, i.e. the change in the two feedbacks seems to be independent from the initial ENAF feedback strength in 20C. We will discuss in the next sections, how the changes in the U and Q net feedback under global warming can be explained by the change in rSST and the associated atmospheric mean-state changes.

Mean state changes
In the area average, the tropical Pacific (120°E-70°W, 15°S-15°N) warms from 1950-1999 to 2050-2099 by 2.5 ± 0.5 K in all climate models. The SST change relative to this area-averaged warming is shown as shading, in Fig. 3a) for all models, and in Fig. 3b, c) for the STRONG and WEAK sub-ensembles, respectively. Relative to the area-averaged warming, the EEP warms strongest, while the far WEP and especially the southern tropical Pacific become colder irrespective of the ensemble (Fig. 3a-c). In the two sub-ensembles the rSST change in EEP and southern tropical Pacific is comparable, while in the WEP the warming is much weaker in STRONG than in WEAK (Fig. 3b, c). The mean state in 20C, represented by the contours in Fig. 3a-c), reveals that the models in the WEAK sub-ensemble exhibit the larger cold tongue bias, as already indicated in Fig. 1e). Thus, models with a larger cold tongue bias tend to simulate a stronger warming in the WEP and therefore a smaller reduction of the SST gradient between the Niño4 and Niño3 regions (Fig. 3d). On the other hand, models with a small cold tongue bias tend to simulate a smaller warming in the WEP and thus a larger reduction of the SST gradient between the Niño4 and Niño3 regions. Figure 4 reveals the equatorial mean states (5°S-5°N) of STRONG and WEAK in 20C and 21C and allows a comparison with observations and reanalysis data. The equatorial mean rSST (Fig. 4a) in 20C is quite close to observed SST in STRONG, especially in Niño3.4 region (see also Table 1), while it is considerably colder in WEAK by − 0.7 K in Niño3.4 region. In 21C, the change in rSST reduces the SST gradient in both sub-ensembles, but more in STRONG than in WEAK. The equatorial mean U (Fig. 4b) over the WEP in 20C is considerably stronger in WEAK than in STRONG, while it is the other way around in the EEP. U differences between the two sub-ensembles are similar in 20C and in 21C and both sub-ensembles show a weakening of U under global warming.
The mean vertical wind at 500 hPa (Ω) (Fig. 4c) in 20C is more realistic in STRONG, as the ascent over the Niño4 region is much too weak in WEAK due to the more pronounced cold tongue bias. In both sub-ensembles, the ascent becomes stronger in 21C over the entire equatorial Pacific, which is due to the stronger warming over the equatorial than over the off-equatorial region. Finally, mean precipitation ( Fig. 4d) in 20C also is more realistic in STRONG than in WEAK, where the STRONG sub-ensemble overestimates the precipitation over the EEP and underestimates the precipitation over the WEP. The WEAK sub-ensemble strongly underestimates the equatorial precipitation over both the CEP and WEP. Under global warming, both subensembles again show a similar increase in mean equatorial In a same as Fig. 1e, but here for the ENSO asymmetry, measured by the difference in skewness of SST between the Nino3 and Nino4 region, on the y-axis; in b same as a but here for 21C; c same as a, but here for the ENSO phase locking to the annual cycle, measured by the standard deviation of SST in NDJ divided by standard deviation of SST in AMJ, on the y-axis; d same as c, but here for 21C 151.5 149.7 (± 10.5) 158.6 (± 13.7) + 9.2 (± 10.5) 139.0 (± 9.6) 148.0 (± 13.0) + 8.3 (± 7.5) Feedbacks U in Niño4 (m/s/K) + 1.43 + 1.05 (± 0.12) + 1.17 (± 0.17) + 0.12 (± 0.13) + 0.63 (± 0.10) + 0.73 (± 0.14) + 0.11 (± 0.12) Ω in CEP (10 - rSST 20C vs. SST grad. change values as indicated by a T-test on a 95% confidence level are indicated by stippling; d same as Fig. 1e), but here the change 21C − 20C in SST gradient between the Nino3 and Nino4 region on the y-axis. The dashed horizontal line marks the rSST in the Nino3.4 region in HadISST precipitation, consistent with the stronger ascent (Fig. 4c). This is in agreement with previous studies (Vecchi and Soden 2007;Cai et al. 2014;Sohn et al. 2019). The mean state and its changes in STRONG and WEAK are also given in Table 1.
Next, we discuss the Pacific Walker Circulation (PWC) in reanalysis data, the two sub-ensembles and the individual models (Fig. 5). The contours in Fig. 5a-c) depict the meanstate PWC in ERA-Interim, the shading shows the 20C mean state in ERA-Interim (Fig. 5a), STRONG (Fig. 5b) and WEAK (Fig. 5c). The PWC is closer to ERA-Interim in STRONG than in WEAK. In particular, the rising branch of the PWC, given by the zero contour over the WEP, is simulated too far west by about 12° in WEAK and quite realistic in STRONG. This is also seen in Fig. 5f), depicting the position of the rising branch of the PWC on the x-axis.
Under global warming, the PWC weakens and shifts to the east in both sub-ensembles (Fig. 5d, e), consistent with the weakened SST gradient (Fig. 4a). The rising branch of the PWC shifts eastward by roughly 10° in both sub-ensembles. The similar change in both sub-ensembles may be explained by the compensating effects between the warming of the rSST in the Niño3.4 region and the reduction of the SST gradient across the equatorial Pacific, as the first is larger in WEAK (Fig. 5h) (Fig. 3). The two sub-ensembles exhibit a similar spread amongst the individual members, ranging from a slight westward shift of 5° to a strong eastward shift of 25° (Fig. 5h). Finally, we find in 21C a similar relation between the rising branch of the PWC and the rSST in the Niño3.4 region ( Fig. 5g), indicating that the rSST in the Niño3.4 region is also quite important for the mean state of the PWC in 21C.
In summary, the stronger SST warming in the equatorial region relative to the off-equatorial regions increases the upward motion and precipitation in the equatorial region  Fig. 1e), but here the longitude of the rising branch of the WC, as measured by the crossing of the zero line of the pressure weighted vertical mean of the zonal stream function, in 20C for ERA Interim and the individual CMIP models on the y-axis; g same as f, but here for the 21C; h same as f, but here for the difference 21C − 20C under global warming. Further, the weakened SST gradient across the equator and the warming of the rSST in the Niño3.4 region goes along with an eastward shift of the PWC and a reduction in the mean U over the Niño4 region.
In the next section, we investigate how the atmospheric mean-state changes under global warming can explain the change in atmospheric feedbacks.

Relation between mean state changes and ENAF changes
As already indicated by Figs. 4 and 5, in 20C the rSST in the Niño3.4 region determines the mean vertical wind strength at 500 hPa in the Niño4 region, which can be seen for the individual models in Fig. 6a). The models of the STRONG sub-ensemble have, due to the warmer rSST in  Fig. 1e), but here for the mean vertical wind at 500 hPa (Ω) in the Nino4 region in 20C on the x-axis; b same as a but here for the 21C; c same as a, but here for the difference 21C − 20C; d-f same as a-c, but here for the Ω feedback in the central equatorial Pacific (CEP; 180°-150° W, 5° S-5° N) on the y-axis; g-i same as d-f, but here for U feedback in the Nino4 region on the x-axis the Niño3.4 region, a stronger ascent in the Niño4 region than the models of the WEAK sub-ensemble. The same relation holds in 21C and for the difference 21C − 20C (Fig. 6b, c). The mean vertical wind is strongly related to the vertical wind feedback in the CEP (180°-150°W, 5°S-5°N) in 20C, in 21C and in the difference 21C − 20C (Fig. 6d-f). This is consistent with Bayr et al. (2020), showing that the position of the rising branch of the PWC strongly influences the position and strength of the vertical wind feedback. Finally, the U feedback in the Niño4 region strongly depends on the vertical wind feedback in the CEP in 20C, in 21C and in the difference 21C − 20C (Fig. 6g-i), as the eastward (westward) shift in convection over the CEP during El Niño (La Niña) causes the weakening (strengthening) of U in the Niño4 region. In 20C as well as in 21C, the models of the STRONG sub-ensemble simulate, due to the stronger mean ascent in the Niño4 region, a stronger vertical wind feedback in the CEP and therefore a stronger U feedback in the Niño4 region than the models of the WEAK sub-ensemble (see also Table 1). Under global warming, both sub-ensembles simulate a quite similar increase in the U feedback of ~ 0.1 m/s/K in the Niño4 region, while the vertical wind feedback over the CEP increases more in WEAK than in STRONG (Fig. 6i). The spread among the individual members is quite similar in the two sub-ensembles.
The Q net feedback in 20C also strongly depends on the atmospheric mean state over the equatorial Pacific (Bayr et al. 2018(Bayr et al. , 2020, mainly on the mean convection and mean cloud cover. Therefore, the Q net feedback over the combined Niño3 and Niño4 region is strongly influenced by the vertical wind feedback strength over the CEP in both 20C and 21C, as well as the difference 21C − 20C (Fig. 7a-c). This can be explained by the strong link between the Q net feedback and the SW feedback over the combined Niño3 and Niño4 region in 20C, in 21C and in the difference 21C − 20C (Fig. 7d-f), which reflects the cloud-cover changes during ENSO events. The LH feedback has a much weaker influence on the Q net feedback in both periods and only is important in the Niño3 region (Fig. 7g-i), where the LH feedback is strongest. As a result of the more realistic mean state, the STRONG sub-ensemble has the more realistic Q net feedback than the WEAK sub-ensemble in 20C, mainly due to the stronger SW feedback (Fig. 7d). Both sub-ensembles show an increase in Q net feedback strength in response to global warming, which in WEAK with − 3.4 W/m 2 /K is larger than in STRONG with − 2.4 W/m 2 /K. The spread within the sub-ensembles is slightly larger in WEAK than in STRONG. Despite the larger increase in WEAK, the absolute Q net feedback in 21C is still much larger in STRONG, with − 13.9 W/m 2 /K, than in WEAK, with − 8.4 W/m 2 /K.
The differences between the two sub-ensembles become clearer when analyzing the spatial patterns of the atmospheric feedbacks. In 20C, all five feedbacks are simulated more realistically in STRONG than in WEAK (Fig. 8) both in terms of amplitude (shading) and significance (regions with explained variance of regression < 0.2 are stippled), even though the U, Q net , SW and LH feedbacks are slightly weaker in STRONG than in ERA-Interim. As described in Bayr et al. (2018Bayr et al. ( , 2020, the stronger cold tongue in WEAK shifts the rising branch of the Walker Circulation to the west, which drives the vertical wind, the U and SW feedbacks too far westward and weakens the atmospheric response (Fig. 8c, f, l). Further, the SW feedback in the Niño3 region is positive in WEAK, while it is clearly negative in ERA-Interim and STRONG (Fig. 8j-l). This can be explained by an overestimation of the low-level stratiform clouds in case of a stronger cold tongue, which dissolve during an El Niño event (Lloyd et al. 2009). This leads to the much weaker Q net feedback in WEAK compared to STRONG (Fig. 8g-i).
In 21C, all five feedbacks become stronger in both subensembles ( Fig. 9), but all feedbacks are still much stronger in STRONG than in WEAK (see also Table 1). Further, the vertical wind, U and SW feedbacks are located farther in the east in STRONG, which is due to the different mean states. The differences 21C − 20C (Fig. 10) reveal some important spatial differences in the global warming response between the two sub-ensembles. Although both sub-ensembles show a similar dipole pattern in the vertical wind feedback change, the increase in ascent is more in the east in STRONG than in WEAK (Fig. 10a, b). This forces the largest zonal U feedback change being slightly more eastward in STRONG than in WEAK (Fig. 10c, d). The Q net , SW and LH feedback changes in the negative pole also are further east in STRONG than in WEAK (Fig. 10e-j).
In summary, due to the different atmospheric mean states the atmospheric feedbacks are much weaker in WEAK than in STRONG in 20C as well as in 21C. Therefore, we next investigate how the ENSO dynamics change in the two subensembles. Previous studies have shown that models with strongly underestimated U and Q net feedbacks under presentday conditions exhibit biased ENSO dynamics that are hard to uncover on first glance due to error compensation Dommenget et al. 2014;Bayr et al. 2019b). The biased ENSO dynamics can be revealed by analyzing the interplay of ocean dynamics and thermodynamics during ENSO growth phase (Bayr et al. 2019b. In reanalysis data, 1 K of SST change during an ENSO event in the Niño3.4 region goes along with a forcing of + 2.3 K/K by the ocean dynamics and a Q net damping of − 1.3 K/K, mainly by the SW feedback amounting to − 0.8 K/K and the LH feedback amounting to − 0.5 K/K (Fig. 11a). The STRONG sub-ensemble behaves quite similar, with a dynamical forcing of + 2.2 K/K and a Q net damping of − 1.2 K/K, that consists mostly of − 0.8 K/K by the SW feedback and − 0.5 K/K by the LH feedback. The WEAK sub-ensemble has quite different ENSO dynamics, as for 1 K of SST change during an ENSO event only a dynamical forcing of + 1.2 K/K is required and a Q net damping of − 0.2 K/K. The SW feedback acts in WEAK as a forcing and amounts to + 0.4 K/K and the LH feedback with − 0.5 K/K acts as a damping. Thus the models of the WEAK sub-ensemble have a hybrid of wind-driven and shortwave-driven ENSO dynamics, and this explains why these models do not realistically simulate important ENSO properties such as the ENSO asymmetry and the ENSO phase looking (Fig. 2) Bayr et al. 2019bBayr et al. , 2021.
We next compare the ENSO dynamics in the two subensembles in 20C and 21C (Fig. 11a, b), and the difference   Fig. 7 a-c same as Fig. 6g-i), but here for Q net feedback in the Nino3 and Nino4 region on the x-axis; d-f same as a-c, but here for the SW feedback in the Nino3 and Nino4 region on the y-axis; g-i same as a-c, but here for the Q net feedback in Nino3 on the x-axis and the LH feedback in Nino3 on the y-axis 21C − 20C (Fig. 11c). In both sub-ensembles, the forcing by ocean dynamics and heat flux damping become stronger under global warming, in WEAK a bit more than in STRONG (Fig. 11c). The stronger heat flux damping in both sub-ensembles is mainly caused by the SW feedback and to a lesser extent by the LH feedback. In 21C, the forcing by ocean dynamics and the damping by Q net is still much stronger in STRONG, with + 2.5 K/K and − 1.5 K/K, respectively, than in WEAK with + 1.7 K/K and − 0.7 K/K, respectively. The SW feedback explains the difference, as it amounts in 21C to − 1.0 K/K in STRONG and 0.0 K/K in WEAK, while the LH feedback is similar in both sub-ensembles amounting to − 0.6 K/K. Thus, the two sub-ensembles also exhibit considerably different ENSO dynamics in 21C. As the models in WEAK have biased ENSO dynamics in 20C, it is questionable if they can properly simulate how ENSO will change in the future (Cai et al. 2021).
There is an ongoing debate about how the ENSO amplitude will change under global warming. Can the change in atmospheric mean state and atmospheric feedbacks explain the large spread in ENSO-amplitude change under global warming (Fig. 12)? The spread cannot be explained by the U feedback change in the Niño4 region, since the correlation between the U feedback change in the Niño4 region and the ENSO-amplitude change is 0.05 (not shown). There is, however, a significant positive correlation of + 0.48 between the ENSO-amplitude change and the east minus west difference of the U feedback change (Niño3 region minus Niño5 region, where the Niño5 region is defined as 130°E-160°E, 0°-10°N, Fig. 12a). This indicates that the ENSO amplitude increases, if the U feedback shifts to the east, i.e. getting weaker over the WEP and stronger over the EEP. This relation becomes stronger in the STRONG subensemble, with + 0.57, and lower in the WEAK sub-ensemble, with + 0.40. There is also a weak but significant negative correlation between the ENSO-amplitude change and the Q net feedback change (Fig. 12b). This is a hint that the Q net feedback is not a driver of ENSO-amplitude change, but reacts to ENSO-amplitude change, as an increase in ENSO amplitude tends to go along with a stronger increase in Q net damping. Further, we find some indications that the changes in the atmospheric mean state seem to play a role in ENSO  (Fig. 12c) and the changes in the mean vertical wind in the Niño5 region (Fig. 12d). The ENSO amplitude tends to increase when the mean ascending in the Niño5 region and the mean U in northern CEP get weaker, which indicates an eastward shift of the PWC, consistent with U feedback change shown in Fig. 12a). The relation between the ENSO amplitude change and the atmospheric mean state changes is considerably stronger in STRONG than in WEAK (Fig. 12c, d).
Finally, regarding the non-linearity of the atmospheric response, we find indications that the ENSO-amplitude change in the Niño3 region is correlated with the change in precipitation in the same region during EP El Niño events, which is considerably stronger in the STRONG subensemble, with 0.72, than in the WEAK sub-ensemble, with 0.35, or all models, with 0.51 (Fig. 12f). EP El Niño events are defined here according to the Trans-Niño-Index (TNI), measuring the normalized SSTA difference between the Niño1+2 region and the Niño4 region (Trenberth and Stepaniak 2001). The enhanced precipitation in the Niño3 region during EP El Niño events indicates if the meridional temperature gradient in the Niño3 region becomes sufficiently weak for the Intertropical Convergence Zone (ITCZ) to move to the equator. This was observed, for example, during the strong EP El Niño events of 1982 and 1997, and this has been suggested to having contributed to their remarkable strength. Under present-day conditions, climate models with weaker ENAF tend to underestimate this feature more than climate models with strong ENAF (Fig. 12e). Further, the models of STRONG category tend to more properly simulate the asymmetry of ENSO under present day conditions (Fig. 2a), and therefore also the difference between EP El Niño and CP-El Niño events. This may explain the relation between the ENSO-amplitude change and precipitation change during EP El Niño events in these models.
In summary, our analyses suggest that ENSO-amplitude change depends on more than one factor and on changes in the non-linearity of ENSO. The climate models' ENSO amplitude under present day conditions also seems to depend on more than one factor (Guilyardi et al. 2020); for example, there is no significant correlation with the U feedback Fig. 9 Same as Fig. 8, but here for 21C in STRONG and WEAK sub-ensemble strength (not shown). Moreover, the climate models of the STRONG sub-ensemble tend to realistically simulate the ENSO dynamics and important ENSO properties under present day conditions. This suggests that predictions with these models could be more trustful than the models of WEAK category, as also hypothesized by recent studies (Cai et al. 2020a(Cai et al. , 2021Hayashi et al. 2020).

Summary and discussion
Here we investigated how and why the ENSO atmospheric feedbacks (ENAF) change in CMIP5 and CMIP6 models under the RCP8.5 and SSP5-8.5 scenarios, respectively. We find that most models predict a strengthening of both the U feedback in the Niño4 region and the Q net feedback Fig. 10 Same as Fig. 9, but here shows the shading the difference 21C − 20C; contour lines show the mean state in 20C as shown in Fig. 8, with the thick line as zero line and dashed (solid) lines show negative (positive) values with an increment of ± 0 .75 Pa/s/K, ± 0.3 m/s/K, ± 3 W/ m 2 /K, ± 4 W/m 2 /K and ± 3 W/ m 2 /K for Ω, U, Q net , SW and LH feedback, respectively; Regions with unsignificant values as indicated by a T-test on a 90% confidence level are indicated by stippling (c) Fig. 11 Estimated interplay of ocean dynamics and thermodynamics in the Nino3.4 region, shown here as normalized SST change by Ocean dynamics, Q net , SW, LW, SH and LH feedbacks per 1 K SST change during ENSO growth phase in ERA Interim, STRONG and WEAK sub-ensembles, in a in 20C, in b in 21C and in c for the difference 21C − 20C; the errorbars indicate the uncertainties estimated by a bootstrapping approach as described in the methods section in the combined Niño3/Niño4 region, in agreement with Dommenget and Vijayeta (2019) and Cai et al. (2020b). The strengthening of both feedbacks under global warming can be linked to the different warming of the relative SST (rSST) in the Niño3.4 region ( Fig. 1e-g), as a stronger relative warming in this region tends to shift the rising branch of the Pacific Walker Circulation (PWC) to the east (Fig. 5). Together with the stronger equatorial warming in comparison to the off-equatorial regions this leads to stronger mean ascent and precipitation in the CEP and EEP (Fig. 4c, d) that Fig. 12 For observations/reanalysis data and the individual CMIP5 and CMIP6 models, in a change (21C − 20C) in standard deviation of SST anomalies in the Nino3.4 region on the x-axis vs. the change in wind-SST feedback in the Nino3 region minus Nino5 region (130°E-160°E, 0°-10°N) on the y-axis; b same as a, but here the Q net -SST feedback in combined Nino3 and Nino4 region on the y-axis; c same as a, but here for the change in mean U in the northern equatorial central Pacific (180°-150°W, 0°-10°N) on the y-axis; d same as a, but here the change in mean vertical wind at 500 hPa in the Nino5 region on the y-axis; e meridional temperature gradient between 2.5°S-2.5°N and 5°N-10°N, both 150°W-90°W, on the x-axis vs. precipitation during EP El Nino events in the Nino3 region in 20C; f change (21C − 20C) in standard deviation of SST anomalies in Nino3 on the x-axis vs. change in precipitation in the Nino3 region during EP El Nino events on y-axis; the regression lines and correlations are given for all models in black; in a-d and f additionally for STRONG subensemble in red and for WEAK sub-ensemble in blue; three, two or one stars behind the correlation value indicate a significant correlation on a 99%, 95% and 90% confidence level increases the vertical wind response during ENSO events over that region (Fig. 6a-f). As both the U feedback and the Q net feedback are strongly correlated with the vertical wind response, the increase of the latter amplifies the two feedbacks (Figs. 6g-i and 7a-c). Further, the change in the Q net feedback is dominated by changes in the SW feedback in the combined Niño3/Niño4 region. The SW feedback also becomes stronger due to the increased convective response under a warming climate (Fig. 7d-f). The change in the LH feedback contributes to the Q net feedback change in the Niño3 region ( Fig. 7g-i), but to a lesser extent than the SW feedback. The processes explaining the changes of the ENAF under global warming are similar to that in Bayr et al. (2020), investigating the ENAF-strength changes between uncoupled atmosphere model (AMIP) simulations and coupled (CMIP) simulations under present-day conditions. We separated the models into two sub-ensembles according to their ENAF strength. Due to more realistic ENAF and ENSO dynamics the models of the STRONG sub-ensemble tend to more realistically simulate important ENSO properties under present day conditions than the models of WEAK category, for example, the ENSO asymmetry and ENSO phase locking (Fig. 2). Under global warming, first, the biased mean state in WEAK in 20C constrains the ENAF projection in 21C, as the models of the WEAK subensemble also have weaker ENAF in 21C (Figs. 1b,9). This results in considerably different ENSO dynamics in WEAK than in STRONG in 21C (Fig. 11b). Therefore the models in the WEAK category tend to have a smaller ENSO asymmetry in 21C in comparison to the models in the STRONG sub-ensemble. The reduced ENSO asymmetry in both periods in WEAK suggests that due to their overly strong cold tongue, these models have problems in simulating the strong EP El Nino events with a large precipitation response over the EEP (Fig. 12e), in agreement with previous studies (Cai et al. 2015(Cai et al. , 2020aBayr et al. 2021). Second, under global warming the models of the STRONG sub-ensemble tend to simulate a stronger reduction of the SST gradient across the equator than the models of the WEAK sub-ensemble (Fig. 3). Despite these differences, the atmospheric meanstate changes in the two sub-ensembles are similar and describe an El Niño-like change, as also reported in previous studies (Cai et al. 2015(Cai et al. , 2020bLian et al. 2018). Nevertheless, we postulate that one should have more confidence in the ENSO projections of the models in the STRONG subensemble, as these exhibit the more realistic mean state, ENAF and ENSO dynamics in 20C, which results in more realistic ENSO asymmetry, phase locking and teleconnections. Especially a realistic ENSO asymmetry was found to be very important to reduce uncertainty in ENSO projections (Cai et al. 2020a(Cai et al. , b, 2021Hayashi et al. 2020). Finally, there is not much difference in the ENSO-amplitude and ENAF changes between CMIP5 and CMIP6, in agreement with previous studies (Planton et al. 2020;Beobide Arsuaga et al. 2021).
There is an ongoing debate how ENSO amplitude will change under global warming (Van Oldenborgh et al. 2005;Meehl et al. 2007;Latif and Keenlyside 2009;DiNezio et al. 2012;Stocker et al. 2013;Kim et al. 2014a;Cai et al. 2015;Chen et al. 2017;Beobide Arsuaga et al. 2021). We find that ENSO amplitude tends to get larger, when the U feedback increases in the EEP and decreases in the far WEP, thus shifting eastward (Fig. 12a). However, this can only explain partly the ENSO-amplitude change. The Q net feedback exhibits a small negative correlation with the ENSOamplitude change (Fig. 12b). This suggests that the Q net feedback change is not a driver but reacts to the ENSOamplitude changes. The amplification of both feedbacks under global warming in most models may also contribute to the unclear signal in ENSO-amplitude change, as the simultaneous increase of an amplifying and damping feedback can cause compensating effects (Lian et al. 2018;Dommenget and Vijayeta 2019;Cai et al. 2020bCai et al. , 2021. Finally, we find some hints that the change in precipitation in the Niño3 region during EP El Niño events could explain the ENSO-amplitude change in the Niño3 region in models of the STRONG sub-ensemble. This provides some evidence for non-linear dynamics being important to explain ENSOamplitude change under global warming, which are better simulated in models with strong ENAF under present-day conditions. In summary, we find some evidence that ENSO amplitude change under global warming is influenced by multiple changes in the atmospheric mean state and feedbacks (Fig. 12). This should go along with changes in the oceanic mean state and feedbacks that are discussed in literature (Kim et al. 2014a;Cai et al. 2015Cai et al. , 2020bCai et al. , 2021.