The influence of ENSO events on the stratospheric QBO in a multi-model ensemble

The Quasi-Biennial Oscillation (QBO) and the El Niño–Southern Oscillation (ENSO) are two dominant modes of climate variability at the Equator. There exist observational evidences of mutual interactions between these two phenomena, but this possibility has not been widely studied using climate model simulations. In this work we assess how current models represent the ENSO/QBO relationship, in terms of the response of the amplitude and descent rate of stratospheric wind regimes, by analyzing atmosphere-only and ocean–atmosphere coupled simulations from a large multi-model ensemble. The annual cycle of the QBO descent rate is well represented in both coupled and uncoupled models. Previous results regarding the phase alignment of the QBO after the 1997/98 strong warm ENSO event are confirmed in a larger ensemble of uncoupled experiments. However, in general we find that a relatively high horizontal resolution is necessary to reproduce the observed modulation of the QBO descent rate under strong ENSO events, while the amplitude response is generally weak at any horizontal resolution. We argue that biases in the mean state and over-dependence on parameterized wave forcing undermine the realism of the simulated coupling between the ocean and the stratosphere in the tropics in current climate models. The modulation of the QBO by the ENSO in a high emission scenario consistently differs from that in the historical period, suggesting that this relationship is sensitive to changes in the large-scale circulation.


Introduction
The El Niño-Southern Oscillation (ENSO) and the Quasi-Biennial Oscillation (QBO) are two well-known patterns in the equatorial regions, and both are major sources of interannual variability, affecting the weather and climate on a global scale.
The ENSO (Philander 1990) is a coupled atmospheric-oceanic phenomenon, in which the sea surface temperature (SST) of the tropical Pacific ocean oscillates between warm (El Niño) and cold (La Niña) states. The largest temperature anomalies are encountered in the eastern Pacific [in a region labelled '3.4', (Trenberth 1997)] but the ENSO has impacts on the global climate system. The cold and warm states are not symmetric, as the warm state features stronger but less persistent anomalies than its cold counterpart (Capotondi et al. 2015). The ENSO displays a marked periodicity between three and seven years, and is modulated by the seasonal cycle, as most events develop in boreal summer and peak in the next winter (Ropelewski and Halpert 1987;Stein et al. 2014). In the following we use the term ENSO for the index of the 3.4 region.
The Quasi-Biennial Oscillation (QBO) is a major pattern of variability of the zonal wind in the tropical stratosphere (Baldwin et al. 2001), and can be described as a quasi-periodic alternation of westerly and easterly wind shear regimes, descending from the upper stratosphere and dissipating around the tropopause level, with a periodicity close to 28 months. Observations of this pattern go back to around 1950, when regular radiosonde soundings began at near-equatorial locations. The basic mechanism of the QBO Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0038 2-020-05131 -7) contains supplementary material, which is available to authorized users. 1 3 was proposed by Lindzen and Holton (1968), and further refined by Holton and Lindzen (1972) and Plumb and Bell (1982), and it involves different types of tropical waves, such as Rossby-gravity, Kelvin and small-scale gravity waves. The generation of these atmospheric waves, which propagate upwards and deposit momentum by interacting with the mean flow, is partly related to the latent heat released by deep convective activity in the tropics (for example, Bergman and Salby 1994). At equatorial latitudes the excitation of gravity waves is dominated by such transient events, while the excitation by orographic barriers plays a lesser role (e.g., Giorgetta et al. 2002).
Previous studies have speculated that the QBO and ENSO can interact within the climate system. The fact that the QBO wind shears can affect tropical convection has been put forward by Gray et al. (1992a, b), Collimore et al. (2003), Liess and Geller (2012), Son et al. (2017), among others. These authors proposed that the low-level shear associated with the descent of the QBO is important in determining the penetration of convective systems. According to Gray et al. (1992a), this could explain the tendency of warm ENSO events to occur during the easterly phase of the QBO, provided that the heat content of Pacific warm pool is large enough. According to this mechanism, the stratospheric state would act to modify the upper tropospheric circulation in the tropics, affecting in turn oceanic conditions.
On the other hand, ENSO events affect the tropical stratosphere, primarily due to structural changes in the shallow branch of the Brewer-Dobson circulation and deep convective patterns (for a review, Domeisen et al. 2019). In particular, the thermal signatures due to ENSO are out-of-phase between the tropical troposphere and stratosphere (Randel et al. 2009). However, these effects are partly masked in observations by other factors, like contemporary volcanic eruptions.
Focusing on the QBO, Maruyama and Tsuneoka (1988) noted the short (6 months) easterly phase of the QBO at 50 hPa coinciding with a moderate warm event in 1986/87. And while the QBO has been remarkably regular for decades, during late 2015 it entered an anomalous state, as the westerly phase was interrupted by a shallow easterly jet, linked to extratropical wave forcing from the boreal hemisphere (Osprey et al. 2016;Newman et al. 2016). This circumstance occurred during one of the largest warm ENSO events ever recorded (Santoso et al. 2017), rekindling speculation on the relationship between these two modes of variability.
By analyzing the observational record, Taguchi et al. (2010) showed that the phase speed and amplitude of the QBO are different during the warm and cold phase of the ENSO. In particular, the phase speed is larger (smaller) during warm (cold) ENSO states; conversely, the amplitude of the QBO is larger (smaller) during the cold (warm) phase. These results, in particular for the QBO amplitude, are however dependent on the period under analysis, as discussed by Geller et al. (2016b). The analysis of Taguchi et al. (2010) has been further extended by Yuan et al. (2014), who have found a consistent modulation of the QBO properties and the cold-point tropopause temperature by the phase of the ENSO in a study of multiple stations in the tropics. The mechanism linking the ENSO to the QBO is not well understood, and could involve changes in convection and related equatorial wave activity due to the ENSO (Taguchi et al. 2010;Geller et al. 2016a).
A notable complication in analyzing the ENSO/QBO relationship is that external forcings play a significant role in both the stratospheric and oceanic variability. One major interference is represented by tropical volcanic eruptions, which affects the tropospheric climate (Christiansen 2008;Driscoll et al. 2012), including the occurrence of ENSO events (Khodri et al. 2017). Eruptions have noticeable impacts on the stratospheric thermal structure (Angell 1997) and can alter the morphology of the QBO, as found in numerical experiments (Aquila et al. 2014).
The current generation of climate models reproduce the observed modes of variability in the tropical ocean (Bellenger et al. 2014) and stratosphere (Butchart et al. 2018), in an increasingly realistic manner compared with earlier times. However, the skill in simulating these patterns may be due to compensating errors (Guilyardi et al. 2009), as models present notable biases in their atmospheric and oceanic components, and these phenomena depend on subgrid processes, like non-orographic gravity wave (NOGW) drag in the case of the QBO (e.g., Alexander et al. 2010).
Recent studies investigated the connection between ENSO events and the QBO from a numerical modelling perspective. Schirber (2015) studied the ENSO impacts on the QBO using a model with an interactive NOGW scheme, but the observed connection was simulated only during the westerly phase of the QBO, and the response to altered upwelling during ENSO events was too small. Geller et al. (2016b) hypothesized that the zonally averaged gravity wave momentum fluxes determine the amplitude and speed of descent of the QBO in climate simulations, and that they are modulated by the ENSO state. While more robust changes are observed for the phase speed, the changes of amplitude are more equivocal. ENSO-related changes of NOGW momentum fluxes, with larger zonally averaged forcings during warm events, are described by Alexander et al. (2017), and by Kang et al. (2018), who also obtained significant changes due to the ENSO in precipitation and stratospheric momentum fluxes within simplified models. Considering instead the results of an atmospheric general circulation model (GCM) integration, Christiansen et al. (2016) analyzed the effects of warm ENSO events on the stratospheric QBO. They noted that the observed tropical zonal wind in the stratosphere is aligned between two and three years following strong warm 1 3 ENSO events. They also found an alignment of the phase of the QBO after the strong warm ENSO event in 1997/98 in ten atmospheric model integrations, as manifested by a reduction of the ensemble spread.
The previous studies cited above documented consistent changes of the wave forcing and the QBO properties due to ENSO events, and have been based on simplified models or single AGCMs. In this work we study the modelled ENSO/QBO interaction in a large ensemble of multiple climate models, considering both atmosphere-only and ocean-atmosphere coupled simulations. The general properties of the QBO, their seasonality and response to ENSO events are therefore studied in different model configurations.
In particular, the influence of the horizontal resolution is analyzed in atmosphere-only simulations, with observed SSTs. The spatial resolution is important for the representation of atmospheric waves (Horinouchi et al. 2003), which are modulated by extreme ENSO states (Li and Lau 2013). It must be noted that the vertical resolution in the stratosphere is also crucial for determining the simulation of stratospheric processes (e.g., Fujiwara et al. 2007). Furthermore the large number of simulations available, performed with the same climate model, allows us to test the results of Christiansen et al. (2016) on the phase locking of the QBO due to warm ENSO events.
In coupled simulations, as opposed to atmosphere-only experiments, the ocean is freely evolving with the atmosphere, allowing for internally consistent interactions. This allows us to explore the changes in a future climate scenario, characterized by changes of the ENSO pattern (Collins et al. 2010). In fact, the frequency and amplitude of extreme ENSO events are projected to increase under higher greenhouse gases (GHG) emissions (Cai et al. 2014). This could also change the likelihood of events similar to the recent disruption of the QBO, in which the 2015/16 warm event played a role (Barton and McCormack 2017;Hirota et al. 2018). It is important to note that forced and natural changes in the atmospheric state have the potential of changing the impacts and teleconnections of ENSO, as already observed during the last decades (Yeh et al. 2018).

Definition of QBO and ENSO indices
In order to concisely describe the variability of the ENSO and the QBO, three time series are calculated for each dataset and used in our analyses.
For the tropical Pacific Ocean, the warming and cooling of the sea surface temperatures is well represented by the Niño 3.4 index, calculated as the average of the SSTs in the region 5 • N-5 • S, 120 • E-170 • W (Trenberth 1997). Before the analysis, the linear trend of the SST field is removed, as done in other studies (Calvo et al. 2017). Anomalies of the time series are then computed by subtracting the long term climatology, and the anomalies are standardized to unit variance. Note that the procedure above is done separately for each sub-period we analyze, as explained in more detail in the next section.
Following Taguchi et al. (2010), the time series are finally smoothed with a 5-months running mean to remove high frequency signals. Warm and cold ENSO events are then identified as the periods exceeding the 90th and 10th percentiles, for at least 6 consecutive months, for each dataset (Wolter and Timlin 2011). The ENSO time series for the reference dataset (ERA-Interim) is reported in Fig. 1, and the results for two arbitrary models can be found in the Supplement ( Fig. S1 and S2).
The QBO is characterized by alternating westerly and easterly zonal wind shears, so that easterlies are found below westerlies aloft, and vice versa. An empirical orthogonal function (EOF) decomposition (e.g., Wallace et al. 1993) allows to reduce the dimensionality of the phenomenon into a small number of principal components (PCs), each associated with a vertical profile of zonal wind (the corresponding EOFs). Before the analysis the zonal mean zonal wind, averaged over the equatorial region (+/− 5 degree latitude), is interpolated onto the same pressure levels (70,50,40,30,20,15,10 hPa), and the seasonal cycle is removed, as in Taguchi et al. (2010). We have verified that the precise choice of the vertical interval, or the use of anomalies from the long-term average instead of deseasonalized quantities, do not change the main results of this study. The first two principal components combined explain typically more than 90 % of the total variance (around 60 and 30 %, respectively, as indicated in the lower left corner of Fig. 2a), and can be used to compactly represent the original bidimensional (time-altitude) variability. While the EOFs for the reanalysis and models are qualitatively similar ( Fig. S3-S6 in the Supplement), the QBO structure is generally shifted upwards in models compared to observations (Schenzinger et al. 2017). This can clearly be seen in Fig. S7 of the Supplement, where EOFs for models and reanalysis are compared. Fig. 1 Time series of the ENSO index from ERA-Interim data. Periods when the index is above the 90th or below the 10th percentiles are shaded. Warm events occur in 1982/83, 1987, 1991/92 and 1997/98. Cold events occur during 1988/89 and 1998 By combining the first two PCs, it is possible to define a polar representation of the QBO (Christiansen et al. 2016), using the two time series as the component of a counterclockwise rotating vector. The quantity A = √ PC 2 1 + PC 2 2 , represents the amplitude of the QBO, while = arctan(PC 1 ∕PC 2 ) indicates the QBO phase. The phase speed (or descent rate, ′ ) is then computed by two-sided finite differencing the unwrapped phase. A 5 months smoothing is finally applied to the resulting QBO indices, as done for the ENSO time series. The results for the QBO decomposition applied to the ERA-Interim zonal wind are shown in Fig. 2, and those for two different models are included in the Supplement.

Description of the models and simulations
In this work we study the impacts of extreme states of the tropical Pacific Ocean on the equatorial stratospheric QBO, as simulated in current comprehensive climate models. In so doing we consider two types of simulations: atmosphereonly (AMIP), in which the evolution of the SST is prescribed from observations, and coupled simulations, where the ocean can interact with the atmospheric component. Given that our focus is the change in the stratospheric QBO, we have selected models that internally generate a QBO-like oscillation of the equatorial zonal mean zonal wind. All the data used for our analyses are monthly averages, as available from the respective archives.
The reference we use is the European Centre for Medium Range Weather Forecasts (ECMWF) atmospheric reanalysis ERA-Interim (Dee et al. 2011), which assimilates observations from different sources into an atmospheric model integration (performed with the Integrated Forecasting System, IFS). The representation of the QBO in ERA-Interim strongly depends on the assimilation of radiosonde observations from equatorial sounding stations ). Here we consider the period 1979-2008 for ERA-Interim, in order to directly compare the reanalysis results with those of the AMIP simulations.
The model simulations used in this study are listed in Table 1. For the AMIP analysis, we consider the EC-EARTH 3.1 model (Hazeleger et al. 2012) in a number of different configurations, where the horizontal resolution is modified. The simulations have been performed for the Stochastic Physics HIgh resolutioN eXperiments (SPHINX) project , starting from a relatively low resolution (T159, 100 km) up to a very fine resolution (T1279, 10 km). In the coupled version of the model, the T255 (60 km) resolution is used. We label the different EC-EARTH model configurations with two letters, as seen in the first column of Table 1. For the atmosphere-only runs, the first letter after the model name indicates the resolution, as C (coarse), L (low), M (medium), H (high) and U (ultra-high), while the coupled configuration is indicated with O (ocean). The second letter indicates whether the parameterizations are deterministic (D) or stochastically perturbed (S). In total there are ten model configurations for the AMIP and two for the coupled simulations, with varying ensemble sizes for each configuration. The stochastic components, which can be activated in the model atmosphere, are the stochastically perturbed parameterized tendencies (SPPT) (Palmer et al. 2009), and the stochastic kinetic energy backscatter scheme (SKEBS) (Shutts 2005). These schemes are only active in the troposphere, and they have been shown to improve the modelled precipitation on short timescales ) and the interannual variability of ENSO . The MPI-ESM-MR (Schmidt et al. 2013) AMIP runs from the fifth phase of the Coupled Model Intercomparison Project (Taylor et al. 2012, CMIP5) project are included, in order to compare with the results of the atmosphere-ocean simulations of the same model. This model is labelled MPI-ESM-MR-A, Fig. 2 The results of the decomposition procedure applied to the quasi-biennial oscillation of ERA-Interim. The two leading empirical orthogonal functions are reported in (a), black line for the first and gray for the second; the explained variance is indicated in the lower left corner. The time series for the amplitude and the phase speed, obtained by combining the leading two principal components, are shown in (b) and (c), respectively to be distinguished from its coupled version. As reported in Table 1, all the atmosphere-only simulations cover the period 1979-2008. For the coupled experiments, we consider models from both the CMIP5 and the SPHINX projects. For the CMIP5, five models are able to internally generate a QBO-like oscillation, namely CMCC-CMS (Manzini et al. 2012), HadGEM2-CC (Martin et al. 2011), MIROC-ESM and MIROC-ESM-CHEM (Watanabe et al. 2011) and MPI-ESM-MR (Schmidt et al. 2013). For the SPHINX project, only the EC-EARTH model with a resolution of T255 and 91 vertical levels is available. Note that the horizontal resolution of the coupled SPHINX simulations is higher than that of most CMIP5 models. We analyze two periods from the CMIP5 and SPHINX datasets [the latter follows the same protocol of historical and high emission scenario of the RCP 8.5 (Riahi et al. 2011)], with one period  from the historical simulation, and the other from the future (2055-2100) to highlight the climate change signal.
As indicated in Table 1, the number of realizations for each model varies, going from just one to ten ensemble members. In any case, all the models we analyze have a large number of vertical levels, and parameterize the effects of unresolved NOGW drag in the middle atmosphere: CMCC-CMS, MIROC, MPI-ESM-MR using the Dopplerspread scheme of Hines (1997a, b), while HadGEM2-CC and EC-EARTH use variants of Warner and McIntyre (2001) spectral scheme. The tuning of the NOGW scheme, different for each model, is crucial for the modelled QBOs. As discussed by Geller et al. (2016a), the QBO period in models is dependent on the magnitude of the wave momentum flux, while its amplitude depends on the phase speeds of the NOGW spectrum.

Climatology and seasonality of the QBO
In order to assess the modelled changes of the QBO induced by ENSO, we begin by describing the climatological and seasonal properties of the QBO, compared against the results for the reference dataset, ERA-Interim.
The averages of the phase amplitude and speed of the reanalysis and AMIP QBOs are reported in Fig. 3a. The model closest to ERA-Interim (denoted by the black star) is MPI-ESM-MR-A, which has a QBO slightly faster and with larger amplitude than the reanalysis. As well known, the QBO period of ERA-Interim is close to 28 months (0.22 radian/month). While the phase speed of the EC-EARTH Table 1 The models and respective simulations analyzed in this study In the first column, the letter following the EC-EARTH models indicates their configuration (as C for coarse, L for low, M for medium, H for high, U for ultra-high and O for coupled), while the last letter marks deterministic (D) or stochastic (S) versions. The two letters in the second column specify the simulation type: the first letter denotes historical (H) or scenario (S) runs, while the second indicates if atmosphere-only (A) or coupled (C). The third column, for the atmospheric model resolution, reports either the spectral triangular truncation ( T n , n being the maximum wave number), or the grid increments in degrees, and the vertical resolution is given by the number m of levels ( L m ). When two numbers are given in the fifth column, they refer to deterministic and stochastic model versions, respectively simulations compare well with the reanalysis, it generally increases with finer spatial resolution, becoming larger than the reference. When stochastic parameterizations are active, the phase speed is typically higher, with the largest value for the ultra-high resolution. On the other hand, the amplitude of the QBOs is generally only half that of ERA-Interim, and grows weakly with increasing resolution. It should be noted that the parameterized momentum fluxes are systematically reduced with increased resolution (their Table 1 Davini et al. 2017), since smaller scales are explicitly resolved at higher resolutions. While on average the phase speed of the SPHINX simulations is close to that of ERA-Interim, their interannual variability is generally larger (not shown). We note that the uppermost level is 0.01 hPa in both EC-EARTH and MPI-ESM-MR-A, and the number of vertical levels is comparable (91 and 95, respectively, equivalent to a resolution around 600 m in the stratosphere), but their horizontal resolutions are very different. The results above indicate that the tuning of the NOGW scheme can be crucial to obtain a modelled QBO with realistic mean phase speed and amplitude.
The results for the coupled simulations are reported in Fig. 3b (note that the period for ERA-Interim is 1979-2008, as in the AMIP simulations). The historical simulations by different models (circles) are characterized by a large spread in both the amplitude and phase speed of the QBO. Except the MPI-ESM-MR model, the simulated amplitude is generally smaller than the reanalysis, and the results for EC-EARTH match the atmosphere-only version of the same model. While the phase speed is in the same range of the AMIP simulations, the amplitudes vary significantly between models. The model which best matches the reanalysis in this case is HadGEM2-CC. Note the CMCC-CMS is an outlier, underestimating both the amplitude and phase speed (and also has a very irregular QBO, not shown).
For the projections, the amplitude is consistently reduced across models, compared to the historical period (note how squares are always on the left of circles). This is in agreement with Kawatani and Hamilton (2013), who argued that the stronger upwelling, expected in a warmer climate (already in the intermediate RCP 4.5 scenario), would reduce the QBO amplitude in the lower stratosphere. Note that the EOF decomposition of the QBO we apply provides a vertically integrated view (as illustrated in Fig. S8 for one of the models). Changes in the phase speed are less consistent than those in the amplitude, but in general, with the exception of the two MIROC models, the descent rate of the QBO is faster in the future, compared to the historical period.
When studying the effects of the warm and cold ENSO states on the stratosphere, the seasonality is an important 1 3 feature of both the QBO and ENSO. Warm ENSO events generally tend to peak during the boreal winter (e.g., Stein et al. 2014) even if exceptions (such in the 2015/16 case) are noted (Santoso et al. 2017), and cold ENSO states tend to last longer. The seasonality is also an important feature of the descent rate of the QBO (the variability of the amplitude is relatively smaller, and is not discussed further). This holds even when deseasonalized equatorial zonal winds are considered. A manifestation of the seasonality is the clustering of east/west and west/east transitions of the QBO (see, e.g., Dunkerton 1990;Baldwin et al. 2001;Christiansen 2010). The mechanism governing the seasonality of the QBO properties is thought to be linked with the variability of the wave forcing (e.g., Peña-Ortiz et al. 2010), convective patterns and upwelling (Dunkerton 2017). While the seasonality of the QBO phase speed is more evident, it also depends less clearly on the phase itself (Christiansen et al. 2016). For this reason longer records than those we consider here should be used for investigating these aspects. Indeed by looking at Fig. 3c, it is clear that the phase speed is characterized by a marked annual cycle in ERA-Interim, with faster descent between March and June. This result is consistent with observations (e.g., Wallace et al. 1993), as the QBO in the reanalysis is mostly due to assimilated radiosondes. The AMIP simulations of EC-EARTH generally display a realistic annual cycle, which is smaller for coarser resolutions and increase in amplitude with the resolution, and is particularly large (from 0.2 to 0.3 radian/ month) in the higher resolution version of EC-EARTH with stochastic physics (US). The annual cycle in MPI-ESM-MR-A is quite different from that of ERA-Interim, since the slower phase speed occurs in boreal summer.
The coupled experiments, reported in Fig. 3d, show more variability in the seasonality of the QBO phase speed compared to the AMIP simulations. The seasonality in the historical period is qualitatively similar to that of the reanalysis, with the exception of CMCC-CMS. The changes emerging from the projections are not conclusive: some models indicate that the annual cycle will remain similar in the future scenario, while others (like HadGEM2-CC) indicate that more marked seasonality will appear.

Response of the QBO to ENSO events
As discussed above, the models considered here reproduce the general properties of the observed QBO fairly realistically. Now we turn our focus to the response of the QBO related to the occurrence of ENSO extremes, by compositing the QBO properties during such events.
From the observational study of Taguchi et al. (2010), faster and slower QBO phase speed are expected during warm and cold ENSO phases, respectively. This pattern is well reproduced in ERA-Interim despite a relatively large uncertainty, due to the small sample size. A realistic modulation of the QBO phase speed caused by ENSO events in the EC-EARTH model has already been described by Christiansen et al. (2016), and an overview of the results for the descent rate and the effects of ENSO is provided in Fig. 4. Note that the climatological averages both in Figs. 4 and 5 are also shown in Fig. 3.
In AMIP runs of EC-EARTH (Fig. 4a), the modulation of the phase speed due to ENSO is realistically reproduced, Fig. 4 The phase speed of the QBO averaged during warm (black) and cold (grey) ENSO events, and +/− one standard error of the mean (defined as the standard deviation divided by the square root of the unadjusted sample size). Climatological averages are indicated with crosses. Statistics are computed by pooling across all months during the selected ENSO phase and available ensemble members. Results for atmosphere-only (a), coupled historical (b) and coupled scenario (c) experiments 1 3 though the models differ in their average values. Increased separation between the composited phases is found when the resolution is finer (going from left to right), in better agreement with observations. Given the dependence of the standard error of the mean on the ensemble size, the uncertainty increases at higher spatial resolution (see Table 1 for the number of realizations). The use of stochastically perturbed parameterizations results in similar separation as that in the ensemble of deterministic simulations. Moreover there is no significant response of the QBO to extreme ENSO events in the MPI-ESM-MR-A model.
For the coupled simulations of the historical period (Fig. 4b), the EC-EARTH simulations generally reproduce the modulation of the phase speed, while the CMIP5 models show little or no difference, with the exception of CMCC-CMS, even if the error bars for the cold and warm phases overlap. No response is found for both models with a relatively fast (e.g., MIROC-ESM) and slow (e.g. MPI-ESM-MR) QBO. In the scenario simulations (Fig. 4c), despite there being some differences compared to the historical period, the CMIP5 models still show little or no modulation, and only the EC-EARTH OD results are similar to the observed response.
The other property of the QBO we analyze during cold and warm ENSO states is the amplitude. According to Taguchi et al. (2010), the amplitude is smaller (larger) during warm (cold) ENSO states, possibly due to related impacts on the precipitation patterns. In Fig. 5 we display the average conditions, and the composite amplitudes during the cold and the warm phases for the three sets of simulations. The differences between the two ENSO phases for the reanalysis are consistent with observational results, with a net change of the amplitude between the two phases, but not as significant as for the phase speed.
In the AMIP simulations of EC-EARTH (Fig. 5a), it is clear that the mean amplitude in the simulations is underestimated in general, as already seen in Fig. 3a. The differences between cold and warm ENSO are not reproduced in AMIP simulations, and the small changes are actually opposite to those observed. Only the ultra-high resolution model with stochastic physics (EC-EARTH US) simulates the observed changes, but with a single realization uncertainties are large. Also the MPI-ESM-MR-A model, which has a more realistic amplitude, does not capture the observed QBO amplitude separation between ENSO phases.
For the coupled models (Fig. 5b) the amplitude changes in the historical period are generally small and not consistent with the reanalysis, except for the CMCC-CMS model, for which uncertainty is large. Interestingly, in the projections (Fig. 5c) several models predict that the amplitude is smaller than average during both ENSO phases, suggesting a nonlinear response. Christiansen et al. (2016) showed how the phase of the QBO is partly synchronized in an EC-EARTH ensemble of AMIPtype experiments following a strong warm ENSO event. This does not occur if climatological SSTs, without the observed ENSO variability, are imposed. They characterized the alignment of the ten member ensemble using the standardized ensemble spread as a measure of coherence.

Phase locking of the QBO in EC-EARTH
Here we repeat the same analysis using the larger (sixty AMIP realizations) EC-EARTH ensemble of SPHINX, which includes models configured from lower to high horizontal resolutions, and with stochastic parameterizations in the troposphere. The coherence is calculated as the spread of the equatorial zonal mean zonal wind across all the realizations, standardized by the interannual standard deviation. The significance is estimated with a bootstrap procedure, using 1000 surrogates, constructed by shuffling blocks of twelve consecutive months, randomly chosen across all the realizations, disregarding the first ten year of each experiment so as to reduce the impact of the initial conditions. The coherence is finally smoothed with a 11 months running average. The coherence is deemed significant at the X % level when it is below the Xth percentile of the surrogate distribution.
The results for the 60-member EC-EARTH ensemble are shown in Fig. 6. We can see that, after the initial period (nearly five years) in which ensemble coherence arises from the common initial conditions (and possibly also from the 1982/83 warm event), two periods of significant coherence (below 0.9) are found around 1990 and again in 2000. Both years, lying beyond the horizon of predictability of the QBO (Scaife et al. 2014), follow warm ENSO events (1987 by around one cycle of the model QBO. This result confirms the findings of Christiansen et al. (2016), and corroborates the hypothesis that warm ENSO events are able to determine the phase of the QBO, at least for this model. Further research is needed to verify if and how similar results can be achieved in a model with a more realistic QBO amplitude.
We also repeated the calculation by clustering models according to their configuration (e.g, considering the 10 coarse resolution and deterministic simulations as a subensemble), in order to test the robustness (not shown). Individual sub-ensembles display different coherence levels after initialization (ranging from 1 to 10 years), and while the vertical coherent patches around 1990 and 2000 are generally present, their significance is more variable (not shown). This difference could be expected from the reduced and variable ensemble size.

ENSO imprint in tropical convection
The response of the QBO properties to both warm and cold ENSO states is likely due to changes in the wave-mean flow interaction in the equatorial stratosphere. It is known that marked SST variability can be conducive to wave generation, ranging from small-scale gravity to planetary waves (e.g., Deckert and Dameris 2009;Scaife et al. 2016). It is however not clear how the occurrence of ENSO events alters the distribution of tropical convection, known to be a source of atmospheric waves. The hypothesis of Geller et al. (2016b), with more widespread convection occurs during warm events (and increased gravity wave momentum fluxes), and more localized but deeper convection during cold events (with a broader phase speed spectrum), has yet to be confirmed (Kang et al. 2018).
To verify the impacts of the ENSO phases on modelled convective patterns, we show in Figs. 7 and 8 the composited differences of the outgoing longwave radiation (OLR) between warm and cold events. Composites are constructed by averaging also over different ensemble members, whenever available. The OLR is a known proxy for the presence of deep convection, and ENSO events can be effectively characterized by changes in this variable (Chiodi and Harrison 2010).  The strong dependence of convective precipitation from parameterized processes (Maher et al. 2018), manifested by the very different climatologies (not shown) and the composites, is directly relevant for simulated wave motions (Horinouchi et al. 2003;Lott et al. 2014). The pattern in ERA-Interim (Fig. 7a) has a negative anomaly (up to −50 W m −2 ) over the Pacific Ocean and a positive anomaly over the Maritime Continent (up to 30 W m −2 ). The positive anomaly has a horseshoe shape with off-equatorial extensions, which are related to subsident motions compensating for the increased convection over the Pacific Ocean (e.g., DeWeaver and Nigam 2004).
The MPI-ESM-MR-A model (Fig. 7b) reproduces the spatial pattern of the reanalysis quite well, though positive anomalies are slightly smaller but more extended in the subtropics, and larger anomalies are found over the Pacific Ocean. Moreover, the OLR anomalies do not extend to the Indian and the eastern Pacific oceans as in the reanalysis. For the EC-EARTH CD simulations (Fig. 7c), the spatial pattern of the reanalysis is again qualitatively well captured, even if the negative anomalies are two times larger. No clear changes of the OLR composites with horizontal resolution can be seen for EC-EARTH (not shown), with the negative anomaly over the Eastern Pacific Ocean generally too large.
The coupled simulations are shown in Fig. 8, with historical (HC) and scenario (SC) experiments side-by-side, so to directly compare the two periods. There are considerable differences between individual models, likely due to their ocean components, which produce ENSO with different properties, and their convective schemes, which respond differently to SSTs.
Considering first the historical period, it can be noted how models such as those MIROC-based, or MPI-ESM-MR, produce relatively small differences between ENSO phases, compared with the reanalysis, over the Pacific and Indian oceans. Even if they share the same atmospheric model (ECHAM5), the results for MPI-ESM-MR and CMCC-CMS are quite different, as the Pacific negative anomalies of the latter are deeper and over a wider area. The composites for the coupled EC-EARTH models are instead similar to their atmosphere-only correspondent.
The differences between the historical (HC) and scenario (SC) simulations are remarkable, and some consistent properties can be noted (Fig. 8). The magnitude of both positive and negative anomalies are larger in most models, possibly as a result of increased water vapor concentration (Held and Soden 2006). The anomalies of MIROC-based models remain nonetheless modest in particular for the scenario (however, also the climatological OLR is very different from that observed, not shown). For most models, the projected ENSO impacts on convective patterns appear to strengthen (with larger and more widespread anomalies), likely leading to the changes in the excitation of atmospheric waves, and increased upwelling Garny et al. 2011). Changes in the wave forcing during warm and cold ENSO phases could be related to the nonlinear change of the QBO amplitude (smaller than average in both cases) shown in Fig. 4c. A more quantitative analysis of the dynamics near the tropical tropopause would be needed to characterize these processes, and should be explored in future studies.

Summary and conclusions
In this work we assess the capabilities of current climate models to reproduce the observed QBO and its modulation due to strong ENSO events. While these climate phenomena are internally generated in state-of-the-art models, they are strongly dependent on parameterized and unconstrained processes. Given the recent and unanticipated QBO disruption in 2016 (Osprey et al. 2016), and the possible influence of the concurrent warm ENSO event (Barton and McCormack 2017; Hirota et al. 2018), it is important to determine the capability of current climate models in reproducing the observed modulation of the QBO by the ENSO, which is still not well understood (Taguchi et al. 2010). In this work we show that the observed ENSO/ QBO relationship in current climate models is generally poorly reproduced, likely as a consequence of the coarse spatial resolution and the reliance on stationary parameterizations. Some improvements for the phase speed are found when the horizontal resolution is increased, at least for a single model with a relatively weak QBO. Robust decrease of the QBO amplitudes is found for the RCP 8.5 scenario, due to changes in the global circulation caused by increased GHG concentrations (Kawatani and Hamilton 2013;Hardiman et al. 2014).
The variability of the ENSO and the QBO are represented by means of widely used indices. The ENSO index is based on the SST anomalies in the tropical Pacific Ocean, while the QBO indices represent the amplitude and the descent rate of the zonal wind shears in the equatorial stratosphere. The procedure is applied on a reanalysis dataset (ERA-Interim) and on a large ensemble of both atmosphere-only and atmosphere-ocean coupled simulations, from the CMIP5 archive and the SPHINX project. Previous numerical studies of such kind have been based on simplified models, or on a single climate model. Here we analyze multiple models, considering both the historical period and a scenario with high emission of GHGs, in which both the ENSO and the QBO are projected to change.
In current GCMs, parameterizations are tuned to achieve a realistic climate. For the QBO, most models achieve a realistic descent rate (including its seasonal variability), as a result of the NOGW drag scheme 1 3 adjustments. The QBO amplitudes are however generally smaller than that of ERA-Interim, particularly in the EC-EARTH model. Note that in most models peak QBO amplitudes are reached at higher altitudes than the reanalysis, as also inferred from the EOF vertical profiles. Consistent decreases of the QBO amplitude in the projection are found across all models, as known from previous studies. Moreover, large differences are found in the projected changes of the seasonal cycle of the QBO phase speed in coupled simulations, which could result from the interaction with the extratropics (e.g., Anstey and Shepherd 2014).
For the atmosphere-only experiments, most simulations analyzed here are performed with the EC-EARTH model. This model is known to simulate a realistic change of the QBO phase speed due to ENSO events (Christiansen et al. 2016), as also found in this study. Indeed, the modulation of the phase speed is realistically reproduced, increasingly better as the spatial resolution of the model is increased. Also the differences in convective patterns, as inferred from the OLR proxy, are realistic for this model. This is not the case for the amplitude of the QBO in the model, which is largely insensitive to the occurrence of anomalous SST in the tropical Pacific Ocean. No modulation is found for the MPI-ESM-MR-A model, which has a coarser horizontal resolution compared to EC-EARTH, neither for phase speed nor amplitude. The mean amplitude and phase speed of the QBO in this model are however very close to that of the reanalysis, likely thanks to a fine tuning of the Hines NOGW scheme (which employs a stationary source spectrum).
In coupled experiments, the occurrence of warm and cold ENSO events depends only on internal model dynamics. Among the coupled models considered here, only EC-EARTH reproduce the observed phase speed modulation and OLR anomalies for the historical period. No appreciable signal is found for the CMIP5 models we analyzed, which have lower spatial resolutions, and QBOs with realistic amplitudes. Despite large biases, the response is qualitatively similar to that of the reanalysis for both the CMCC-CMS and the MPI-ESM-MR models. An interesting feature in several coupled models is the modulation of the QBO amplitude, due to cold and warm ENSO events, in the scenario experiments. In both phases, the amplitude of the QBO is smaller than its average throughout the whole scenario period. This is concomitant with the weakening of the QBO, documented in previous studies (Kawatani and Hamilton 2013) and further supported here. We also find evidence of stronger convective anomalies during both ENSO phases in the scenario, compared to the historical period. We argue that, together with the increased tropical upwelling under the scenario Hardiman et al. 2014), this may be another indication of the dependence of the ENSO and the QBO relationship on the mean state. The projected weakening of the QBO in the lower stratosphere likely affects the ENSO feedbacks, leading to qualitative changes in the future. We note in passing that recently Garfinkel et al. (2018) described a nonlinear response of the stratosphere to strong ENSO events, for temperature and water vapor.
The underestimation of the QBO amplitude (by about a factor two) in the EC-EARTH model can also help interpreting the realistic response of the QBO phase speed to ENSO events, described by Christiansen et al. (2016). By repeating their analysis, we found that a sixty member AMIP ensemble of EC-EARTH becomes coherent about after two years following the strong warm event in 1997/98, possibly as a result of changes in the stratospheric wave forcing due to warmer SSTs. This alignment is also present in smaller subsets of the sixty member ensemble. We argue that the model's weak QBO could be strongly affected by the wave forcing related to anomalous convection.
Large uncertainties persist on the mechanism by which the ENSO can modulate the QBO, and how this relationship could change over time. For better understanding this process, it will be necessary to further investigate the changes in surface climate due to ENSO events. It is worth mentioning that external factors, like volcanic eruptions, may affect the results presented here, but it is not possible to isolate them systematically in these simulations. It will be important to analyze the components of the momentum forcing of the stratosphere (with a focus on the role of NOGW drag), considering the latest reanalyses and additional model experiments, such as those undertaken for the key intercomparison projects described by Gerber and Manzini (2016) and Butchart et al. (2018). otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.