The ENSO amplitude asymmetry is typically measured by skewness, which is a measure of the extent to which a distribution differs from a normal distribution:

where m_k is the kth moment,

where x_i is the ith observation, the overbar denotes the mean, and N the number of observations.

Here, skewness is calculated from monthly data.

The ocean heat budget can be expressed as:

Where T, u, v, and w are temperature, zonal, meridional, and vertical velocities, respectively. The overbar and prime denote climatological mean and anomalies, respectively. The first bracket represents the linear terms, the second bracket is the nonlinear terms and the residual term, R, denotes thermodynamic and subgrid-scale contributions. The heat budget terms are estimated from monthly mean anomalies with respect to the long-term monthly climatology.

ENSO asymmetry is assessed in a number of long climate simulations of the Met Office Hadley Centre general circulation model HadGEM-GC3.1 (Williams et al., 2018) in three different resolutions: low (LL), medium (MM), and high (HH). LL and MM are CMIP6 preindustrial control simulations, at N96 (atmosphere midlatitude spacing of 135 km)/ORCA1 (1° ocean with enhanced meridional resolution near the equator) and N216 (60 km)/ORCA025 (1/4°) resolution, respectively (Menary et al., 2018). Each run length is 500 years. HH has a resolution of N512 (25 km)/ORCA12 (1/12°) and forcing representative of the year 1950 (Roberts et al., 2019). The run length is 150 years. All configurations have 85 levels in the atmosphere and 75 levels in the ocean. To investigate possible reasons for the weak asymmetry in these models, we use a PPE based on a similar model, HadGEM-GC3.05, at MM (N216/ORCA025) resolution. The objective of a PPE is to represent model uncertainty by perturbing model parameters within realistic ranges. There are 19 members, the “standard variant” (STD) unperturbed member (member 0) and 18 members with combinations of perturbations applied to atmosphere and land surface parameters. No perturbations are applied to the ocean. We use 200-year simulations with CMIP5 greenhouse gas and radiative forcing for the year 1900. Heat and fresh-water flux adjustments are derived for STD and applied to all PPE members. In addition, for the 18 perturbed members, a globally uniform, member specific, heat flux adjustment is applied (Yamazaki et al., 2021). The methodology behind the choice of the perturbed parameters and the values they take is described in detail by Sexton et al. (2021), and the identifier which links the ensemble members to the parameter perturbations is shown in brackets: 0 (00000), 1 (00090), 2 (00605), 3 (00834), 4 (01113), 5 (01554), 6 (01649), 7 (01843), 8 (01935), 9 (02089), 10 (02123), 11 (02242), 12 (02305), 13 (02335), 14 (02491), 15 (02753), 16 (02832), 17 (02868), 18 (02914).

WWBs are identified here as events where the daily 10 m zonal wind anomaly, averaged from 2.5°N-2.5°S, exceeds 3 ms⁻¹ over a longitude range of at least 10° and lasts for at least 3 consecutive days in the equatorial Pacific (130°E-80°W). We define a WWB index as the sum of the maximum anomaly for each event over a given time period to give a measure of accumulated strength. In recent literature, a number of different criteria have been used for defining WWBs, leading to varying results in terms of numbers of WWBs identified (Santoso et al., 2017). The criteria used here gives an observed average WWB count of just over 12 per year which is within the range identified by other studies (e.g., Puy et al., 2016).

Model data are compared with monthly mean SST from HadISST1.1 for 1873 to 2020 (Rayner et al., 2003), monthly zonal wind stress from ERA5 for 1950 to 2020 (Hersbach et al., 2020), and monthly mean depth of the 20° isotherm estimated from the ocean temperature dataset, EN4, for 1950 to 2020 (Good et al., 2013). WWB activity is calculated using daily 10 m winds from the ERA-Interim reanalysis for 1980 to 2017 (Dee et al., 2011). Monthly mean ocean potential temperature and currents are from two products, the ocean reanalysis used in the GloSea5 seasonal forecast system (MacLachlan et al., 2015) and a recent ocean reanalysis based on the GO5 ocean model (Megann et al., 2014), both for the period 1993 to 2017. In the scatter diagrams, the time-period used for calculating HadISST skewness is chosen to be the same as that used for the corresponding reanalysis variable.

Assessment of ENSO characteristics for HadGEM-GC3.1 at MM resolution shows that this model performs favorably with observations across a range of performance metrics (Williams et al., 2018) and in historical simulations, the PPE shows a range of plausible ENSO variability (Yamazaki et al., 2021). In this study, standard deviation for Nino3.4 (170-120°W, 5°N-5°S) across the PPE ranges from 0.51°C to 0.97°C. El Niño and La Niña composite SST anomaly patterns in HadGEM-GC3.1 are also reasonably realistic and similar at MM and HH resolutions, although for LL the ENSO pattern extends somewhat too far west (Menary et al., 2018; Roberts et al., 2019). However, compared with the observations, these El Niño and La Niña patterns are too symmetrical, reflecting a lack of extreme El Niño events in the models.

Amplitude asymmetry, measured here by SST skewness is weaker than in observations in all our model simulations and versions, irrespective of resolution, as shown in Figure 1. Because skewness in the more commonly used Niño3 region (150-90°W, 5°N-5°S) is weak in the PPE, we focus on a region in the far east Pacific (FEP) (110-80°W, 5°N-5°S) that maximizes skewness in the observations (Su et al., 2010; Santoso et al., 2017).

The spatial pattern of SST skewness for the most positively skewed (SK+), standard variant (STD), and least positively skewed (SK–) members of the 19-member PPE (Figures 1A–C) can be compared with the observations (Figure 1D). STD shows a similar pattern of skewness but is much weaker than observed, with an FEP skewness of 0.34 compared with 1.02 for observations. SK- shows weakly negative skewness along the equator with an FEP skewness of −0.11. For SK+, our most skewed simulation, the westward extent of the region of positive skewness is similar to that observed but again the SK+ member skewness for the FEP region is much weaker than observed, 0.65 compared with 1.02, and the pattern is more equatorially confined.

Skewness for the FEP region for all PPE members shows diversity across the ensemble in the mean, and also a variety of multidecadal variability (Figure 1E). However, no member has a mean skewness as high as observed, or indeed a maximum skewness in any 50-year period as high as the observed mean. We also show model data from different resolution HadGEM-GC3.1 models with lower (LL), the same (MM), and higher (HH) ocean and atmosphere resolution as the PPE. These all show similar characteristics to STD and lie within the PPE range, suggesting that lack of ENSO asymmetry is insensitive to resolution in HadGEM-GC3.1

Observations are most strongly positively skewed in the far east Pacific, decreasing to become weakly negatively skewed in the region near the dateline and in a horseshoe-like pattern in the west Pacific and extending somewhat into the South Pacific Convergence zone and the Intertropical Convergence zone. We note that PPE member SK+ has unrealistically high negative skewness in these regions and that in the PPE there is a tendency for the more positively skewed members to also have unrealistically high negative skewness in the west Pacific.

In the tropical Pacific there is a close relationship between mean equatorial SST, zonal wind stress and thermocline depth in the mean state (Philander, 1990). As expected, the PPE reflects this, with significant correlations between these variables across individual member simulations.

The relationship between these key variables and skewness also shows significant correlations, with higher SST (r = 0.54, p < 0.02), weaker zonal wind stress (r = 0.57, p < 0.02) and shallower 20°C isotherm in the west pacific (r = −0.44, p < 0.1) all being related to higher skewness in the model (Figure 2). The PPE members with deeper thermocline in the western Pacific tend to have a somewhat shallower thermocline in the eastern Pacific (Figure 2E). Biases in the coupled model mean state can therefore clearly impact the model skewness and in principle could be responsible for the lack of ENSO asymmetry. However, as the observations lie well within the range of the modeled mean state for each of these assessed variables, the model mean state is not systematically biased in this respect. A realistic representation of these variables by the model does not in itself therefore necessarily lead to realistic asymmetry and we reject the hypothesis that these mean state errors are responsible for the lack of ENSO asymmetry.

The relationship between WWBs and skewness is examined using a simple index as a measure of accumulated wind burst activity (WWB index, see section Methods) The index is calculated from daily 10 m zonal winds that have previously been filtered to remove timescales longer than a season (90 days).

We find there are large differences in the annual mean WWB index across our model versions, ranging from an average of 18.3 to 108.8 ms⁻¹ per year (Figure 3A). However, there is only a weak (not significant) positive correlation (r = 0.28) between FEP skewness and WWB activity. We also note that although the two members with largest skewness have relatively high levels of WWB activity, we see that this is not a prerequisite for the PPE skewness. Indeed, the member with the lowest WWB index has similar moderate skewness to the member with the highest WWB index. Finally, we note that the WWB index for the reanalysis lies well within the range of the PPE simulations, suggesting that the model versions are not systematically biased with respect to annual mean WWB activity.

To test the importance of the modulation of wind burst activity by SST in the PPE we look at the difference in WWB activity between El Niño and La Niña events in the January to November period prior to ENSO events (Figure 3B). El Niño/La Niña events are defined here as DJF Niño3.4 anomalies exceeding a threshold of +/- 0.8°C. As has been previously demonstrated (e.g., Puy et al., 2016) there is a significant modulation of WWBs by ENSO, with more WWBs during El Niño. However, the majority of PPE members show only a weak difference, significantly underestimating the observed modulation. We also look at this relationship for two separate periods, during boreal spring, early in the ENSO growth cycle, and in boreal summer-autumn, when SST anomalies are established and growing in magnitude. We find a similar relationship for both periods (not shown), suggesting that the rather weak difference in WWB activity in the PPE may result in very large initial El Niño growth in the PPE being less likely than observed, and also, in the second part of the year, that the model may be underestimating the observed dependence of WWBs on SST. However, because of the lack of extreme El Niño events in the PPE this latter result may also simply reflect this lack of model skewness. While there is no relationship (r = 0.2) across the PPE with skewness, the difference between the model and observations suggests possible model error in HadGEM-GC3.1 in this respect.

Equatorial non-linear advection for the combined ocean reanalysis product shows warming in the region of the thermocline and in the surface mixed layer, largely concentrated to the east of the dateline (Figure 4A). These features agree well with the analysis shown by Hayashi et al. (2020). Non-linear advection for the PPE, combined over all models, shows a similar pattern, but with muted warming in the ocean subsurface (Figure 4B).

Following Hayashi et al. (2020), area averages are calculated for the region (180-100°W, 1°S-1°N, 50–150 m), encompassing the maximum subsurface non-linear dynamic warming (NDW) in the reanalysis. There is a significant relationship (r = 0.74, p < 0.001) between the FEP skewness and the NDW for the model simulations and ocean reanalyses (Figure 4C). The correlation over PPE members alone is also significant (r = 0.41, p < 0.1). Unlike our other hypothesized mechanisms for the underestimated ENSO asymmetry, there is a clear systematic model bias in the strength of the NDW and all of the PPE members underestimate the NDW relative to the ocean reanalysis. The NDW is therefore a strong candidate for the lack of asymmetry in the model, consistent with other studies on NDW (Hayashi et al., 2020; Chen et al., 2021).

To consider the contribution of NDW to the growth of ENSO, the NDW for the mean of April to November prior to the peak of ENSO is calculated for all El Niño and La Niña events. Overall, there is a positive contribution during both El Niño and La Niña for both the combined ocean reanalysis and the model members (Figures 4D,E) This acts to weaken cold (La Niña) events and strengthen warm (El Niño) events, contributing to amplitude asymmetry. However, the total warming in the reanalysis is over twice the magnitude of that in the model (see the right-hand bars in Figures 4D,E).

We now consider the three component terms of the NDW (Equation 2) to try to better understand the origin of the weak NDW in the model. The zonal non-linear warming term (-u'∂T'/∂x) and meridional term (-v'∂T'/∂y) both make a positive contribution to the warming, whereas the vertical (-w'∂T'/∂z) term is negative and acts to partially offset the horizontal terms. Comparing the PPE with the ocean analysis we see that it is the zonal component of non-linear advection which largely accounts for the difference in total magnitude (left hand bars in Figures 4D,E). It is noteworthy that there is also a significant relationship (r = 0.48, p < 0.05) between skewness and NDW for the zonal component across PPE members (not shown).

To further explore the reason for the difference in zonal NDW, we decomposed the zonal NDW for the observational reanalysis and the model (Figures 5A,B), into the zonal anomalous temperature gradient (Figures 5C,D) and anomalous zonal current (Figures 5E,F), although this separation is an approximation due to the correlation between the strength of T' and u' fluctuations. It can be seen that the anomalous eastward surface current anomalies extend too deep in the PPE compared with the reanalysis and in particular, that the anomalous westward currents in the region of the equatorial undercurrent are much weaker in the model (Figures 5E,F). This latter difference suggests that the ocean circulation response to ENSO is too weak in our model and provides a focus for future model development.