Boreal winter stratospheric climatology in EC-EARTH: CMIP6 version

The performance of the European Consortium Earth-system model (EC-EARTH) in the boreal winter stratosphere is comprehensively assessed for the first time, in particular its version 3.3 that contributes to CMIP6. A 100-year long simulation with prescribed climatological boundary conditions and fixed radiative forcing, representative of present-day climate, is used to evaluate the simulation of the climatological stratospheric circulation and to identify model biases. Results show that EC-EARTH has a large issue with the vertical distribution of stratospheric temperature from the tropics to mid-latitudes, seemingly linked to radiative processes of ozone, leading to a biased warm middle-upper stratosphere. Associated with this model bias, EC-EARTH simulates a stronger polar vortex at upper-stratospheric levels while the Brewer-Dobson circulation at middle/lower levels is weaker than reanalysis. The amplitude of the climatological planetary waves is overall underestimated, but the magnitude of the background wave injection from the troposphere into the stratosphere is overestimated, related to a weaker polar vortex at lower-stratospheric levels and, thus, a less effective wave filtering. This bias in the westerly flow could have a contribution from parameterized waves. The overestimation of background wave driving is maximum in early-winter, and may explain the overestimated frequency of sudden stratospheric warmings at this time, as compared to reanalysis. The spatial distribution of wave injection climatology has revealed a distinctive role of the climatological planetary waves: while large-scale waves (wavenumbers 1–2) dominate the eddy heat flux over the North Pacific, small-scale waves (wavenumbers 3–4) are responsible for the doubled-lobe structure of the eddy heat flux over Eurasia. EC-EARTH properly simulates this climatological feature, although overestimates its amplitude over central Eurasia.


Introduction
The stratosphere is now recognized as a fundamental part of the atmosphere, not as a passive, subordinate layer of the more chaotic and variable troposphere. This recognition relies on its chemically and dynamically active role driving Froila M Palmeiro fm.palmeiro@meteo.ub.edu 1 the former depends on the conditions of the latter. During winter, upward planetary wave propagation is limited by the strength of the zonal wind according to its wavenumber (Charney and Drazin 1961), and when wave breaking occurs, the deposited momentum modifies the background conditions and thus wave propagation (e.g. Matsuno 1970;Christiansen 1999;Kodera et al. 2000).
Another important phenomenon modulating the vortex variability is the Quasi-Biennial Oscillation (QBO; Baldwin et al. 2001 and references therein). The QBO is described by alternating phases of descending westerlies and easterlies from the stratopause to the lower stratosphere in the tropics, and can modulate the extratropical circulation via the Holton-Tan effect (Holton and Tan 1980). Responding to the wave propagation criteria, more waves are deflected to the pole during easterly QBO phases, so the polar vortex becomes weaker and warmer, while the opposite occurs during westerly phases (e.g. see Anstey and Shepherd 2014 for review). In addition, since both are wave-driven, there is also an expected influence of tropical upwelling on the QBO (Dunkerton 1997).
A correct representation of all these stratospheric phenomena is fundamental in global climate models (GCMs) to simulate realistic dynamics. It has been shown that biases in stratospheric processes have an impact on stratosphere-troposphere coupling (Gerber et al. 2012;Charlton-Pérez et al. 2013), tropospheric climatology (Shaw et al. 2014) and variability (Gerber et al. 2010), and constitute a source of uncertainty in future climate projections (Manzini et al. 2014). The Stratosphere-troposphere Processes And their Role in Climate (SPARC) project, part of the World Climate Research Program (WCRP), includes several initiatives to detect model biases in the stratosphere and advance towards their elimination: from the more general Dynamical Variability activity (DynVar 1 ; e.g. Gerber et al. 2012) or the Network on Assessment of Predictability (SNAP 2 ; e.g. Tripathi et al. 2015), to the more processoriented QBO initiative (QBOi 3 ; e.g. Butchart et al. 2018) among others.
The present study presents the first comprehensive assessment of the stratospheric climatology in the state-ofthe-art EC-EARTH model and complements the study of Palmeiro et al. (2020) who reported stratospheric variability in its previous version 3.1. This model version (v3.3) is the one contributing to the Coupled Model Intercomparison Project phase 6 (CMIP6) and other international projects and initiatives, such as QBOi. The simulation analysed here, and described in Sect. 2, constitutes the control 1 https://www.sparc-climate.org/activities/dynamical-variability/.
3 https://www.sparc-climate.org/activities/quasi-biennial-oscillation/. downwelling and impacts associated with the stratospheric polar vortex variability. Stratosphere-troposphere coupling is of especial interest in subseasonal and seasonal climate forecasting, since the impact of the slowly varying stratosphere can be used as a predictor for tropospheric circulation beyond a week (e.g. Christiansen 2005;Sigmond et al. 2013;. The transition layer between the troposphere and the stratosphere, so-called the upper troposphere -lower stratosphere (UTLS) region, is dynamically active and modulates atmospheric composition, chemistry, and the radiation budget of the Earth (e.g. Forster and Shine 1997;Shepherd 2007). In the tropical UTLS, the troposphere and the stratosphere are well separated by the level of the lowest temperature in boreal winter, which is called the cold point tropopause (CPT). The CPT controls the amount of water vapour entering the tropical stratosphere, which is fundamental in the radiative balance (Highwood and Hoskins 1998). The main entry of chemical species to the stratosphere occurs through the tropical UTLS and constitutes the ascending branch of the meridional residual circulation, also known as the Brewer-Dobson circulation (e.g. Brewer 1949;Holton et al. 1995;Plumb 2002). This circulation can be differentiated in a shallow branch with maximum upwelling in the lower stratosphere and downwelling from subtropical to middle latitudes, and a deep branch whose upwelling extends also along the upper stratosphere and downwelling takes place from middle to polar latitudes (Plumb 2002;Birner and Bönisch 2011). Since the residual circulation is mainly wave-driven, it becomes maximum during winter, when waves can propagate into the stratosphere, and particularly in the more dynamically active Northern Hemisphere (Charney and Drazin 1961). Given the major role of the Brewer-Dobson circulation in global air mass transport, it is a key driver of changes in ozone concentration and radiative heating.
The high-latitude top-down signal is mainly characterised by the strength of the stratospheric polar vortex (Perlwitz and Graff 1995), and the largest tropospheric impact from the stratosphere is expected after extreme events: anomalous vortex intensifications (Limpasuvan et al. 2005;Baldwin et al. 2001) and primarily, sudden stratospheric warmings (SSWs; see Baldwin et al. 2021 for review). During winter, wave activity from the lower stratosphere propagates upwards through westerly background flow while increasing in amplitude until they break. This can alter the vortex strength and eventually trigger SSWs (Matsuno 1971;Charlton and Polvani 2007), although pre-existing conditions in the stratosphere also play a role in the evolution of the vortex (Chen and Robinson 1992;de la Cámara et al. 2017). In fact, the interaction between wave propagation and the background flow is a complex paradigm since parameterization consists of a zonally-symmetric momentum flux that is continuously launched from the mid-troposphere to simulate the effect of gravity waves, and it has been modified from the original implementation in cycle 36r4 (Orr et al. 2010) to be resolution-dependent (see Davini et al. 2017 for details). Orographic gravity waves are parameterized following Lott and Miller (1997) for blocking and hydrostatic gravity wave drag and Beljaars et al. (2004) for turbulent orographic form drag (see ECMWF 2010 for details). Although Döscher et al. (2022) provide information of how stratospheric composition and chemistry are configured in EC-EARTH3.3, neither dynamics nor climatological aspects of the stratospheric circulation are reported. At present this gap calls for attention now that the stratosphere is properly resolved in EC-EARTH.
EC-EARTH3.3 contributes to most CMIP6-endorsed MIPs 4 , particularly relevant here -the Dynamics and Variability MIP (Gerber and Manzini 2016) and the Polar Amplification MIP (Smith et al. 2019), but it also participates in other international initiatives such as the Multi-Model Large Ensemble Archive (Deser et al. 2020) or the SPARC's QBOi. With the focus on the model performance of stratospheric climatology, we make use of the time-slice "Experiment 2" from QBOi phase-I. It consists in a 100year long atmosphere-only simulation with prescribed, repeated seasonal cycle for sea surface temperature (SST) and sea ice concentration (SIC), whose climatology is computed over 1988-2007, and fixed radiative forcing, aerosols and chemical constituents (including ozone) at year 2002, which is well removed from any explosive volcanic eruption and ENSO was in a neutral phase; hence, there is no interannual variability or secular change in the applied forcings (see Butchart et al. 2018). Note that both boundary conditions and radiative forcing/atmospheric composition are taken from the CMIP6 input dataset, instead of CMIP5 as in the original QBOi protocol. Also, it is worth noting that this "Experiment 2" represents the control simulation of the QBOi El Niño/La Niña sensitivity experiments (Kawatani et al., in preparation).
Outputs from QBOi simulations are stored at JASMIN, provided by the British Atmospheric Data Centre (BADC). Here, monthly means of zonal wind (u), air temperature (T) and geopotential height at 28 pressure levels from 1000 to 1 hPa are analysed. Daily-mean variables are only available in a subset of pressure levels; here 300,250,200,100,70,50,30,20 and 10 hPa are analysed. The latter have been used to compute the residual (mass) streamfunction in the Transformed Eulerian-Mean (TEM) framework, which in pressure coordinates can be defined as follows (e.g. Andrews et al. 1987;Abalos et al. 2015): simulation for the last exercise of QBOi Phase-I (i.e. ENSO experiments; Kawatani et al. in preparation) and will be the benchmark of Phase-II. The present analysis aims to contribute to a further understanding of common GCM biases in the stratosphere, and to serve as a reference for future multi-model analyses where EC-EARTH will be taking part.

Model, simulation, and methods
This study focuses on the atmospheric component of the European Consortium coupled climate model EC-EARTH, which is based on the Integrated Forecast System (IFS) atmosphere model of the European Centre for Medium-Range Weather Forecasts (ECMWF). Details about the other main components of the CMIP6 version of EC-EARTH -EC-EARTH3.3, i.e. ocean (Nucleus for European Modelling of the Ocean -NEMO3.6), land (Hydrology Tiled Scheme for Surface Exchanges over Land -H-TESSEL) and sea-ice (Louvain-la-Neuve Sea Ice Model -LIM3), as well as the coupler (Ocean Atmosphere Sea Ice Soil -OASIS3-MCT), can be found in Haarsma et al. (2020; developed for HighResMIP) and Döscher et al. (2022). The standard configuration of IFS in EC-EARTH3.3, based on cycle 36r4, is at horizontal spectral resolution T255 (triangular truncation at wavenumber 255), corresponding approximately to 0.7° in longitudelatitude (~ 80 km), with 91 vertical levels up to 0.01 hPa (L91); hence, the model is considered as "high top". For comparison, the standard configuration of IFS in the CMIP5 version of EC-EARTH -EC-EARTH2.3 (Hazeleger et al. 2010(Hazeleger et al. , 2012, based on cycle 31r1, was at T159, approximately 1.125° in longitude-latitude (~ 125 km), with 62 vertical levels up to 5 hPa; thus, considered as a "low top" model. Developed at ECMWF as a weather model, IFS is tuned and improved for climate purposes by the EC-Earth Consortium (see Döscher et al. 2022 for review). Examples of major advances from CMIP5 to CMIP6 follow. While EC-EARTH2.3 only accounted for the direct effects of prescribed aerosol concentration, EC-EARTH3.3 also includes a representation of the indirect effects via the aerosol impact on clouds (see Wyser et al. 2020 for details). Note that interactive chemistry (particularly, stratospheric chemistry) is not yet implemented in the standard configuration of EC-EARTH (van Noije et al. 2021); so there is no coupling between chemistry and model dynamics/physics. Now, EC-EARTH includes a non-orographic gravity wave (NGW) parameterization, formulated by Scinocca (2003), which, in contrast to previous IFS cycles (before 35r3) that used Rayleigh friction, can spontaneously generate QBOlike variability (e.g. Christiansen et al. 2016); this NGW of ~ 6 K at middle-upper levels. Note that this thermal bias is persistent throughout the seasonal cycle (see Supplementary Material). The temperature in the stratosphere mainly relies on radiative process (e.g. Vallis 2017), hence it is likely due to a problem with the heating linked to ozone absorption between 20 and 50 km (e.g. Edwards 1982), as reported previously for ECMWF-IFS (Hogan et al. 2017).
To discard any possible effect of comparing a simulation with fixed radiative forcing (see Sect. 2) against the ERA-Interim reanalysis with time-varying radiative forcing, the model bias has also been assessed using CMIP6 AMIP runs ( Fig. 11 in Appendix 3), which shows the same issue with the tropical stratospheric temperature. Note that the differences in ozone concentration between EC-EARTH, prescribed from CMIP6 input dataset, and ERA-Interim climatology are minor and inconsistent with the warm bias in Fig. 1f (see Supplementary Material), pointing at the ozone radiation scheme as the cause of the model systematic error (see discussion in Sect. 4).
where a is the Earth radius, ϕ is latitude, g is the gravity, p is pressure, H is the scale-height (7 km), and S = HN 2 /R is the stability parameter, with N 2 the squared Brunt-Väisälä frequency and R the gas constant for dry air (287m 2 /s 2 K).
[v] is the Eulerian-mean meridional wind and [v*T*] the eddy heat flux, where * indicates perturbation from the zonal-mean and [] denotes a zonal-mean. Note that [v*T*] is proportional to the vertical component of the Eliassen-Palm flux, thereby related to vertical propagation of wave activity (e.g. Andrews et al. 1987;Vallis 2017).
The analysis is focused on boreal winter, namely December-February (DJF). When indicated, the model performance is compared to ERA-Interim reanalysis (Dee et al. 2011) covering the period 1979-2019. The statistical significance of the model biases is assessed with a Student's t-test for equal means at 95% confidence level considering each winter as an independent sample. Figure 1 shows vertical cross-sections of zonal-mean zonal wind (left column) and temperature (right column) in boreal winter, illustrating that EC-EARTH (first row) properly simulates the main features of the observed circulation (ERA-Interim; second row); although clear model biases are present (bottom row). In the troposphere, the subtropical jet is slightly shifted to the north (Fig. 1e). In the stratosphere, which is the target of the study, there are two different behaviours. EC-EARTH largely overestimates the climatological westerly wind at upper-stratospheric midlatitudes, where the core of the polar vortex is stronger (~ 6 m/s) and shifted poleward as compared to reanalysis (Fig. 1, left column). Linked to this bias, the meridional extension of the warm tropical stratopause is also biased towards higher latitudes (~ 6 K; Fig. 1, right column). On the other hand, EC-EARTH underestimates the climatological westerly wind (~ 3 m/s) in the lower stratosphere at polar latitudes ( Fig. 1e), consistent with a warm bias (~ 2-4 K; Fig. 1f) in thermal wind balance. The model bias in the strength of the polar vortex at these low levels may affect the vertical propagation of planetary-scale waves from the troposphere into the stratosphere.

Zonal-mean structure
The most striking aspect of the model bias in zonal-mean temperature is probably at low latitudes (Fig. 1f). There appears to be an important issue with the vertical distribution of tropical stratospheric temperature, with a warm bias At 30 hPa, the climatological flow encircles hemispherically at higher latitudes (Fig. 2d,e). The model underestimation of the polar vortex strength (Fig. 1e) is apparent over both the Eurasian and northern North Pacific regions (Fig. 2f). At 3 hPa, the climatological westerly wind is stronger (Fig. 2a,b), and the model overestimation of the polar vortex strength at midlatitudes is especially noticeable over the eastern North Pacific and the North Atlantic-European region.

Brewer-Dobson circulation
The Brewer-Dobson circulation is a meridional overturning, residual mass circulation in the stratosphere that transports air from the tropical tropopause, linked to tropical upwelling, towards extratropical latitudes, including polar downwelling (see Shepherd 2007 for review). Due to limited availability of daily data (see Sect. 2), the analysis can only assess the Brewer-Dobson circulation at middle and lower levels. EC-EARTH properly simulates the observed upward and downward motions of the meridional circulation ( Fig. 3a,b), although model biases are found. In particular, the strength of the residual circulation is underestimated from low to midlatitudes and overestimated at subpolar latitudes (Fig. 3c). Separating the residual streamfunction into the Eulerian-mean component (Fig. 3, middle row) and

Polar vortex
In the extratropical stratosphere, vertical motions are relatively slow (much slower than in the troposphere), the stratification is high and the Rossby number small, thereby the circulation is well described by quasi-geostrophic approximation and the flow can be considered of large-scale and quasi-horizontal (see Andrews et al. 1987, Vallis 2017 for a review on geostrophic scaling). To better characterise the model bias in zonal wind (Fig. 1e), Fig. 2 shows longitude-latitude maps at three vertical levels: 100 hPa (bottom), 30 hPa (middle) and 3 hPa (top); representative of the lower, middle and upper stratosphere, respectively.
At 100 hPa, reminiscent of the upper-tropospheric circulation ( Fig. 2 g,h), EC-EARTH depicts a strengthening (weakening) at the poleward (equatorward) side along the North African-Asian jet and in the North Pacific (Fig. 2i), implying a northward shift of the subtropical jet as compared to ERA-Interim, which is consistent with the zonalmean analysis (Fig. 1e). In the North Atlantic, instead, EC-EARTH yields a weaker and, particularly, less tilted eddy-driven jet than reanalysis (Fig. 2i) Andrews et al. 1987). EC-EARTH accurately simulates the magnitude and extension of [v*T*] in the lower-middle stratosphere (Fig. 4a,b). It shows a slight northward shift from 50 to 10 hPa, as compared to reanalysis (Fig. 4c), but is not statistically significant due to the large variability at those levels. EC-EARTH overestimates the eddy heat flux around 100 hPa at high latitudes, which may counteract the Eulerian-mean component of the meridional circulation (Fig. 3f), contributing to the biased downwelling (Fig. 3c). There is also a positive model bias in the upper-tropospheric mid-latitudes (Fig. 4c), consistent with the biased jetstream (Fig. 1e).

Wave activity
The total wave injection from the troposphere into the stratosphere is estimated as the zonally-averaged 100 hPa the eddy component (not shown), reveals that the former dominates the model bias (Fig. 3f): both, the direct cell in the tropical-subtropical stratosphere and the indirect cell at mid-latitudes are weaker than in reanalysis (Fig. 3d,e). The model bias can be traced back to an underestimation of the zonal-mean meridional wind ([v]; Fig. 3 g,h), particularly the poleward flow in the subtropics, with negative difference values from 10 hPa to the tropopause (Fig. 3i). Note that EC-EARTH also shows a marked bias of the divergent flow in the equatorial upper troposphere-lower stratosphere (Fig. 3c,f,i). These biases in the Brewer-Dobson circulation are consistent with the biases in temperature as discussed below (Sect. 4).  Note the irregular contours in panels a,b,d,e (± 200, ± 150, ±100, ± 75, ±50, ± 30, ±20, ± 10, ±5, ± 3, ±1) and c,f (± 100, ± 50, ±30, ± 10, ±5, ± 3, ±1); the dotted contours in panels g,h correspond to ± 0.1 m/s. Statistically significant areas at 95% confidence level are shaded in panels c,f,i Scandinavia and Siberia (cf. Figure 7 g,h, Fig. 6a,b). Note that both large-scale waves (WN1-2) and small-scale waves (WN3-4) contribute to the EC-EARTH overestimation of the wave forcing in early winter (Fig. 5c). WN5 mainly adds to the model overestimation of v*T* climatology at v*T* over 45-75 N ( Fig. 5a; e.g. Nishii et al. 2009;Smith and Kushner 2012), and compared to its high-latitude contribution averaged over 60-75 N (Fig. 5b); the latter is also considered because stratospheric planetary-scale wave propagation maximizes at high latitudes (Shaw and Perlwitz 2013;Shaw 2015, 2018). The seasonal cycle of the climatological wave activity flux is overall well simulated by EC-EARTH (blue line), although it is overestimated, as compared to reanalysis (black line), in early winter, December/mid-January, and slightly underestimated in late winter, mid-January/February (Fig. 5a). The former appears to be associated with an overestimation of the high-latitude eddy heat flux (Fig. 5b), in agreement with the analysis of [v*T*] (Fig. 4c). This overestimation of the climatological wave activity in early winter could explain the higher frequency of SSW occurrence in the model at this time of the year (see Fig. 9 in Appendix 1), given the potential relationship between the seasonality of these two fields . Figure 6 shows the climatology of v*T*, whose spatial distribution depicts the key regions of tropospheric wave driving, linked to the climatological stationary wave pattern (e.g. Plumb 1985): poleward warm air advection over eastern Eurasia-northwestern North Pacific and equatorward cold air advection over central Eurasia-northeastern Europe; note a weak southward advection of warmer air over northern Canada (e.g. Newman and Nash 2000). EC-EARTH correctly simulates the two areas contributing to positive zonal-mean eddy heat flux (Fig. 6a,b; cf. Figures 4a and  b and 5a), although some differences are found (Fig. 6c). In particular, while the model simulates a weaker and eastward-shifted centre of action around the Aleutian Islands, it yields a weaker but westward-shifted centre of action in the North Atlantic-European sector. On the other hand, EC-EARTH overestimates the amplitude of v*T* over central Eurasia (Fig. 6c).
Decomposing the eddy heat flux into wavenumbers (WNs; Fig. 7) provides valuable information on the waves that contribute to the wave injection. Long planetary-scale waves (WN1 and WN2) are associated with most of the climatological v*T* over the North Pacific and Eurasia (Fig. 7a,b), and represent the largest contributors to the model bias around the Aleutian Islands (Fig. 7c). Interestingly, WN3 contributes to the amplitude of the wave driving in the North Pacific sector (Fig. 7d,e) and shapes the tilt of the model bias toward lower latitudes (Fig. 7f). WN3 appears also to elongate the centre of action of the climatological eddy heat flux over Eurasia (Fig. 7d,e), better defining the shift bias in the North Atlantic and the amplitude bias over central Eurasia (Fig. 7f). On its part, WN4 seems to be the key contributor to the characteristic double lobe of the wave injection over Eurasia, with maxima over . The high-latitude eddy heat flux has been decomposed into WN1-2 (dashed) and WN1-4 (solid) in panel c. Time-series are smoothed with a 7-day running-mean. Statistically significant differences at 95% confidence level are indicated with stars in both time-series deviations of the geopotential height field decomposed into wavenumbers. More evident for large-scale waves (WN1, WN2; first, second row) than for small-scale waves (WN3, WN4; third, fourth row), they all depict a westward tilt with height associated with upward vertical propagation (e.g. Andrews et al. 1987). In the stratosphere, WN1 is the 100 hPa over central Eurasia (Fig. 7 L), although its contribution to the stratospheric circulation is expected to be minor due to the Charney-Drazin wave filtering in the vertical propagation.
Finally, the vertical structure of the climatological stationary waves is assessed. Figure 8 shows zonal-mean Fig. 7 As Fig. 6, but after decomposing v* and T* into wavenumbers (WNs): 1-2 (a-c), 1-3 (d-f), 1-4 (g-i) and 1-5 (jl). Statistically significant areas at 95% confidence level are shaded in panels c,f,i,l Fig. 6 DJF climatology of 100-hPa eddy heat flux (v*T*; m·K/s) in a EC-EARTH and b ERA-Interim, and its difference (c). Statistically significant areas at 95% confidence level are shaded in panel c that EC-EARTH simulates it with a westward shift in the upper stratosphere, which leads to an apparent (but not statistically significant) turning of the vertical tilt (Fig. 8c). On the contrary, the model bias in the lower stratosphere displays a statistically significant eastward shift (Fig. 8,  first row), which is consistent with the eastward shift of v*T* around the Aleutian Islands noted above (Figs. 6c  and 7c).
The model bias of WN2 mainly relies on an amplitude underestimation in the lower stratosphere (Fig. 8f). The wrong performance of WN3's amplitude is maximum in only wave that EC-EARTH simulates with correct amplitude (Fig. 8a,b), yet weaker than ERA-Interim (Fig. 8c). Compared to reanalysis, the model systematically underestimates the amplitude of WN2 (Fig. 8d,e), WN3 (Fig. 8 g,h) and WN4 (Fig. 8j,k).
A close inspection to climatological WN1 (Fig. 8a,b), where the positive zonal-eddy component projects on the stratospheric Aleutian High and the negative one on the displaced polar vortex (e.g. Nigam and DeWeaver 2003) forming a couplet that propagates into the mesosphere (Harvey and Hitchman 1996;Harvey et al. 2002), reveals Fig. 8 Longitude-pressure cross-section of DJF geopotential height (m) climatology in EC-EARTH (left) and ERA-Interim (middle), and its difference (right), decomposed into wavenumbers (WNs): 1 (a-c), 2 (d-f), 3 (g-i) and 4 (j-l). Statistically significant areas at 95% confidence level are shaded in panels c,f,i,l cycle), according to recent progress with ECMWF-IFS (Hogan et al. 2017;Shepherd et al. 2018), which will be incorporated in the next generation of EC-EARTH (version 4) as it will use IFS cycle 43r3 where these improvements became operational. The model bias in the low-latitude stratospheric temperature has also an impact on the zonalmean meridional wind ([v]; Fig. 3i) via the Coriolis term of the momentum equation in the meridional direction, that involves [u]. In particular, positively biased [u] is associated with a weaker [v] in the subtropical low-to-middle stratosphere, which in turn leads to a weaker Brewer-Dobson circulation (Fig. 3c,f). While the residual circulation is mainly driven by the torque induced by wave dissipation, recent work has highlighted the important role of anomalous heating in forcing changes in the Brewer-Dobson circulation (Ming et al. 2016). We argue that this mechanism could be contributing to the residual circulation bias in EC-EARTH. Indeed, the negative bias in tropical upwelling throughout the low-to-middle stratosphere is consistent with the positive bias in the background static stability (Fig. 1f), linked to the biased radiative heating by ozone absorption (as shown for ECMWF-IFS; Hogan and Bozzo 2018), via the TEM thermodynamic equation (e.g. Andrews et al. 1987). On the other hand, resolved wave drag does not feature large biases (see small eddy heat flux bias in Fig. 4c), but parameterized wave drag could contribute to the bias in residual circulation. This contribution might be assessed in EC-EARTH4 (CMIP7), which will presumably have a much weaker thermal bias linked to ozone heating rate.
EC-EARTH also shows a warm bias in the lower stratosphere over the polar cap (Fig. 1f), but, in this case, it does not appear related to radiative processes of ozone (Hogan et al. 2017;Shepherd et al. 2018). Associated with this positive temperature bias, the model simulates a stronger downwelling (Fig. 3c) and a weaker westerly wind of the polar vortex at lower-stratospheric high-latitudes (Fig. 1e), which suggests it is a consequence of too-strong wave driving -resolved (see below) and parameterized (see Appendix 4). EC-EARTH presents a polar cold bias in the lowermost stratosphere (~ 200 hPa; Fig. 1f), a common problem in GCMs that is related to an excessive transport of water vapour from the troposphere and its infrared thermal emission (Gates et al. 1999;Stenke et al. 2008;Hogan et al. 2017), leading to a biased subtropical upper-tropospheric jetstream (Fig. 1e); but this is out of scope.
The amplitude of the climatological planetary waves in the stratosphere is overall underestimated by EC-EARTH at lower levels (~ 10-200 hPa; Fig. 8), but the magnitude of the background wave injection from the troposphere into the stratosphere (i.e. hemispheric-average 100-hPa v*T*) is overestimated (Fig. 4c), particularly in early winter -December-January (Fig. 5). This behaviour could be the lowermost stratosphere (~ 200 hPa) and below the middle stratosphere (Fig. 8i). The amplitude underestimation of WN4 is restricted to the lower-lowermost stratosphere (Fig. 8 L).

Summary and conclusions
This study presents the first comprehensive assessment of boreal winter stratospheric climatology in EC-EARTH, particularly version 3.3 that contributes to CMIP6 and SPARC's Quasi-Biennial Oscillation initiative (QBOi). Indeed, the simulation reported here corresponds to the 100-year control "Experiment 2" of QBOi with prescribed climatological SST and SIC over 1988-2007 and fixed radiative forcing at year 2002; note that the boundary conditions (taken from AMIP) and the radiative forcings are from CMIP6, not from CMIP5 as in the original QBOi protocol (Butchart et al. 2018). Boreal winter stratospheric variability in EC-EARTH was comprehensively assessed by Palmeiro et al. (2020), although with a previous version of the model, EC-EARTH3.1. Appendices 1 (Figs. 9) and 2 ( Fig. 10) report the EC-EARTH3.3 performance of the key stratospheric variability phenomena in this simulation, i.e. sudden stratospheric warmings (SSWs) and the QBO, that are overall realistic. Richter et al. (2020) provide QBO metrics for EC-EARTH3.3 in a CMIP6 multi-model analysis. The modulation of the QBO in EC-EARTH3.3 (Fig. 10) by El Niño and La Niña is currently under investigation in the framework of QBOi (Christiansen et al., Kawatani et al.; in preparation). On the other hand, the influence of El Niño/La Niña on SSWs in a previous version of the model, EC-EARTH3.2, has been recently assessed . Also, inspired by the preliminary results of Haarsma et al. (2020), there are ongoing efforts in the EC-Earth consortium to address the impact of horizontal resolution and radiative forcing on SSWs in the CMIP6 version of the model reported here. In the following, we summarize our main findings on the model climatology in the stratosphere, which may affect variability and predictability (e.g. Hogan et al. 2017;Polichtchouk et al. 2019).
EC-EARTH shows a large issue with the vertical distribution of zonal-mean stratospheric temperature from the tropics to mid-latitudes, associated with a biased warm middle-upper stratosphere (~ 6 K; Fig. 1f). Linked to the enhanced meridional temperature gradient, and in thermal wind balance, the model simulates a stronger and northward-shifted core of the polar vortex at upper-stratospheric mid-latitudes, as compared to reanalysis ([u]; Fig. 1e). This issue with the zonal-mean temperature appears related to ozone heating rate and can be alleviated by modifying the ozone radiation scheme (e.g. timestep, zenith angle, diurnal variability (Hardiman et al. 2017;Haase and Matthes 2019). Implementing fully-interactive stratospheric chemistry is indeed one of the current strategies in climate modelling (Collins et al. 2017;Morgenstern et al. 2017).
A prospect note follows. Good representation of stratospheric circulation and stratosphere-troposphere coupling is expected to provide predictability in subseasonal and seasonal climate forecasting, and careful assessments of the dynamics involved in stratospheric variability of GCM-based forecast systems are underway, e.g. in the tropical stratosphere ) and in the polar stratosphere (Portal et al. 2022). The results shown here recommend that the dynamics involved in stratospheric climatology should also be assessed in forecast mode, for example in the Eulerian-mean [u] (linked to extratropical wave driving) and [v] (related to the Brewer-Dobson circulation), since biases may affect predictions once the models are initialised (see Saurral et al. 2021).

Appendix 1. SSWs in EC-EARTH3.3
See Fig. 9 Fig. 9 November to March climatological distribution of SSWs per decade in a [-10,10]-day window around the SSW date in EC-EARTH (blue) and ERA-Interim (black). Time-series are smoothed with a 7-day running-mean. Note that there is no statistically significant difference between the two curves at 95% confidence level explained by the negative model bias in zonal-mean zonal wind in the lower stratosphere ([u]; Fig. 1e), since a weaker westerly flow is less effective for the Charney-Drazin wave filtering in the vertical propagation. A contributor to this bias could be the parameterized non-orographic gravity wave driving that might be excessive in this model version (see Appendix 4) and can induce changes in the mean flow (de la Cámara and Lott 2015). Further, the overestimated background wave injection could be linked to the overestimated occurrence of SSWs in early winter (Fig. 9), according to its potential impact on the SSW seasonal cycle (e.g. Palmeiro et al. 2020). Hence, attention is needed in model validation to assess the simulation of realistic climatological eddy heat flux seasonality.
The analysis of the spatial distribution of v*T* climatology has brought new insight into the role of climatological planetary waves on the background wave injection. While large-scale waves (WN1-2) dominate the eddy heat flux over the northern North Pacific, small-scale waves (WN3-4) are responsible for the distinctive doubled-lobe structure of the eddy heat flux over Eurasia (Fig. 7). EC-EARTH properly simulates this climatological feature, although overestimates its amplitude over central Eurasia (Fig. 6c). Small-scale planetary wave activity over Eurasia is prominent during wintertime, at seasonal (e.g. Liu et al. 2014;Smoliak and Wallace 2015), intraseasonal (e.g. Bueh and Nakamura 2007;García-Serrano et al. 2017) and submonthly (e.g. Palmeiro et al. 2020Palmeiro et al. , 2022 timescales. It is also worth stressing that small-scale planetary waves propagate vertically reaching middle-upper stratospheric levels (Fig. 8), in agreement with the teleconnection pathway of Eurasian wave activity and the NAO via changes in the polar vortex strength (e.g. Kuroda and Kodera 1999;Takaya and Nakamura 2008).
Recommendations also emerge from this study. Although EC-EARTH spontaneously simulates QBO-like variability, its NGW parameterization follows a source-unrelated approach (see Sect. 2). The absence of source mechanisms in NGW parameterizations limits their potential calibration with the growing number of in-situ and satellite observations and is a cause of systematic errors (de la Cámara et al. 2016). Implementing a source-related parameterization, where the amplitude of NGWs depends on the resolved dynamics in the model, may help to reduce biases in the tropical and extratropical stratospheric climatology and variability (Berner et al. 2017). Another way to alleviate systematic errors in EC-EARTH may be to add chemistry coupling, thereby providing a suitable link between chemical reactions and stratospheric temperature and circulation via the radiative heating associated with the constituents. It has been shown that chemistry-climate models substantially reduce biases in stratospheric transport, climatology and

Appendix 3. Zonal-mean temperature in transient simulations
See Fig. 11.   Fig. 11 Vertical cross-section of DJF zonal-mean temperature (K) climatology in EC-EARTH (top) and its difference with respect to ERA-Interim (bottom); to be compared with Fig. 1b,

Appendix 4. From EC-EARTH3.1 to EC-EARTH3.3
To complement Palmeiro et al. (2020), who focused on boreal winter variability in a previous version of the model, EC-EARTH3.1. Figure 12 shows DJF climatology of zonalmean zonal wind and temperature in EC-EARTH3.1 using a similar model configuration (see Sect. 2 in Palmeiro et al. 2020 for details), as well as their difference with respect to ERA-Interim. The warm bias in the tropical-midlatitude stratosphere is a persistent feature in both model versions, which is largely due to the ozone radiation scheme in the atmospheric component of EC-EARTH, IFS (Hogan et al. 2017;Shepherd et al. 2018). This thermal bias is expected to be reduced in EC-EARTH4 (see Sect. 4) On the other hand, there has been a clear improvement of the stratospheric polar vortex strength in the middleupper stratosphere, which was largely overestimated in EC-EARTH3.1 (Fig. 12) as compared to EC-EARTH3.3 (Fig. 1e). Associated with this considerably biased westerly wind in EC-EARTH3.1, and in thermal wind balance, there is a large negative temperature bias in the upper stratosphere (Fig. 12), which is mostly absent in EC-EARTH3.3 (Fig. 1f). A too strong/too cold stratospheric polar vortex suggests a lack of dynamical variability, mainly wave activity (e.g. Charlton-Perez et al. 2013;Shaw et al. 2014). A hypothesis for this improvement is the difference of the non-orographic gravity wave (NGW) parameterization, with an increase in the launched momentum flux of these waves from subtropical to polar latitudes in EC-EARTH3.3 with respect to EC-EARTH3.1 ( Fig. 13; see Davini et al. 2017 for details). Acknowledgements This work has been supported by the Spanish GRAVITOCAST project (ERC2018-092835). FMP was partially supported by the Spanish ATLANTE project (PID2019-110234RB-C21). JG-S was supported by the "Ramón y Cajal" programme (RYC-2016-21181). MR was supported by a PREDOCS-UB fellowship (2019 − 258) and the "Ayudas para la Formación de Profesorado Universitario" programme (FPU20/03517). MA acknowledges funding from the Spanish STEADY project (CGL2017-83198-R). Red Española de Supercomputación is acknowledged for awarding computing resources (RES projects AECT-2019-2-0019 and AECT-2020-3-0009). Technical support at BSC (Computational Earth Sciences group) is sincerely acknowledged. The authors are thankful to Yolanda Sola (UB) for her help during the review process and two anonymous reviewers for their comments, which helped to improve the scope of the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Data Availability Observational data of the ERA-Interim reanalysis were retrieved from ECMWF (https://www.ecmwf.int/). Model data from QBOi are stored at JASMIN, provided by the British Atmospheric Data Centre (BADC). Model data from CMIP6 AMIP are published at the Earth System Grid Federation (ESGF).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated 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://creativecommons. org/licenses/by/4.0/.

Fig. 12
Vertical cross-section of DJF zonal-mean zonal wind (m/s; left) and temperature (K; right) climatology in EC-EARTH3.1 (top) and its difference with respect to ERA-Interim (bottom); to be compared with Fig. 1b,f   Fig. 13 Amplitude of the launched momentum flux (Pa) in the nonorographic gravity wave parameterization, depending on latitude at T255 resolution, in EC-EARTH3.1 (red) and EC-EARTH3.3 (blue)