1 Introduction

The El Niño/Southern Oscillation (ENSO) is the dominant mode of tropical climate variability on interannual timescales and driven by a complex interplay between amplifying (positive) and damping (negative) coupled ocean-atmosphere feedbacks (e.g. Jin et al. 2006). ENSO is associated with extreme weather such as heavy precipitation events and droughts in the tropical Pacific region and beyond (e.g. Philander 1990; Yeh et al. 2018). During its warm (cold) phase, El Niño (La Niña), the sea surface temperature (SST) is warmer (colder) than normal in the eastern and central equatorial Pacific. These warm (cold) SST anomalies are related to anomalously high (low) upper-ocean heat content and eastward (westward) advection of warm (cold) water. The onset of El Niño events are often triggered by westerly wind bursts over the western equatorial Pacific that force eastward propagating downwelling oceanic Kelvin waves that warm the SST in the central and eastern equatorial Pacific, where the thermocline is shallow (e.g. Lengaigne et al. 2004; Fedorov et al. 2015; Neske and McGregor 2018; Timmermann et al. 2018). This initial warming is amplified by the Bjerknes feedback loop, i.e., the equatorial SST warming in the east weakens the zonal wind over the western Pacific, which then results in a deeper thermocline that reinforces the initial positive SST anomaly in the central and especially eastern equatorial Pacific.

The growth of SST anomalies is damped by the heat flux feedback consisting mainly of the negative shortwave-SST feedback and the negative latent heat flux-SST feedback (Lloyd et al. 2009). The negative shortwave-SST feedback is caused by the eastward shift of atmospheric deep convection during an El Niño event and dominates the heat flux feedback over the western and central equatorial Pacific, where the largest shift in convection is observed (Lloyd et al. 2011). The negative latent heat flux-SST feedback relates to an increased evaporation due to warmer SSTs during El Niño events and is strongest over the eastern equatorial Pacific, where the SST change is largest (Lloyd et al. 2011). The same amplifying and damping feedbacks operate during La Niña events, but the anomalies have opposite signs. In association with the longitudinal shift of the rising branch of the Walker Circulation over the western and central Pacific, there is a pronounced asymmetry in the shortwave feedback between El Niño and La Niña (Lloyd et al. 2012; Bellenger et al. 2014; Bayr et al. 2018). The shortwave feedback varies both in amplitude and location, thus is stronger and operates more eastward during El Niño than during La Niña (Bayr et al. 2018). The latent heat flux feedback over the eastern Pacific is mostly linear (Lloyd et al. 2011).

ENSO events typically grow during boreal summer and autumn, reach their maximum around the end of the year and decay in the following boreal spring (Timmermann et al. 2018). This seasonal phase locking of ENSO can be explained by the seasonality of the amplifying and damping feedbacks. During ENSO growth in boreal summer and autumn, the positive coupled feedbacks reach their maximum while the thermodynamic damping is at its minimum. This can be explained by the seasonal outcropping of the thermocline in the east (Galanti et al. 2002), the more equatorward position of the Intertropical Convergence Zone (Tziperman et al. 1998) and the reduction in the negative cloud feedbacks (Dommenget and Yu 2016). During the ENSO decay phase in the following spring, the thermodynamic damping reaches its maximum while the positive coupled feedbacks are at their minimum (e.g. Tziperman et al. 1998; Galanti et al. 2002; Dommenget and Yu 2016; Wengel et al. 2018).

One important aspect of ENSO is the asymmetry between El Niño and La Niña, i.e., that SST anomalies during El Niño events are usually stronger and located further to the east than during La Niña events (Takahashi et al. 2011; Dommenget et al. 2013; Capotondi et al. 2014; Timmermann et al. 2018). The causes of this asymmetry are still under debate. Important contributors are: nonlinearities in the wind-SST feedback (Frauen and Dommenget 2010; Karamperidou et al. 2017) and in the shortwave-SST feedback (Lloyd et al. 2009, 2011, 2012; Bellenger et al. 2014), stronger subsurface-surface coupling during El Niño (Meinen et al. 2000; Lübbecke and McPhaden 2017), larger dynamical ocean response per unit zonal wind stress change during El Niño (Im et al. 2015), nonlinear time-averaged ocean response to tropical instability waves (Imada and Kimoto 2012), and state-dependent stochastic noise (Levine et al. 2016; Hayashi and Watanabe 2017). Teleconnections from the tropical Atlantic and Indian Oceans (Okumura and Deser 2010; An and Kim 2018) and the extratropics also may contribute to ENSO asymmetry (Li et al. 2007; Wu et al. 2010).

State-of-the-art climate models still exhibit severe deficits in simulating important ENSO properties such as the seasonal phase locking (Bellenger et al. 2014; Wengel et al. 2018) or the asymmetry between El Niño and La Niña (Dommenget et al. 2013; Zhang and Sun 2014; Timmermann et al. 2018). These deficiencies can be partly explained by a cold equatorial SST bias present in many climate models, which shifts the models into a La Niña-like mean state with a too western position of the rising branch of the Walker Circulation. The wind-SST feedback and heat flux-SST feedback are strongly underestimated in these models (Kim et al. 2014; Wengel et al. 2018; Bayr et al. 2018). As the position of the rising branch of the Walker Circulation determines the strength of both atmospheric feedbacks (Bayr et al. 2020), models that underestimate the wind-SST feedback also tend to underestimate the heat flux-SST feedback, with error compensation between these two feedbacks (Guilyardi et al. 2009; Bellenger et al. 2014; Kim et al. 2014; Bayr et al. 2018, 2019). Bayr et al. (2019) investigated in detail how this error compensation hampers the simulated ENSO dynamics in climate models: they analyzed within the Bjerknes Stability (BJ) index framework (Jin et al. 2006) the positive thermocline, zonal advection and Ekman feedbacks and the negative thermodynamic damping. All four feedbacks are weaker if the atmospheric feedbacks are underestimated, but the resulting total BJ index is not too different from that in models with strong atmospheric feedbacks. To quantify the error compensation, they also calculated the contribution of the atmospheric heat fluxes to the SST tendency during the ENSO growth phase. They report that in many climate models ENSO is a hybrid of ocean-driven and shortwave-driven ENSO dynamics, in equal proportions in the most biased models. This already was suggested by previous studies (Dommenget 2010; Dommenget et al. 2014; Bayr et al. 2018). Bayr et al. (2019) quantified for the first time the contribution of the shortwave feedback to the SST tendency.

This study is part of a series of studies that focus on the ENSO atmospheric feedbacks: Bayr et al. (2018) described the relation between the equatorial SST bias, the mean-state Walker Circulation, and the ENSO atmospheric feedback strength. Bayr et al. (2019) revealed the effect of the error compensation of the atmospheric feedbacks on simulated ENSO dynamics and the hybrid nature of ENSO dynamics in biased models. Finally, Bayr et al. (2020) showed that also in uncoupled models the atmospheric mean state is important for atmospheric feedback strength and explained the difference in atmospheric feedback strength between uncoupled and coupled simulations.

The aim of this study is to better understand the interplay of thermodynamics and ocean dynamics during the ENSO growth phase in observations, reanalysis products and in climate models participating in the Coupled Model Intercomparison Project phase 5 (CMIP5). We employ the offline slab ocean SST method introduced in Bayr et al. (2019) and apply it more accurately by using a temporally and spatially varying mixed layer depth (MLD) instead of a constant MLD. This enables a better understanding of the role of the heat fluxes in ENSO dynamics, as the MLD considerably varies with the seasons and during ENSO events, especially in the eastern equatorial Pacific where the MLD is shallow and therefore has a stronger influence on the SST. Further, the importance of the atmospheric feedbacks for the seasonal phase locking and the asymmetry of ENSO is investigated.

This study is organized as follows: in Sect. 2, we describe the data and methods. Section 3 provides the results pertaining to the interplay of thermodynamics and ocean dynamics during the ENSO growth phase in the different datasets. In Sect. 4, we investigate the seasonal ENSO phase locking and in Sect. 5, the ENSO asymmetry. The results are summarized and discussed in Sect. 6.

2 Data and methods

We use observed SSTs for the period 1958–2018 from HadISST (Rayner et al. 2003), zonal wind stress and heat fluxes for the period 1958–2001 from ERA40 (Uppala et al. 2005) and for the period 1979–2018 from ERA-Interim reanalysis (Simmons et al. 2007). Observed heat fluxes for the period 1979–2018 are from the Wood Hole Oceanographic Institute data set, also referred to as OA Flux data set (Yu et al. 2008). OA Flux uses shortwave and longwave fluxes from ISCCP (Rossow and Schiffer 1999) for the period 1979–2001 and from CERRES (Wielicki et al. 1996) for the period 2002–2018. The mixed layer depth (MLD) is calculated from ocean temperature from SODA reanalysis (Carton and Giese 2008) for the period 1958–2001 and from observations using HadEN4.2 (Good et al. 2013) for the period 1979–2018. The MLD is defined as the depth at which the ocean temperature has decreased by 0.2 K compared to the surface. Further, we set the minimum MLD to 10 m, as with the 0.2 K criterion some models have MLDs of even less than 1 m. As we divide by the MLD to compute the offline slab ocean temperature, such small MLDs would give unrealistically large temperature changes. The 10 m limit can be justified by the fact that observed MLDs in the tropical Pacific less than 10 m are very rare. Further, most ocean components of the CMIP5 models only have one vertical level in the upper 10 m and the depth of the first layer determines the possible minimum MLD. Finally, the shortwave radiation warms the upper 50 m of the ocean independent of the MLD, as e.g. two thirds of the surface solar irradiance penetrates below 10 m (Hieronymi et al. 2012).

From CMIP5, we analyze the twentieth century historical simulation (1900–1999) for a set of 38 models. Only those models are considered for which the heat fluxes, zonal wind stress and ocean temperatures are available (see Fig. 3 for a list of models). The MLD was calculated in the same way as for observations using the 0.2 K difference criterion. All model data is interpolated on a regular 2.5° × 2.5° grid. We use monthly-mean anomalies relative to seasonally varying climatology and the data from each calendar month is detrended separately.

ENSO events are defined according to Trenberth (1997) as time periods in which the five-month running mean of Niño3.4 SST is above the ± 0.5 standard deviation threshold for at least 6 consecutive months. We determine the time of the maximum for each event separately (Fig. 1a), and long ENSO events like the 1998/1999 La Niña with two maxima are considered as two events if there are more than 11 months between the two peaks. The Niño3 region is defined 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, and the Niño4 region as 160°E–150°W and 5°S–5°N.

Fig. 1
figure 1

a Five month running mean of SST anomalies in the Niño3.4 region from 1979 to 2018 in HadISST; the black horizontal lines mark the threshold for ENSO events, which has to be passed for at least 6 consecutive months; green circles mark the beginning of an ENSO event, yellow the growth period, red the maximum, cyan the decay phase and blue the end; b in black the SST evolution (dSST) in Niño3.4 region averaged over all El Niño events from 6 months before the maximum till the maximum, normalized by the dSST at the maximum; in gray the evolution of the Offline Slab Ocean SST (dSSTslab), calculated by the integration of the net heat flux anomalies to the local mixed layer depth, as described in Eq. (1); the blue arrow indicates the SST change by ocean dynamics (dSSToc), estimated from the difference between dSST and dSSTslab in Eq. (3)

The seasonal phase locking of ENSO is measured by the phase locking index proposed by Bellenger et al. (2014). It is the ratio between the Niño3.4-averaged standard deviation of the SST in the months November, December, January (NDJ, high variability months) divided by the average standard deviation in the months April, May, June (AMJ, low variability months). Phase locking of the net heat flux, shortwave and latent heat flux feedback is defined as the difference between the average feedback in February, March, April (FMA, strongest damping) minus the average feedback in September, October, November (SON, weakest damping) as shown in Sect. 4.

The thermodynamic contribution to the SST change during ENSO growth is calculated by integrating the net heat flux Qnet, similar to Drews and Greatbatch (2016). Here, the integration is performed over the 6 months preceding the maximum of each ENSO event, as 6 months is the average growth length in observations (Fig. 1a, b):

$$dSST_{slab} = \frac{1}{{c_{p} \cdot \rho }}\int_{t = - 6}^{t = 0} {\frac{{Q_{net} }}{H}dt}$$
(1)

where cp= 4000 J kg−1 K−1 is the specific heat capacity at constant pressure of sea water, ρ = 1024 kg m−3 the average density of sea water, H the temporally and spatially varying MLD in meters, and t time in months. We normalize the SST change (dSST) and the SST change caused by the heat fluxes (dSSTslab) by dSST at the maximum (t = 0), yielding an SST change per Kelvin warming (Fig. 1b) and enabling comparison of the contributions in the individual climate models irrespective of their ENSO amplitude. We also perform the integration with the shortwave (SW), longwave (LW), sensible (SH) and latent heat (LH) flux in order to estimate their contribution to the SST tendency. As a simplification, we integrate the full SW anomalies to the dSSTslab even if the MLD is shallower than the SW penetration depth. The effect of this simplification is negligible, as such shallow MLDs are quite rare and the heat flux feedback in these regions (mainly the eastern equatorial Pacific) is dominated by the LH feedback. Finally, as the observed SST change is the balance of the SST changes due to ocean dynamics (dSSToc) and net heat flux (dSSTslab), we can estimate the contribution of the ocean dynamics from the results above (Fig. 1b):

$$dSST \approx dSST_{oc} + dSST_{slab}$$
(2)
$$dSST_{oc} \approx dSST - dSST_{slab}$$
(3)

In this way, it is possible to quantify the contributions of thermodynamics and ocean dynamics during ENSO growth and to describe their interplay.

To estimate the uncertainty in the SST change by the heat fluxes, we apply a bootstrapping approach. We randomly choose 1000 times only two thirds of the ENSO events and calculate the SST change. We show the uncertainty as error bars, indicating the 90% quantile of the estimated values. In a similar manner, we estimate the uncertainty of the seasonal variation of the SST variability and the mean MLD, but here by randomly choosing 1000 times two thirds of the years.

3 ENSO growth in observations and CMIP5 models

We apply Eq. (1) to the observations and reanalysis products over the Niño3.4 region. This yields a damping of − 1.6 ± 0.2 K by the net heat flux per Kelvin change in SST (dSSTslab(Qnet), Fig. 2). The damping is primarily caused by a negative shortwave (dSSTslab(SW), − 0.9 ± 0.1 K/K) and latent heat flux-SST feedback (dSSTslab(LH), − 0.7 ± 0.1 K/K), while the contributions of the longwave radiation and sensible heat flux are small (Fig. 2). According to Eq. (3), the ocean dynamics (dSSToc) would yield an amplifying contribution of about + 2.6 K ± 0.2 K per 1 K SST change, as the rest is damped away by heat fluxes (Fig. 2). The error bars in Fig. 2 indicate the uncertainty of each feedback calculated by the bootstrapping approach described above.

Fig. 2
figure 2

Average SST change in the Niño3.4 region at the maximum of all ENSO events in an observation and reanalysis ensemble, due to ocean dynamics (blue bar), net heat flux (gray bar), shortwave radiation (yellow bar), longwave radiation (red bar), sensible heat flux (cyan bar) and latent heat flux (green bar); All values are normalized per 1 K total SST change in the Niño3.4 region; The errorbars indicate the 90% confidence interval estimated by a bootstrapping approach, as described in the methods section

An analysis similar to that shown in Fig. 2 for observations is performed for the CMIP5 models. In agreement with previous studies by Lloyd et al. (2009, 2011, 2012) and Bayr et al. (2018, 2019), we find that in the Niño3.4 region differences in the net heat flux (Qnet) feedback among the CMIP models is mostly due to differences in the shortwave (SW) feedback (Fig. 3a). Despite a marked spread in MLD in the Niño3.4 region (Fig. 3c), we find a strong correlation between dSSToc and dSSTslab(SW) in the models, amounting to − 0.73 (Fig. 3b). This strong correlation highlights the error compensation between the underestimated wind-SST forcing and the underestimated heat flux-SST damping in many CMIP5 models (Fig. 3d). Bayr et al. (2019) showed that models with strongly underestimated wind-SST feedback can have ENSO statistics that at first glance are not too different from observations but due to very different ENSO dynamics. These models exhibit hybrid wind-driven and shortwave-driven ENSO dynamics, as the weaker ocean dynamical forcing is compensated by an erroneously positive SW feedback that acts as an additional forcing (Fig. 3b). In the most strongly biased models, the SW feedback contributes to the SST change with a similar magnitude as the ocean dynamics (Fig. 3b). The strength of dSSTslab(SW) is dominated by the strength of the SW feedback, with a correlation of 0.93 (not shown), which explains why dSSTslab(SW) and MLD are uncorrelated (Fig. 3c). The three outliers in Fig. 3b (model number 14, 15 and 33) have a stronger than average latent heat flux (LH) feedback and shallower than average MLD. Therefore, the dSSToc is stronger in these models than in models with a similar dSSTslab(SW), due to a stronger dSSTslab(LH). We find that the model spread in the LH feedback over the Niño3.4 region is generally much smaller than the spread in SW feedback (10 W/m2 in the former in comparison to 30 W/m2 in the latter) and that there is no significant correlation between the LH feedback and the Qnet feedback (not shown). Thus, the error compensation between the SW feedback and the wind feedback is very clear in the Niño3.4 region during the ENSO growth phase, even though the models differ substantially in MLD.

Fig. 3
figure 3

For observations/reanalysis and CMIP5 models in a net heat flux feedback in the Niño3.4 region on the y axis vs. the shortwave feedback in Niño3.4 region on the x axis; b SST change during ENSO events by ocean dynamics on the y axis vs. SST change by shortwave radiation on the x axis, both for the Niño3.4 region; c same as b, but here the mean mixed layer depth during ENSO growth phase on the y axis; d same as a, but here for the zonal wind stress-SST feedback in the Niño4 region on the x axis; e average strength of wind and heat flux feedback, both normalized by the observed strength, on the x axis vs. relative SST bias in the Niño4 region on the y axis; The black line is the regression line between the CMIP5 models with the correlation given in the upper right/left corner, where 1, 2, 3 stars are indicating a significant correlation on a 90%, 95% and 99% confidence level, respectively; the color of the numbers indicate the sub-ensembles with STRONG (red) and WEAK (blue) ENSO atmospheric feedbacks

To elucidate the effect of the error compensation on the simulated SST change in the Niño3.4 region during the ENSO growth phase, we define two sub-ensembles: STRONG (indicated in red in Fig. 3) and WEAK (indicated in blue). The sub-ensembles are defined according to their combined atmospheric feedback strength, with STRONG > 0.55 and WEAK < 0.55 of the observed atmospheric feedback strength. The combined atmospheric feedback strength is defined as the average of the wind and heat flux feedback, both normalized by the observed value. We use the combined atmospheric feedback strength for defining the sub-ensembles, as it is a measure of the error compensation between the two feedbacks and enables separating the models into sub-ensembles in which dSSTslab(SW) acts either as damping (STRONG sub-ensemble) or erroneously as forcing (WEAK sub-ensemble) (Fig. 3b). The atmospheric feedback strength depends on the relative SST bias in the Niño4 region (Fig. 3e), as described in more detail in Bayr et al. (2018). The relative SST is defined with respect to the tropical Pacific (120° E–80°W, 30°S–30°N) area-mean SST and used to account for the different mean temperatures in the climate models, as tropical convection is more strongly related to relative SST than to absolute SST (Johnson and Xie 2010; Bayr and Dommenget 2013; Izumo et al. 2019).

Compared to observations, the STRONG sub-ensemble has quite realistic ENSO dynamics, with only 0.5 K/K weaker dSSToc and dSSTslab(Qnet) (Fig. 4). From STRONG to WEAK, dSSTslab(Qnet) decreases from − 1.1 K/K to − 0.2 K/K, with dSSTslab(SW) contributing most to the difference between the two sub-ensembles. The dSSToc decreases similarly from + 2.1 K/K to + 1.2 K/K, while dSSTslab(SW) changes from a − 0.8 K/K damping to a + 0.6 K/K forcing, causing the hybrid nature of ENSO dynamics in models belonging to the WEAK sub-ensemble, with a large shortwave-driven component. The dSSTslab(LH) has a quite similar strength in both CMIP5 sub-ensembles (− 0.5 K/K and − 0.6 K/K in STRONG and WEAK, respectively) than in observations (− 0.7 K/K).

Fig. 4
figure 4

Same as Fig. 2, but here additionally for sub-ensembles of CMIP5 models with STRONG and WEAK atmospheric feedbacks, as indicated by the colors in Fig. 3); The values shown here are the averages over the boxes shown in Fig. 5; The errorbars indicate the 90% confidence interval estimated by a bootstrapping approach, as described in the methods section

Performing the analysis at each grid point provides the patterns (Fig. 5) that are shown at the peak of the ENSO event (t = 0). As we will first focus on the strength of the dynamics in observations and CMIP5 sub-ensembles and later on the asymmetry between El Niño and La Niña, we show here the results for both El Niño and La Niña events together. Further, as the dSSTslab (in K) depends on both the heat flux strength (in W/m2) and the MLD, we show additionally the equatorial heat flux strength for Qnet, SW and LH (in W/m2) and the equatorial MLD in Fig. 6. Thus, for clarification, when we refer, for example, to the Qnet feedback, the feedback in W/m2 is meant, while dSSTslab(Qnet) is the effect of the heat flux feedback on SST in K, with temporarily and spatially varying MLD considered.

Fig. 5
figure 5

a observed SST change [dSST] at the maximum of the ENSO event; b SST change by ocean dynamics [dSSToc], estimated as residual in Eq. (3); c SST change by the net heat flux [dSSTslab(Qnet)], as calculated in Eq. (1); d same as c but here for the shortwave radiation [dSSTslab(SW)]; e same as c but here for latent heat flux [dSSTslab(LH)]; fj same as ae, but here for the CMIP5 STRONG sub-ensemble; ko same as ae but here for the CMIP5 WEAK sub-ensemble; Note that all figures are normalized by dSST in Niño3.4 (black box) for a better comparison; pt same as ae, but here the difference between STRONG and WEAK sub-ensembles (second minus third row); The hatching in pt indicates significant differences on a 90% confidence level calculated by a Student’s t test

Fig. 6
figure 6

Heat flux feedbacks and mean mixed layer depth along the equator (5° S–5° N), in a net heat flux regressed on SST in Niño3.4 region in observation/reanalysis data; b same as a but here for the STRONG sub-ensemble; c same as a but here for the WEAK sub-ensemble; df same as ac but here for the shortwave radiation; gi same as ac but here for the latent heat flux feedback; jl same as ac but here for mean equatorial mixed layer depth; The solid line mark the values for all months, the dashed line for El Niño months and dashed-dotted for La Niña months; the light gray area marks the spread between the different models for “all months” case; The dark gray solid line in the second and third column are the observed values for the “all months” case for comparison

In observations, the SST change in the 6 months before the maximum of the ENSO event (dSST, Fig. 5a) exhibits its maximum to the west of the largest SST anomaly (not shown). The dSSTslab(Qnet) has its maximum over the eastern Pacific (Fig. 5c) and its pattern is similar to dSSToc (Fig. 5b). Largest dSST is in Niño3.4, because there the ocean dynamical forcing is strong while heat flux damping is weak. The strongest heat flux damping is observed over the eastern part of the basin. The dSSTslab(LH) dominates the damping in the eastern part of the Niño3 region (90°W–120°W, 5°S–5°N, Fig. 5e), while dSSTslab(SW) dominates the damping in the Niño3.4 and Niño4 regions (Fig. 5d). The dSSTslab(Qnet) (Fig. 5c) is much smaller in the western than in the eastern Pacific even though the SW feedback in the western Pacific (≈ − 20 W/m2/K, Fig. 6d) is stronger than the LH feedback over the eastern Pacific (≈ − 13 W/m2/K, Fig. 6g). This can be explained by the shallower mixed layer in the eastern (20–30 m, Fig. 6j) relative to that in the central and western Pacific (40–60 m, Fig. 6j). Therefore, the minimum of dSSTslab(SW) in the western Pacific, amounting to − 2 K/K, is half as strong as the minimum in dSSTslab(LH) in the eastern Pacific (Fig. 5d, e). The LW radiation is mostly of opposite sign with regard to the SW radiation but much weaker, and the SH contributes very little to dSSTslab(Qnet) (not shown).

We conclude from the above analyses: first, the Qnet feedback is largest over the Niño3.4 region (− 20 W/m2) due to strong SW and LH damping (Fig. 6a, d, g), but dSSTslab(Qnet) is relatively weak there due to the deep MLD (Figs. 5c, 6j). Second, while dSST is largest in the Niño3.4 region, dSSToc is not (Fig. 5a, b). The situation is quite different in the eastern part of the Niño3 region: the Qnet damping is weaker (between − 20 and − 10 W/m2/K, Fig. 6a), but due to a much shallower MLD (Fig. 6j), the dSSTslab(Qnet) has its maximum in the east (Fig. 5c). The dSSToc also has its maximum in the eastern part of the Niño3 region (Fig. 5b). However, dSST is much smaller there than in the Niño3.4 region due to stronger dSSTslab(Qnet) damping in the Niño3 region (Fig. 5a, c). In summary, the interplay between thermodynamics and ocean dynamics differs substantially between the Niño3.4 region and eastern part of the Niño3 region, and the difference in the MLD between these two regions is an important factor for this difference.

The two CMIP5 sub-ensembles on average simulate the ensemble-mean MLD quite realistically but exhibit a large spread between the models (grey shading in Fig. 6k, l) as already indicated in Fig. 3c). The spread in MLD is largest over the central Pacific and slightly larger in STRONG than in WEAK. The dSST is extending further to the west in WEAK than in STRONG (Fig. 5f, k), which becomes obvious in the difference between STRONG and WEAK (Fig. 5p). This more westward extension in the dSST pattern in WEAK in comparison to STRONG can be related to a decrease of dSSTslab(Qnet) over the central Pacific, which almost vanishes in WEAK (Fig. 5h, m, r). The pattern of the estimated dSSToc (Fig. 5g, l, q) is quite realistic in STRONG but much weaker in WEAK. In WEAK, there is a secondary maximum in the west (Fig. 5l).

The differences in dSSTslab(Qnet) mainly result from dSSTslab(SW) (Fig. 5i, n, s), as its minimum is much further west in WEAK than in STRONG. Further, an area with a positive SST change emerges in the eastern Pacific in WEAK. The dSSTslab(LH) is a bit stronger in WEAK than in STRONG along the equator (Fig. 5j, o, t), but the difference in dSSTslab(Qnet) is clearly dominated by the difference in dSSTslab(SW) (Fig. 5r, s, t). When we compare the heat flux feedbacks and MLDs in the two sub-ensembles (second and third column in Fig. 6), we find that the difference in dSSTslab(SW) is mainly caused by the SW feedback (Fig. 6e, f), as the mean MLD is quite similar in both sub-ensembles (Fig. 6k, l). Further, the LH feedback is quite similar in the two sub-ensembles (Fig. 6h, i), so that the difference in the Qnet feedback is dominated by the difference in SW feedback (Fig. 6b, c). Finally, the changes in the SW feedback can be linked to a westward shift of the Pacific Walker Circulation from STRONG to WEAK due to a larger equatorial cold SST bias (Fig. 3e), which leads to the positive SW feedback over the eastern Pacific in WEAK (Bayr et al. 2018, 2020). The influence of the error compensation between the underestimated wind and heat flux feedbacks on ENSO dynamics is largest in the Niño3.4 region, where the SST (in observations) is less damped in comparison to the eastern part of the Niño3 region and the largest differences in dSSTslab(SW) and dSSToc occur. Therefore, only a weak dSSToc forcing is required in models of the WEAK category to produce a sizeable SST anomaly.

In summary, ENSO dynamics in the Niño3.4 region differ substantially among the CMIP5 models, with only half of the models simulating a realistic interplay between thermodynamic and ocean dynamical feedbacks during the ENSO growth phase, while the other half exhibits a hybrid of ocean-driven and shortwave-driven ENSO dynamics due to an equatorial cold SST bias (Fig. 3). Further, the too westward extension of the ENSO-related SST anomaly pattern observed in many models can be explained by the too weak heat flux damping and too strong ocean dynamical forcing over the central and western equatorial Pacific.

4 Seasonal ENSO phase locking

The realistic simulation of the seasonal ENSO phase locking is still a challenge (Bellenger et al. 2014; Wengel et al. 2018). In observations, there is a distinct peak in the occurrence of maximum SST anomalies during ENSO events in December (Fig. 7a). In the CMIP5 models the distribution is much more spread over the calendar year (Fig. 7b, c). In STRONG, the probability of ENSO peaking in boreal spring is quite low but it is considerably higher in WEAK (Fig. 7b, c). The seasonal ENSO phase locking over the Niño3.4 region is simulated more realistically in STRONG than in WEAK (Fig. 8a). In agreement with Wengel et al. (2018), the seasonal phase locking can be explained in observations by the maxima of the wind feedback in boreal autumn (Fig. 8d) and of the heat flux damping in boreal spring (Fig. 8b). The latter is mainly caused by the SW feedback (Fig. 8e), but also to a minor part by the LH feedback (Fig. 8h). Observations (purple dot) and reanalysis data (yellow and orange dot), however, do not agree on the seasonal variation of the LH feedback (Fig. 8i). The seasonal cycle of the MLD enhances the effect of the seasonal cycle of the SW and LH feedbacks, as the mixed layer is shallowest in boreal spring when these heat flux feedbacks are strongest (Fig. 8g).

Fig. 7
figure 7

Probability that ENSO events peak in the individual calendar month, in a for observations, in b for the STRONG sub-ensemble and in c for the WEAK sub-ensemble

Fig. 8
figure 8

a Standard deviation of SST in the Niño3.4 region for each calendar month, normalized by the annual mean standard deviation, for observation, STRONG and WEAK sub-ensemble; b same as a but here net heat flux of Niño3.4 region regressed on SST of Niño3.4 region for each calendar month; c phase locking index of SST in Niño3.4 (defined as mean std in NDJ divided by mean std in AMJ) on the y axis vs. phase locking index of net heat flux in Niño3.4 region (defined as mean std in FMA minus mean std in SON) on the x axis; d same as b but here the zonal wind stress in Niño4 region regressed on SST in Niño3.4 region; e same as b but here for the shortwave radiation in Niño3.4; f same as c but here for phase locking index of the shortwave radiation (same definition as for Qnet) on the x axis; g same as a but here for the mean mixed layer depth in Niño3.4; h same as b but here for the latent heat flux in Niño3.4; i same as c but here for phase locking index of the latent heat flux (same definition as for Qnet) on the x axis; the errorbars mark the uncertainty on a 90% confidence level, in a, g estimated by a bootstrapping approach as described in the methods section and in b, d, e, h estimated by the uncertainty of the regression slope; the colors in c, f, i mark the members of the sub-ensembles with STRONG and WEAK ENSO feedbacks; the yellow, orange and purple circles are the values of ERA Interim, ERA40 and OA Flux, respectively; the correlation between the individual CMIP5 models is given in the upper right corner, with 3 stars indicating a significant correlation on a 99% confidence level

The seasonality of the wind feedback is quite similar in the two sub-ensembles, but the averaged strength is larger for models in STRONG relative to those in WEAK (Fig. 8d). With respect to the Qnet and SW feedbacks, the seasonality is pronounced in STRONG while there is less seasonality in WEAK (Fig. 8b, e). The seasonality of the LH feedback and MLD is similar in both sub-ensembles (Fig. 8g, h). Significant correlations exist in the model ensemble between the seasonal phase locking of SST and that of the Qnet feedback (− 0.71, Fig. 8c), the SW feedback (− 0.59, Fig. 8f), and the LH feedback (− 0.53, Fig. 8i). This suggests that the heat flux damping plays an important role in the ENSO phase locking, as models with a stronger seasonality in the heat flux feedbacks tend to exhibit a stronger SST phase locking. Those are not the only factors that influence the SST phase locking, as some models exhibit weak SST phase locking in the presence of a strong phase locking of the Qnet and SW feedbacks (e.g. models 14 and 22). Further, other models have a reasonable SST phase locking despite weak seasonality in the Qnet, SW and LH feedbacks (e.g. model 21), arguing for an importance of, for example, the wind-SST feedback.

5 ENSO asymmetry

ENSO asymmetry, i.e., the differences between the pattern of El Niño and La Niña events, is investigated next. In observations, the maximum of dSST is more in the east and zonally wider during El Niño than during La Niña (Fig. 9a, f), resulting in a difference with an insignificant positive pole in the east that is surrounded by a significant horseshoe-type pattern of opposite sign (Fig. 9k). To account for the skewness, i.e., the asymmetry in amplitude between El Niño and La Niña events, both patterns have been normalized by the respective dSST averaged over the Niño3.4 region (black box in Fig. 9a, f), which is + 1.17 K and − 0.89 K for El Niño and La Niña, respectively. The dSSToc in the eastern Pacific is larger during La Niña than during El Niño and extends slightly more westward (Fig. 9b, g). The difference between El Niño and La Niña reveals that dSSToc during El Niño is smaller in the western, northwestern and in the far eastern equatorial Pacific, and larger in the central equatorial Pacific (Fig. 9l). The differences in the eastern Pacific, however, are not significant.

Fig. 9
figure 9

Same as Fig. 5a–e, but here separately for El Niño (1. row), La Niña (2. row), and the difference El Niño − La Niña (3. row) in observations/reanalysis; Note that all figures in the 1. and 2. row are normalized by dSST in Niño3.4 (black box) for El Niño and La Niña separately to yield the contribution of each component per 1 K dSST change in Niño3.4; The hatching in the difference plots indicates significant differences on a 90% confidence level calculated by a Student’s t test

The dSSTslab(Qnet) pattern derived from observations is similar in the two ENSO phases (Fig. 9c, h), but with differences in amplitude. The difference pattern reveals a tripole structure at the equator, with stronger dSSTslab(Qnet) damping over the far eastern and western equatorial Pacific during La Niña relative to El Niño and weaker damping in the central Pacific (Fig. 9m). Only the difference in the central Pacific is significant. The dSSTslab(SW) during El Niño affects a larger area and the center extends more eastward than during La Niña (Fig. 9d, i), resulting in a quite strong asymmetry with a dipole structure of the differences at the equator (Fig. 9n). The asymmetry in dSSTslab(SW) is mainly caused by the more eastward position and strengthening of the SW feedback during El Niño and more westward position and weakening during La Niña (dashed and dashed-dotted lines in Fig. 6d, respectively), as the MLD asymmetry would cause a dSSTslab(SW) asymmetry of opposite sign (dashed and dashed-dotted lines in Fig. 6j).

The dSSTslab(LH) is similar in both ENSO phases with its maximum in the eastern Pacific, but it is stronger by more than 1 K/K during La Niña (Fig. 9e, j, o). This asymmetry is mainly caused by the asymmetry of the MLD (dashed and dashed-dotted lines in Fig. 6j), as the LH feedback only exhibits little asymmetry. The LH feedback would even result in an asymmetry of the opposite sign, as it is a bit stronger during El Niño than during La Niña (dashed and dashed-dotted lines in Fig. 6g). The negative pole in the difference pattern in the far western equatorial Pacific is of similar strength as the positive pole in the east (Fig. 9o) and caused by a reversal of the sign of dSSTslab(LH), i.e., the LH feedback acts as a weak damping during El Niño but as a weak forcing during La Niña. The other two heat flux components do not contribute significantly to the asymmetry of the Qnet (not shown). In summary, Fig. 9 illustrates that both dSSTslab(SW) and dSSTslab(LH) contribute to ENSO asymmetry, the former in the central equatorial Pacific and due to an asymmetry in the SW feedback and the latter in the eastern equatorial Pacific and due to differences in the MLD. More specifically, the stronger dSSTslab(SW) damping during El Niño over the central equatorial Pacific prevents the SST anomalies from extending as far west as during La Niña. The stronger dSSTslab(LH) damping during La Niña over the eastern equatorial Pacific results in weaker SST anomalies there compared to El Niño. Thus, both damping terms jointly contribute to the east-west asymmetry of the SST pattern of ENSO.

In the CMIP5 models, the El Niño and La Niña patterns of the five SST tendencies individually do not differ much from the combined tendency patterns shown in Fig. 5. The models thus underestimate the asymmetry in comparison to observations, which is even more the case in WEAK than in STRONG (not shown). This becomes evident in the difference plots between El Niño and La Niña (Fig. 10). The STRONG sub-ensemble (second row) exhibits a similar difference pattern compared to observations (first row), but with a weaker amplitude. The WEAK sub-ensemble has an even smaller asymmetry in all five SST tendencies (third row), which is further illustrated in the difference plots (fourth row). The reduction in the asymmetry of the dSSTslab(Qnet) from STRONG to WEAK (Fig. 10h, m, r) in the central and western equatorial Pacific is caused by a reduced asymmetry of dSSTslab(SW) (Fig. 10i, n, s) and in the eastern Pacific by a reduced asymmetry in dSSTslab(LH) (Fig. 10j, o, t). Further, the reduced asymmetry in dSSTslab(SW) from STRONG to WEAK (Fig. 10i, n, s) can be explained by the smaller longitudinal shift of the SW feedback in response to SST change (dashed and dashed-dotted lines in Fig. 6e, f).

Fig. 10
figure 10

The difference El Niño − La Niña, in ae for observations/reanalysis (same as Fig. 9k–o), in fj for the CMIP5 STRONG sub-ensemble, in ko for the CMIP5 WEAK sub-ensemble and in pt the difference between STRONG and WEAK sub-ensembles (2. row minus 3. row); The boxes in d, i, n and e, j, o mark the region of highest asymmetry in observations/reanalysis, as used as a measure for the asymmetry of dSSTslab(SW) and dSSTslab(LH) in Fig. 11; The hatching indicates significant differences on a 90% confidence level calculated by a Student’s t test

The weaker asymmetry in the MLD over the eastern equatorial Pacific in WEAK (dashed and dashed-dotted lines in Fig. 6k, l) can explain the reduced asymmetry in dSSTslab(LH) (Fig. 10j, o, t). Finally, there is a weaker asymmetry in dSST (Fig. 10f, k, p) and dSSToc (Fig. 10g, l, q) in WEAK relative to STRONG. Thus, in models with smaller atmospheric feedback strength ENSO asymmetry as well as the asymmetries in the heat fluxes are more strongly underestimated than in models with stronger atmospheric feedbacks.

To analyze this finding in more detail, we define an index of ENSO asymmetry as the difference in skewness between Niño3 and Niño4. This index reflects that El Niño events are stronger than La Niña events and exhibit their maximum further to the east, resulting in observations in a positive skewness in Niño3 and a negative skewness in Niño4. In support of the above results, there is a significant relationship between the ENSO asymmetry and the overall strength of the atmospheric feedbacks (defined as the strength of the combined normalized wind and heat flux feedback) (Fig. 11a). This can be explained by the more realistic ENSO dynamics in the models in STRONG (Bayr et al. 2019). Further, we find indications, that the asymmetry in both dSSTslab(SW) and dSSTslab(LH) is important for ENSO asymmetry, as together they explain the tripole structure in the asymmetry in dSSTslab(Qnet) (Fig. 9m). Due to a shallower MLD during La Niña compared to El Niño, the dSSTslab(LH) is stronger during La Niña and contributes to the positive skewness of SST anomalies in the Niño3 region (Fig. 11c, d). The nonlinearity of dSSTslab(LH) and the asymmetry of the MLD, here defined as the difference between El Niño and La Niña in the Niño3 region, are simulated more realistically in models of the STRONG category than in models of the WEAK category (Fig. 11c, d). We also find a significant correlation between the asymmetry of dSSTslab(SW) and the ENSO asymmetry (Fig. 11b), as previously postulated by Lloyd et al. (2012) and Bellenger et al. (2014). The asymmetry of dSSTslab(SW), defined here as the difference between a western box (130°E–160°E, 5°S–5°N) and an eastern box (120°W–170°W, 5°S–5°N) (black boxes in Fig. 10n), dominates the asymmetry of dSSTslab(Qnet) over the western and central equatorial Pacific (Fig. 9m). This can be explained by the longitudinal shift of the rising branch of the Pacific Walker Circulation during ENSO events (Bayr et al. 2014). The above analyses strongly suggest that both the SW and the LH feedbacks in combination with the MLD asymmetry are important contributors to ENSO asymmetry. This does not rule out the importance of other factors such as the nonlinearity of the wind-SST feedback (e.g. Frauen and Dommenget 2010; Karamperidou et al. 2017), and several other oceanic and remote processes that were mentioned above.

Fig. 11
figure 11

For observation/reanalysis and CMIP5 models in a ENSO asymmetry measured by the skewness of SST in Niño3 minus skewness of SST in Niño4 on the x axis vs. atmospheric feedback strength (average strength of the normalized wind and heat flux feedback) on the y axis; b same as a but here on the y axis the asymmetry of dSSTslab(SW) measured by the difference [(130°E–160°E, 5°S–5°N) minus (120°W–170°W, 5°S–5°N)] between El Niño and La Niña in dSSTslab(SW) as indicated by the boxes in Fig. 10d, i, n; c the nonlinearity of dSSTslab(LH) measured by the difference between El Niño and La Niña in dSSTslab(LH) in Niño3 as indicated by the box in Fig. 10e, j, o on the y axis vs. the skewness of SST in Niño3 region on the x axis; d same as c but here the difference in mixed layer depth between El Niño and La Niña in the Niño3 region on the y axis; The color of the numbers indicates the members of the STRONG (red) and WEAK (blue) sub-ensembles; The yellow, orange and purple circles are the values of ERA Interim, ERA40 and OA Flux, respectively; the black line is the regression line between the CMIP5 models and the correlation is given in upper left corner, with 3 stars indicating a significant correlation on a 99% confidence level

6 Summary and discussion

In this study, we have investigated the interplay of thermodynamics and ocean dynamics during ENSO growth phase in observations/reanalysis products and CMIP5 models. We find that while in observations, 1 K SST change in the Niño3.4 region is composed of an ocean dynamical SST forcing of + 2.6 K and a thermodynamic damping of − 1.6 K, most CMIP5 models underestimate both the SST forcing by ocean dynamics and the thermodynamic damping. This is due to an error compensation between the underestimated wind and heat flux feedbacks, which affects the SST tendency most strongly in the Niño3.4 region. In this region, the SW damping, which has a major influence on the thermodynamic damping, is most strongly underestimated. Moreover, climate models with strongly underestimated wind-SST and heat flux-SST feedbacks do not exhibit “purely” ocean dynamics-driven ENSO, as observed, but instead a hybrid of ocean dynamics-driven and shortwave-driven ENSO, as described in Dommenget (2010), Dommenget et al. (2014) and Bayr et al. (2019).

Our results are based on computing the offline slab ocean SST during the ENSO growth phase. This enables to estimate the contribution of the ocean dynamics to SST change during ENSO growth from the residual of the total SST change and the SST change by the heat fluxes. In particular, we could quantify the contribution of the SW feedback to SST growth, which is erroneously positive in half of the CMIP5 models. In the most biased models, the SW feedback contributes as much to the SST growth in the Niño3.4 region as the ocean dynamics. The situation is different in the eastern part of the Niño3 region, where the SST is strongly damped by a combination of a strong LH feedback and a shallow MLD. Although the LH feedback is weaker in models with weak atmospheric feedbacks, it is still much more robust across the models than the SW feedback (Lloyd et al. 2011; Ferrett et al. 2017). Therefore, the eastern part of the Niño3 region is less affected by the error compensation.

The method applied in this study, which is based on the offline slab ocean SST tendencies due to heat fluxes, also helps to understand the failure of many climate models to simulate realistically the ENSO-related SST anomaly patterns. In these models, the patterns extend too far to the west in comparison to observations, which can be explained by the too weak heat flux damping over the western equatorial Pacific. Our results further show that due to the biased SW feedback in the models of the WEAK category, the SST anomalies are weakly damped over the central equatorial Pacific. Therefore, only a weak forcing is sufficient to generate sizeable SST anomalies.

Further, we show that the strength, seasonality and asymmetry of the atmospheric feedbacks are essential for realistically simulating the phase locking of ENSO to the seasonal cycle and the asymmetries between El Niño and La Niña. More specifically, models with stronger atmospheric feedbacks also more realistically simulate seasonality and asymmetry of the heat flux feedbacks than models with weaker atmospheric feedbacks, as well as the seasonal phase locking and the asymmetry of ENSO (Figs. 8, 11). This can be traced back to the equatorial cold SST bias in models with weaker atmospheric feedbacks, which hampers proper coupling of the Walker Circulation with upper-ocean heat content due to a biased atmospheric mean state (Bayr et al. 2019). Finally, from the heat flux feedbacks it is the SW and LH feedbacks that contribute the most to the seasonal phase locking and the asymmetry of ENSO. This again emphasizes the importance of the atmosphere for ENSO asymmetry postulated by Frauen and Dommenget (2010), as both the wind forcing and the heat flux damping strongly contribute to the asymmetry. We note that other factors too contribute to ENSO asymmetry, as discussed in Sect. 1.

As previously shown, the erroneously positive SW feedback observed in some models is caused by biases in atmospheric physics (Li et al. 2015; Ferrett et al. 2018) and can be amplified by a cold equatorial SST bias. The latter leads to an overestimation of low-level stratiform clouds that dissolve if SST is rising during El Niño (Lloyd et al. 2009, 2011, 2012; Dommenget et al. 2014; Bayr et al. 2018). The cold equatorial SST bias is a common problem in current climate models and its causes are still under discussion. Possible contributors are too strong equatorial trade winds, too sensitive ocean vertical mixing (Guilyardi et al. 2009) and biases in the cloud albedo in the subtropics that influence the equatorial SST via the subtropical cells (Vannière et al. 2013, 2014; Burls et al. 2017).

The main purpose for applying the offline slab ocean SST method was the quantification of the error compensation of the atmospheric feedbacks by estimating the contributions of the heat fluxes and ocean dynamics to the SST tendency (dSSTslab and dSSToc, respectively) during the ENSO growth phase. Using a temporally and spatially varying MLD, instead of a constant MLD of 50 m as in Bayr et al. (2019), provides more accurate estimates of dSSTslab and dSSToc especially over the far eastern equatorial Pacific, where the MLD is shallow and differs considerably between El Niño and La Niña phases. Nevertheless, the main results of Sects. 3 and 4, in which we focus on the Niño3.4 region, are very similar when a constant MLD of 50 m is used (not shown). However, with regard to the ENSO asymmetry (Sect. 5) a varying MLD is important to consider, because the asymmetry of dSSTslab(LH) strongly depends on the MLD asymmetry between El Niño and La Niña while the LH feedback itself exhibits a weak asymmetry of opposite sign. The calculation of the MLD itself is subjected to uncertainties, as there are several ways to estimate the MLD. Further, the mean state and asymmetry of the MLD varies considerably among the climate models (Wang and McPhaden 1999; Vialard et al. 2001).

More advanced measures such as the BJ index (e.g. Jin et al. 2006; Wengel et al. 2018) provide detailed information about all positive and negative feedbacks that drive equatorial Pacific SST anomalies. However, such measures require 3D ocean data and are quite intense to compute. The offline slab ocean SST method applied here has the advantage that it only needs 2D data, at least when assuming constant MLD, and is easier to implement in advanced software packages dealing with large data sets. This method is part of the “CLIVAR ENSO 2020 metric package” to measure the error compensation of ENSO atmospheric feedbacks in CMIP models (Planton et al. 2020). We note again that using a constant MLD is a major simplification: the observed MLD in the eastern equatorial Pacific is much shallower than that in the central and western equatorial Pacific. Moreover, the MLD is deeper during El Niño than during La Niña in the eastern equatorial Pacific. Therefore, applying the varying MLD is more accurate, especially over the eastern equatorial Pacific.

To conclude, we have presented a method, the offline slab ocean SST method, which enables quantifying the influence of error compensation in ENSO atmospheric feedbacks that is observed in many climate models. We show that the error compensation has its largest impact in the Niño3.4 region. Finally, stronger ENSO atmospheric feedbacks lead to a more realistic simulation of the seasonal phase locking and the asymmetry of ENSO primarily due to a better simulation of the seasonal cycle and asymmetry of the SW and LH feedbacks and the MLD.