The advective Brewer-Dobson circulation in the ERA5 reanalysis: variability and trends

The stratospheric Brewer-Dobson circulation (BDC) is an important element of climate as it determines the transport and distributions of key radiatively active atmospheric trace gases, which affect the Earth’s radiation budget and surface climate. Here, we evaluate the inter-annual variability and trends of the BDC in the ERA5 reanalysis and inter-compare with the ERA-Interim reanalysis for the 1979–2018 period. We also assess the modulation of the circulation by the Quasi-Biennial 5 Oscillation (QBO) and the El Niño-Southern Oscillation (ENSO), and the forcings of the circulation by the planetary and gravity wave drag. A comparison of ERA5 and ERA-Interim reanalyses shows a very good agreement in the morphology of the BDC and in its structural modulations by the natural variability related to QBO and ENSO. Despite the good agreement in the spatial structure, there are substantial differences in the strength of the BDC and of the natural variability impacts on the BDC between the two reanalyses, particularly in the upper troposphere and lower stratosphere (UTLS), and in the upper 10 stratosphere. Throughout most regions of the stratosphere, the variability and trends of the advective BDC are stronger in the ERA5 reanalysis due to stronger planetary and gravity wave forcings, except in the UTLS below 20 km where the tropical upwelling is about 40 % weaker due to a weaker gravity wave forcings at the equatorial flank of the subtropical jet. In the extratropics, the large-scale downwelling is stronger in ERA5 than in ERA-Interim linked to significant differences in planetary and gravity wave forcings. Analysis of the BDC trend shows a global acceleration of the annual mean residual circulation with 15 an acceleration rate of about 1.5 % per decade at 70 hPa due to the long-term intensification in gravity and planetary wave breaking, consistent with observed and future climate model predicted BDC changes.

hemispheric deep branch of the BDC (e.g., Bönisch et al., 2011;Diallo et al., 2012;Monge-Sanz et al., 2012;Chabrillat et al., 2018;Ploeger et al., 2019), consistent with BDC trends derived from trace gas observations, including CO 2 , SF 6 and N 2 O (Stiller et al., 2012;Hegglin et al., 2014;Ray et al., 2010Ray et al., , 2014Haenel et al., 2015). The southern deep branch is strongly modulated by the changes in ozone depletion (e.g. Lin and Fu, 2013;Polvani et al., 2018), therefore, its changes will be affected by the ozone recovery. These findings are also consistent with the advective BDC trends found in reanalyses by 5 Abalos et al. (2015). An updated time series of BDC changes since 1976 derived from in-situ observations of CO 2 and SF 6 in the northern hemisphere middle stratosphere (Engel et al., 2017) show no significant long-term BDC trend between 27 and 32 km, consistent with the mean age of air trends derived from the ERA-Interim and with some climate models (e.g., Garfinkel et al., 2017), but inconsistent with other known reanalyses (Chabrillat et al., 2018;Ploeger et al., 2019). Linz et al. (2017) reported that many climate models are not yet able to reproduce the strength of the BDC as determined from observations. 10 The fact that there is no clear observational evidence of northern hemispheric deep branch changes challenges both the validity of climate model predictions and observational uncertainties in the northern hemisphere middle stratosphere. According to Eichinger et al. (2019) notable source of uncertainty in reanalyses and for future climate projections lies in the strength of the BDC changes. Robust knowledge of the natural variability of the BDC on seasonal to decadal timescales, in turn, is a prerequisite for the reliable detection and attribution of long-term anthropogenically-forced trends. 15 Commonly used as a basis for comparisons of the predicted BDC changes in climate models (Dietmüller et al., 2017;Polvani et al., 2018), the reanalysis data sets provide the best knowledge of the past and present atmospheric state by combining models with observations. Recently, the European Centre for Medium-Range Weather Forecasts (ECMWF) released its fifth generation of atmospheric reanalysis: the ERA5 global reanalysis (Hersbach et al., 2020). Built to replace the ERA-Interim reanalysis (Dee et al., 2011), this newly available high-resolution reanalysis is expected to be a milestone for meteorological analysis as 20 it includes improvements in the representation of atmospheric processes compared to the previous generations of reanalyses.
Hence, it is of key importance to evaluate the representation and characteristics of the BDC in ERA5 reanalysis as well as its consistency with the ERA-Interim reanalysis.
In the present study, we seek to objectively evaluate the representation of the advective BDC in the ERA5 reanalysis in comparison with the ERA-Interim reanalysis, using different diagnostics, including the residual circulation transit time, the residual 25 vertical velocity, the residual mass stream-function and its modulation by the QBO and ENSO signals (Butchart, 2014;Abalos et al., 2015;Ploeger et al., 2015b, a). We describe the reanalysis data sets and methods used in this study, including the Transformed Eulerian Mean and the statistical approaches, in Section 2. The mean climatology of the upwelling/downwelling mass flux and the inter-annual variability from the ERA5 reanalysis are discussed and compared with the ERA-Interim reanalysis in Section 3.1. The residual stream-function and transit times are shown in Section 3.2. Section 3.3 presents the effects of dif- 30 ferent variability modes on the zonal mean wind, temperatures, residual vertical velocity and mass stream-function estimated using a statistical prediction model. Section 3.4 shows the modulation of the advective BDC by the planetary and gravity wave drag. Section 4 presents the trend in the advective BDC together with the trends in planetary and gravity wave drag. Finally, Section 5 provides further discussions and conclusions.

Description of the reanalyses
The wind and temperature data used in this study are from the ERA5 and the ERA-Interim reanalyses, provided by the ECMWF.
The ERA5 reanalysis is based on the Integrated Forecasting System (IFS) and benefits from a decade of progress in model physics, dynamics, and data assimilation. Compared to ERA-Interim, ERA5 has better temporal resolution (6-hourly versus 1-5 hourly) and better horizontal (80 km versus 31 km) and vertical (60 levels versus 137 levels) resolutions and extends higher into the middle atmosphere (0.1 hPa versus 0.01 hPa or 65 km versus 80 km). The ERA5 has replaced the ERA-Interim reanalysis from 31 August 2019. Due to stratospheric temperature biases for 2000-2006 time period exhibited in the first version of ERA5, the ECMWF has published the ERA5.1 to improve upon the cold bias in the lower stratosphere seen in ERA5 (Simmons et al., 2020). In addition, a warm bias higher up persists for much of the period from 1979 (Hoffmann et al., 2019). Apart from 10 higher spatial and temporal resolution, key improvements of ERA5 compared to ERA-Interim are better ability to resolve features like hurricanes and tropical cyclones, better representation of the tropospheric circulation; and assimilation of data from many recent satellite instruments (Li et al., 2020;Hoffmann et al., 2019). Potential limitations of reanalyses, including ERA5 and ERA-Interim, are non-physical trends and variability due to changes in the observing system such as the introduction of COSMIC radio occultation data in 2006, which affects the variability of temperatures near the tropical tropopause. 15 In this study, the wind and temperature fields from both reanalyses have been interpolated onto a 1 • × 1 • longitude and latitude grid, and are extracted from the analysis available at 6 h interval. The difference between ERA-Interim and ERA5 is that the assimilation cycles start at 06:00 UT and 18:00 UT for ERA5 and at 00:00 UT and 12:00 UT for ERA-Interim. In addition, the dynamical fields in ERA5 are archived hourly as tendencies calculated over one hour, and no longer as accumulation over every 3 h hours as in ERA-Interim. 20

Metric of advective Brewer-Dobson circulation
Commonly used as a proxy for the advective BDC, the Transformed Eulerian Mean (TEM) residual circulation is derived from the standard formula of zonally-averaged zonal momentum equation, in latitude and log-pressure coordinates (φ , z), which is given by Andrews et al. (1987) where a is the Earth's radius, φ is latitude. z = −H · ln( p p s ) is the the log-pressure height (vertical coordinate), H is a constant height defined as the scale height (= R · T s /g) for the log-pressure coordinate taken as 7 km, T s (= 240 K) is chosen as standard reference temperature, R is the gas constant for dry air, P s (= 1013 hPa) is chosen as standard reference pressure, (u, v, w) are respectively the zonal mean, meridional and vertical wind velocity, ℑ corresponds to the total zonal momentum forcing and Ω is the Earth's rotation rate. The Coriolis frequency is f = 2 · Ω · sin φ . Here, and in the following, the zonal and the deviations 30 4 https://doi.org/10.5194/acp-2020-881 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License. from the zonal mean are respectively indicated by an over-bar and a prime. The partial derivative with respect to z (vertical) and φ (meridional) direction are indicated by ∂ ∂ z (.) and ∂ ∂ φ (.) respectively. As the advective BDC is a the Lagrangian mean circulation, the Eulerian mean velocity is not a sufficient and direct diagnostic metric. However, a useful proxy for the Lagrangian mean circulation under time-averaged conditions is provided by the residual mean meridional circulation (v * , w * ) (Andrews andMcIntyre, 1976, 1978;Dunkerton et al., 1978;Holton, 1990). The 5 latitudinal and vertical components of the residual mean meridional circulation (v * , w * ) and the mean mass stream-function of the residual circulation, ψ * (φ , z), are given as follows Introducing the equations (2, 3) into the equation (1) will lead to the TEM momentum equation (Andrews and McIntyre,10 1978) where X urGW D is the residual mean nonconservative forcing unresolved by the model and DF is the normalised EP-flux divergence. The X urGW D contribution consists of parameterized gravity wave drag and further imbalances of the reanalysis momentum budget that are caused by the data assimilation system, and that can be interpreted as another contribution of 15 unresolved gravity waves (e.g., Alexander and Rosenlof, 1996;McLandress et al., 2012;Ern et al., 2014). The density ρ = ρ o · exp(−z/H) is the stratification of the atmosphere. θ is potential temperature (hence, w = −(H/p) · (d p/dt), where p = p o · exp(−z/H) is the pressure). The Eliassen-Palm (EP) flux, F, decreases exponentially with height due to the scaling factor ρ sc = exp(−z/H) = ρ/ρ o . DF can be decomposed into resolved planetary (X PW D ) and resolved gravity (X GW D ) wave drags and is written as following with F = {F φ , F z } the EP-flux vector, with respective latitudinal and vertical components DF represents an important part of the forcings of the mean mass flux circulation from the dissipation of resolved planetary 25 and gravity waves. DF can be decomposed into wave numbers ranging from 1 to 180 following previous studies (Ern et al.,5 https://doi.org/10.5194/acp-2020-881 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License.
2014; Alexander and Rosenlof, 1996). The resolved planetary wave drag, X PW D , is estimated as the integration over zonal wave numbers ranging between 1-20 in equation (5). The total gravity wave drag, X GW D , is the sum of the small contribution of explicitly resolved gravity waves with zonal wave numbers larger than 20 in equation (5), parameterized gravity wave drag and the momentum budget imbalances induced by data assimilation. In the reanalyses, the missing wave drag can be considered equal to the unresolved wave drag by the ECMWF model grid, therefore, contributing to the gravity wave forcings in the 5 momentum equation (eq. 4). Thus, the sum of this missing wave drag and the integrated model-resolved wave drag between the 21 and 180 zonal wave numbers gives the estimate of the total gravity wave drag (Alexander and Rosenlof, 1996). For additional details please see Ern et al. (2014Ern et al. ( , 2015Ern et al. ( , 2016. If we consider that X urGW = 0 and DF = 0, representing the zonal forcing due the wave activity such as planetary and gravity waves, then, to maintain the steady conditions, the BDC has to stay non-zero Andrews andMcIntyre (1976, 1978); Eliassen 10 and Palm (1961); Charney and Drazin (1961). Hence under this approximation, the circulation resulting from the wave forcings is the BDC. Further, the mass stream-function of the residual circulation, ψ * (φ , z), is obtained by integrating vertically the total forcing term from the equation (4). In the steady cases and under the downward control principle (e.g. Haynes et al., 1991;Rosenlof and Holton, 1993;Garcia and Boville, 1994), it can be written as 15 where the angular momentum per unit mass is defined by m(φ ) = a · cos φ · (u + a · Ω · cos φ ). At a line, where angular momentum is constant, the mass stream-function is integrated along a contour φ (z).
The tropical upwelling mass flux F is given by This upward tropical mean mass flux, that is the mass per unit time, is integrated for the whole set of latitude circles where 20 the vertical residual component is positive (i.e. w * > 0).
In addition to the mean mass flux and the residual mean mass stream-function diagnostics, we also calculated the residual circulation transit time (RCTT) using a two dimensional backward trajectory model driven by the residual mean meridional circulation (v * , w * ). Used a metric of the advective BDC strength, the RCTT is defined as transit time of an air parcel that is only advected by the residual circulation through the stratosphere. For additional details about the method of the RCTT 25 estimates please see Ploeger et al. (2015b).

Statistical prediction model
For appropriately quantifying the QBO-and ENSO-induced variability in the advective BDC, different metrics of the circulation, including the temperature, zonal wind, w * and ψ * are analyzed using an established statistical prediction model. This regression analysis model has been described in details and applied in our numerous previous studies (Diallo et al., 2012(Diallo et al., , 2017(Diallo et al., , 2018Tao et al., 2019). The statistical model is based on the principle that the monthly zonal mean time-series of the BDC metrics can be decomposed into a sum of different contributions, including a linear trend, QBO, ENSO, seasonal cycle and a residual. For a given metric of the advective BDC strength, BDC metric (herein U, T, w * , ψ * and wave drags), at any latitude and altitude position in the stratosphere, the regression model can be written simply as following The estimated QBO and ENSO amplitude variations are normalized by the standard deviation of the QBO and ENSO predictors, which the QBO index at 50 hPa and the Multivariate ENSO Index (Wolter and Timlin, 2011). The statistical prediction model estimates an uncertainty using a Student's t-test (Zwiers and von Storch, 1995;Friston et al., 2007). Note that the contributions of volcanic aerosol and solar cycle are voluntarily omitted for simplifying the description as these terms are not 10 diagnosed in this study.
3 Climatological advective BDC and its modulations

Inter-annual variability of the Upwelling
Evaluating the large-scale ascent and descent of air mass flux is a prerequisite for disclosing possible biases in the morphology and strength of the advective BDC. Therefore, we compare the tropical upwelling and extratropical downwelling of the air 15 mass flux from the ERA5 to the ERA-Interim using the residual vertical velocity (w * ). Figure 1a-i show the annual mean and seasonal variations of the w * estimated from the two reanalyses, together with the associated differences for the 1979-2018 time period. Overall, there is a remarkably good agreement between the two reanalyses in the main structure of tropical upwelling and extratropical downwelling derived from the annual and seasonal mean w * climatology ( Fig. 1a, b). However, substantial differences in the strength of the tropical upwelling and extratropical downwelling arise in three distinct regions 20 of the stratosphere (Fig. 1c). In the tropical pipe region associated with large-scale upwelling (e.g., Neu and Plumb, 1999;Ray et al., 2014), the ERA5 reanalysis exhibits a weaker annual mean tropical upwelling in the UTLS, consistent with the well-known faster ERA-Interim tropical upwelling than observations (Dee et al., 2011;Seviour et al., 2011;Ploeger et al., 2012). In the mid-latitude surf zone region associated with strong large-scale stirring, and poleward and downward transport (e.g., McIntyre and Palmer, 1984), the ERA5 reanalysis also shows a weaker downwelling. Conversely in the polar vortex 25 region, the ERA5 shows a stronger large-scale downwelling of the air mass than the ERA-Interim, suggesting a stronger advective BDC in ERA5. The maximum difference in the annual mean extratropical downwelling is about 0.5 mm.s −1 and occurs in the southern hemisphere. The residual velocity also exhibits seasonal variations in these three distinct regions of the stratosphere with even larger differences between the two reanalyses in the upper stratosphere above 30 km ( Fig. 1d-i). As the strength of advective BDC varies seasonally between the hemispheres, these large differences in the w * are displaced toward show the annual mean climatology from the ERA5 (a), the ERA-Interim (b) and the associated differences (c). The panels (d-f) show the DJF mean from the ERA5 (d), the ERA-Interim (e) and its associated differences (f). The panels (a-i) show the JJA mean from the ERA5 (g), the ERA-Interim (h) and its associated differences (i). Grey vertical line indicates the zero w * contours from the ERA5 and ERA-Interim, respectively. The climatological tropopause from the ERA-Interim is indicated as the llack dashed horizontal line. Unit: mm · s −1 . with the stronger ERA-Interim upwelling (Fig. 1). Above 20 km, ERA5 reveals stronger tropical upwelling than ERA-Interim.
In the northern hemisphere, the downwelling differences between the two reanalyses are negligible below 30 km and reach 10 30 % above (Fig. 2c, f). The errorbars, which correspond to 15 times the standard deviation from the mean differences, indicate that the large differences in upwelling and downwelling correspond to the regions where the maximum variance occurs. The observed differences in the upwelling and downwelling are likely due to a combination of several factors, including a higher resolution, higher model top, new non-orographic gravity wave parameterization scheme, and the progress made in the more realistic representation of the sponge layer in ERA5 and will be discussed later (Sect. 3.4). 15 The seasonal variations in the w * and its associated differences between the two reanalyses are also assessed at 70 hPa for the 1979-2018 time period (Fig. 3a-c). The 70 hPa pressure level is commonly used as the reference level for model inter-comparisons of the upwelling strength of the air mass entering the stratosphere (e.g., Butchart et al., 2010). The tropical upwelling and extratropical downwelling patterns also agree remarkably well between the two reanalyses even in the variations of the upwelling zero-line. The w * exhibits a 6-month phase shift between the northern and the southern hemisphere mid-20 https://doi.org/10.5194/acp-2020-881 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License. latitudes, resulting from the correlation of the tropical upwelling annual cycle with the strongest descent located in the winter hemisphere ( Fig. 3a, b). At 70 hPa, a seasonal variation of the differences in the w * between the two reanalyses is also shown in Figure 3c. The weaker tropical upwelling in ERA5 varies seasonally between 0 and 20 • N with a maximum of about 0.2 mm . s −1 . In the mid-latitude, the large-scale downwelling is also weaker and varies seasonally with a maximum occurring during the boreal summer (June-July-August). In the polar vortex region, ERA5 shows a stronger large-scale downwelling than 5 ERA-Interim, as already evident from Fig. 1-2. The mass flux differences maximize in the polar vortex region of the southern hemisphere. These seasonal variations in the tropical upwelling and the extratropical downwelling results from the pumping action of the wave forcings (Dunkerton et al., 1981;Randel, 1993;Holton et al., 1995;Plumb and Eluszkiewicz, 1999).

Stratospheric residual circulation
In addition to the w * analyses, we also evaluate the consistency and uncertainty in the residual circulation mass stream-function pipe and the mid-latitude surf zone than in ERA-Interim, but stronger in the polar regions ( Fig. 4a-c). The maximum difference in the residual circulation also occurs in the UTLS below 20 km, consistent with the w * differences. The seasonal variations in the residual circulation also agree well in the structure, but not the strength during the winter and summer between the two reanalyses ( Fig. 4d-i). During the boreal winter (December-January-February), the differences in the mass stream-function between the two reanalyses are negative between the 0 and 20 • N, which extends from the UTLS into the upper stratosphere. 20 This negative difference indicates a weaker circulation in ERA5 and consistent with the weaker upwelling found in Fig. 3.
However, the differences in the mass stream-function between the two reanalyses are positive between the 40 • S and 0 and extend toward the southern hemisphere, leading to a stronger extratropical descent of the residual circulation in ERA5 than in ERA-Interim ( Fig. 4d-f). The residual circulation is clearly stronger in ERA5 than in ERA-Interim. During the boreal summer (June-July-August), these discrepancies in the residual circulation are reversed ( Fig. 4e-i). These dominating seasonal features 25 of the differences in the residual circulation suggest that the significant improvements in ERA5 likely induce a horizontal shift of upwelling in the UTLS compared to ERA-Interim.
For further insights into the circulation differences between the two reanalyses, we also evaluate the annual and seasonal mean variations of the residual circulation transit time (RCTT), which is defined as the integrated time-scale of air mass transport by the pure residual circulation (e.g., Birner and Bönisch, 2011;Garny et al., 2014;Ploeger et al., 2015a, b). Note 30 that the RCTT represents the integrated residual circulation effect, whereas the w * is a local quantity. The RCTT shows a very good agreement in the morphology of the circulation between the ERA5 and ERA-Interim, e.g. the shorter transit time in the tropics (faster BDC) and longer transit time in the extra-tropics (slower BDC) for 2017-2018 period (Fig. 5a, b). However, the RCTT also shows clear differences in the residual circulation, consistent with the discrepancies in residual mass stream-  function and residual vertical velocity previously discussed. In the tropical pipe (e.g., Neu and Plumb, 1999;Ray et al., 2014) and mid-latitude surf zone (e.g., McIntyre and Palmer, 1984), the annual mean RCTT in ERA5 exhibits longer transit time (below about 0.5 year) associated with the slower tropical ascent of the transition and shallow branch of the BDC (Lin and Fu, 2013;. In the northern hemisphere polar regions, the integrated residual circulation from the ERA5 shows a shorter transit time (below about 0.5 year), consistent with a faster large-scale downwelling of the advective BDC in ERA5 5 than in ERA-Interim in the northern hemisphere. In the southern hemisphere, the residence time in ERA5 tends to be longer in ERA5 than in ERA-Interim. Notice that the RCTT is based just on a two year averaged, therefore, the circulation differences might not be representative of the full reanalysis period. However, the RCTT discrepancies are consistent with the differences in w * and ψ * . This remark is even more relevant when analyzing the differences in the seasonal variations, which show almost very similar differences as shown in the annual mean RCTT (Fig. 5d-i). The differences vary seasonally and the maximum 10 difference is found in the southern hemisphere near the polar vortex ( Fig. 5f, i).

Natural variability related to QBO and ENSO
To further understand the linkage of differences between the two reanalyses with the impact of natural modes of climate variability, we analyze the representation of the QBO and ENSO variability. One of the major mode of variability in the ascending branch of the BDC on seasonal to intrannual is the QBO (Lindzen and Holton, 1968;Plumb and Bell, 1982). 15 Composed of alternating westerly and easterly zonal wind shears, the QBO propagates downward from the tropical middle stratosphere into the troposphere with a period of ∼ 28 months. Both reanalyses agree remarkably well in the downward propagating QBO phases (Fig. 6a, b). The depiction of the QBO westerly and easterly phases from the lower to the upper stratosphere (from about 15 to 50 km) in ERA5 is very similar to the ERA-Interim for the 1979-2018 time period. The QBO disruption in January 2016, which was associated with the development of an easterly phase in the center of the westerly 20 phase (Osprey et al., 2016;Newman et al., 2016), is clearly visible in both reanalyses. Apparent differences are also observed in the equatorial transitions in the eastward and westward zonal mean zonal wind and in the strength of QBO westerly and easterly phases between the ERA5 and ERA-Interim reanalyses. ERA5 exhibits stronger westerly and easterly shear of the QBO compared to the ERA-Interim (Fig. 6c). In addition to being stronger, the westerly QBO phase in ERA5 extends further into the troposphere below the tropopause, while the ERA-Interim QBO westerly phase does not propagate that far downward. 25 According to previous findings, the QBO westerly phase is sensitive to many model details, including the parametrized nonorographic gravity wave drag, and the vertical and horizontal resolution Geller et al., 2016;Polichtchouk et al., 2017Polichtchouk et al., , 2018Shepherd et al., 2018). The amplitude of the westerly phase of the QBO and its frequency increase with increasing non-orographic gravity waves. In the ERA5 a non-orographic gravity wave drag parameterization is used, which is different from the ERA-Interim, likely related to the observed differences (Orr et al., 2010;Scaife et al., 2002;Lott et al., 30 2012; Richter et al., 2014). In ERA-Interim, the Rayleigh drag used to be a substitute for the non-orographic gravity wave drag. When a Warner and McIntyre (2001)  as sources of gravity waves are better represented, such as convection, which is important for gravity wave driving the QBO (Lindzen and Fox-Rabinovitz, 1989). is partly driven by waves trapped in the equatorial region (Wallace et al., 1993;Baldwin et al., 2001;Ern and Preusse, 2009;Ern et al., 2014). The QBO modulates the stratospheric residual circulation in the meridional and vertical direction by affecting temperature and tropical upwelling (Plumb and Bell, 1982;Collimore et al., 2003;Niwano et al., 2003;Punge et al., 2009). The QBO westerly shear is partly induced by equatorially trapped Rossby-gravity waves and the QBO easterly shear is induced by Kelvin waves (Holton and Lindzen, 1972;Plumb, 1977). Figure 7a, b show that the QBO westerly amplitude variation 10 occurs between 15 and 25 km, the QBO easterly amplitude variation between 25 and 35 km and again QBO westerly amplitude variation above 35 km. Note that the region of tropical westerly shear between 15 and 20 km is associated with anomalously warm tropical tropopause temperatures below 20 km (Fig. 7c, d) because it reduces the tropical upwelling and enhances the horizontal transport and mixing of stratospheric trace gases poleward, consistent with anomalously cold temperature in the extratropical lower stratosphere (Plumb and Bell, 1982;Trepte and Hitchman, 1992). The strength of the QBO westerly wind and shear amplitude variations near the tropical tropopause region in ERA5 corroborates its weaker tropical upwelling in the UTLS compared to ERA-Interim (Figs. 1-3). Conversely, the region of tropical easterly shear between 20 and 30 km is associated with anomalously cold tropical temperatures between 20 and 30 km, leading to enhanced tropical upwelling by the 5 secondary meridional circulation associated with the QBO (Choi et al., 2002). The tropical upwelling is anti-correlated with the tropical temperature above the tropopause and its strength modulates trace gas mixing ratios by advecting tropospheric air into the stratosphere (Randel et al., 2006;Diallo et al., 2018;Ray et al., 2020). The strength of the QBO modulation of the circulation is stronger in ERA5 than in ERA-Interim as shown in both QBO-induced temperature and zonal mean wind amplitude variations for the 1979-2018 period. Near the polar vortex region, the QBO effect on temperatures in ERA5 is larger 10 than in ERA-Interim likely related to the Holton-Tan effect (e.i. the QBO modulates the Arctic polar vortex) (Holton and Tan, 1980;Garfinkel et al., 2012;Lu et al., 2014). In the upper stratosphere above 30 km, large differences in the strength are also visible in both temperatures and zonal mean winds associated with the differences in the QBO westerly shear and wind regions, consistent with the differences in w * (Fig. 1).
In addition to the analysis of secondary meridional circulation, we analyzed the QBO-and ENSO-induced variability in 15 the stratospheric residual circulation from variations in w * and ψ * . Considered as one of major modes of BDC modulation , the coupled atmosphere-ocean phenomenon, ENSO, is defined as extreme sea surface temperature (SSTs) changes in the tropical Pacific Ocean. These drastic changes in SSTs encompass with severe surface climate and weather conditions (e.g., Bjerknes, 1969;Cagnazzo and Manzini, 2009;Wang et al., 2016). ENSO has two phases known as El Niño (anomalously warm SSTs) and La Niña (anomalously cold SSTs), and its occurrence varies between 2 and 8 years (Philander, 20 1990; Baldwin and O'Sullivan, 1995). El Niño modulates the UTLS by warming the upper troposphere and cooling the tropical lower stratosphere. The latter has been associated with an acceleration of the ascending branch of the BDC (Randel et al., 2009;Calvo et al., 2010;Konopka et al., 2016).
To quantify the QBO-and ENSO-induced variability in the w * and ψ * , the statistical analysis (eq. 10) is performed by explicitly including ENSO and QBO predictors to isolate their modulations of the BDC. The obtained QBO and ENSO ampli- The results of the QBO-induced variability in the w * show a very good agreement between the ERA5 and ERA-Interim reanalyses (Fig. 8a, b). Both reanalyses depict similar patterns of the QBO-induced changes in the w * . The structural changes in w * indicate a clear decrease in the tropical upwelling below 20 km induced by the QBO westerly shear between 15 and 20 km 30 and consistent with QBO-induced anomalously warm tropical temperature changes in tropical UTLS below 20 km (Fig. 7c, d).
Conversely, the tropical upwelling increases between 20 and 30 km due to the QBO easterly shear in that region, consistent with the QBO-induced anomalously cold tropical temperatures between 20 and 30 km (Fig. 3c, d). Another decrease in the tropical upwelling is induced by the QBO westerly shear above 30 km, consistent with the structure of QBO-induced temperature changes (Fig. 7c, d). The w * variations induced by the tropical QBO westerly shear region at altitudes 15-20 km, QBO easterly phase at altitudes 20-30 km and tropical westerly QBO shear above 30 km are consistent with previous findings regarding QBO modulations (Plumb and Bell, 1982;Trepte and Hitchman, 1992;Collimore et al., 2003;Niwano et al., 2003;Punge et al., 2009). Also note that the positive w * anomalies induced by the QBO easterly shear in the tropical middle stratosphere below 30 km, and extending toward the extra-tropics are remarkably well captured in both reanalyses. The QBO westerly shear-induced negative w * anomalies in the middle and upper stratosphere above 30 km agree fairly well in both reanalyses.

5
The two main differences are the strength of the QBO-induced variability with the ERA-Interim, showing weaker region of QBO easterly and westerly effects. In addition, the positive w * anomalies in the upper tropical stratosphere above 40 km are missing in the ERA-Interim.
Our analysis of the QBO-induced variability in the residual stream-function circulation is also consistent with the QBO impact on w * (Fig. 8c, d). The QBO-induced morphological changes in the residual stream-function show a good agreement 20 and 30 km. Also, the residual stream-function shows a stronger QBO modulation in the ERA5 than in the ERA-Interim ( Fig. 8a, b).

5
The ENSO-induced variations in w * agree well between the two reanalyses in their morphology. The structural changes in w * in the UTLS region are characterized by positive anomalies in the tropics and negative anomalies in the extra-tropics (Fig 9a,   b). The El Niño-like condition enhances the tropical upwelling between 15 and 20 km. These ENSO-induced variations in the w * from the two reanalyses agree well. During El Niño-like condition, the strengthening of the tropical upwelling (positive anomalies of the w * in the tropics) increases upward transport of young air from the troposphere into the stratosphere (Calvo period. Unit: mm · s −1 × SD(proxy)). The climatological tropopause from the ERA-Interim is indicated as the black dashed horizontal line.
Zonal mean climatology of w * is overplotted as dashed grey lines.
downward transport of old stratospheric air into the polar regions, competing with the projected acceleration of the BDC in a warming climate (e.g., McLandress and Shepherd, 2009;Lin and Fu, 2013;Butchart, 2014;Hardiman et al., 2014).
However, there are differences in the magnitude of the ENSO-induced variation in w * . The w * changes related to ENSO show a globally stronger tropical upwelling and extratropical downwelling in the ERA5 than in the ERA-Interim, except in the tropical UTLS below 20 km where ERA5 exhibits weaker anomalies than ERA-Interim. This specific difference in the UTLS 5 suggests that because of the high vertical resolution in ERA5, the tropospheric ENSO-induced variation is more confined below the tropopause barrier in ERA5 than in ERA-Interim, consistent with larger inter-annual differences in the reanalyses at levels below the cold point (Tegtmeier et al., 2020). The ERA-Interim w * positive anomalies in the inner tropics extend from the upper troposphere into the lower stratosphere while in ERA5 this is not visible. Additional differences occur in the tropical middle stratosphere between 20 and 30 km where ERA-Interim shows larger negative anomalies than the ERA5. In the northern hemisphere upper stratosphere, ERA5 also reveals larger w * negative anomalies than the ERA-Interim likely due to the difference in wave activities (Randel et al., 2002(Randel et al., , 2008. Insights on the ENSO-induced structural changes in the advective BDC become clearer from the analysis of the residual stream-function (Fig. 9c, d). The structural changes in the residual circulation induced by the ENSO also show a good agreement between the ERA5 and ERA-Interim reanalyses, consistent with the w * changes induced by the ENSO. In both reanalyses, 5 ENSO strengthens the residual circulation of the BDC, particularly the shallow branch (Fig. 9c, d) Yang et al., 2014). The changes in the deep branch are less evident , but stronger in ERA5 than in ERA-Interim. This spatial changes in the residual circulation is consistent with positive (negative) stream function changes in the northern hemisphere (southern hemisphere). The apparent differences between the two reanalyses are associated with the strength of the ENSO-induced acceleration of the BDC. In the inner tropical UTLS below 20 km, ERA5 exhibits weaker circulation anomalies than the ERA-Interim, consistent with the differences between the two reanalyses in the tropical upwelling (Fig. 9a, b). In the middle and upper stratosphere, the strengthening of the residual circulation is stronger in ERA5 than in ERA-Interim. The strengthening of the deep branch in response to ENSO does not extend as far upward in ERA-Interim as it does in ERA5, consistent with a less evident ENSO-induced changes in the deep BDC branch . As all differences in the circulation and in the natural variability point out to the possible discrepancies in the BDC forcings, we analyze the contribution of the planetary 15 and gravity wave in the following (Sect. 3.4).

Planetary and gravity wave forcings
To better understand the contribution of wave forcings to the circulation and the natural variability differences between the ERA5 and ERA-Interim, we evaluate the net forcing, the planetary and the gravity wave driving the BDC. conditions (∂ u/∂t = 0) and negligible meridional and vertical gradients of the zonal wind, and assuming that the right hand side of equation (4) is only given by wave forcing ℑ, this equation reduces to ℑ = − f * v * . This means that the meridional wind contribution δ v * that is attributed to zonal wave forcing can be written as: δ v * = −ℑ/ f . Correspondingly, positive wave drag 25 indicating eastward forcing would weaken the BDC (δ v * is equatorward), and negative wave drag indicating westward forcing would weaken westerly winds, therefore, accelerating the BDC (δ v * is poleward). (Please note that the Coriolis parameter, f , changes its sign at the equator, and this will result in δ v * switching its sign at the equator if the sign of ℑ does not switch).
Overall, the morphology of the mean wave drag climatology agrees well between the ERA5 and ERA-Interim, including the positive wave drags in the easterly zonal mean wind regions and negative wave drag in the westerly zonal mean wind 30 regions (Fig. 10a-f). However, significant differences occur in the amplitude of the wave drag ( Fig. 10g-i). The gravity wave drag contributes the most to the UTLS differences in the net forcings of the BDC between the ERA5 and ERA-Interim, while the planetary and gravity wave contribution is stronger in the upper stratosphere ( Fig. 10g-i). In the UTLS below 20 km, the previously reported weaker tropical upwelling in ERA5 compared to ERA-Interim is due to a weaker gravity wave breaking at the equatorial flank of the subtropical jets in the ERA-Interim (Fig. 10g, i). The contribution of the planetary waves to the tropical upwelling differences is less evident. The breaking of synoptic and small-scale waves near the tropical tropopause 21 https://doi.org/10.5194/acp-2020-881 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License.
layer and the upper flank of the subtropical jets are the primary forcing sources of the transition and shallow branches of the BDC (Plumb, 2002;Shepherd and McLandress, 2011;. Thus, this difference in the strength of gravity wave breaking in the UTLS explains the weaker ERA5 tropical upwelling as observed in the residual vertical velocity, residual mass stream-function and the RCTTs. In the upper stratosphere above 30 km, the differences in the net forcings of the BDC between the two reanalyses are also related to the contribution of the planetary and gravity wave breaking (Fig. 10g-i). While 5 the planetary wave drag is stronger in the extra-tropics between 30 and 40 km in ERA5 than in ERA-Interim, it is weaker above 40 km. In the tropical upper stratosphere above 40 km, the gravity wave drag in ERA5 is significantly weaker than in ERA-Interim. These differences are likely due to better vertical and horizontal resolution, higher model top and changes to the non-orographic gravity wave drag parametrization scheme in ERA5 Polichtchouk et al., 2018). In ERA-Interim the non-orographic gravity wave drag was handled by the Rayleigh drag, but in ERA5 the scheme is based on 10 non-orographic gravity waves parametrization from Warner and McIntyre (2001). So ERA5 does not have Rayleigh drag on any longer. The energy deposition from upward propagating gravity waves is a significant contributor to the thermodynamic budget at these altitudes. In addition, note that the sponge layer is also not always applied equally to different variables, and is sometimes enhanced in the vicinity of the equator presumably to control equatorial inertial instability.
Furthermore, the apparent differences in wave drag between the ERA5 and ERA-Interim reveal seasonal variations and gravity wave dissipation are also consistent with the discrepancies in the BDC between the ERA5 and ERA-Interim. In the 25 UTLS region, the gravity wave breaking at the equatorial flank of the subtropical jet causes the differences between the ERA5 and ERA-Interim tropical upwelling (Fig. A2a-i). The observed maximum of the seasonal variations in the tropical upwelling at 70 hPa between 0 and 20 • N is due to a weaker gravity wave breaking in the ERA5 UTLS. In the upper stratosphere, the gravity wave breaking in the southern hemisphere winter has a stronger contribution than the planetary waves and is also stronger in ERA5 than in ERA-Interim, suggesting a stronger southern hemispheric deep branch (Fig. A2g-i). This is also consistent 30 with the reported hemispheric asymmetry in the w * , RCTT, and natural variability induced variations. As the strength of the BDC trends is also driven by the net wave breaking, we expect the BDC trends to be impacted by the stronger wave drag in the ERA5 reanalysis, there, we estimate the trends in the following (Sect. 4).
For better insights in the impact of stronger gravity and planetary wave drag on the BDC in ERA5, we also investigate the structural changes in the advective BDC and wave drag for the 1979-2018 time period (Fig. 11a-d). The ERA5 reanalysis shows an acceleration of the BDC, consistent with previous reanalysis studies (Diallo et al., 2012;Abalos et al., 2015;Miyazaki et al., 2016) and climate model projections (e.g., Butchart et al., 2010;Hardiman et al., 2014). The annual mean trend of the previous studies (Sigmond and Shepherd, 2014;Abalos et al., 2015;Eichinger et al., 2020). The hemispheric asymmetry in the BDC trends is due to the asymmetry in wave breaking, which is stronger in the northern hemisphere than in the southern hemisphere as indicated by the negative net forcing trend (Fig. 11b). We also note that the southern hemispheric circulation cell exhibits an apparent structural difference compared to the northern hemispheric circulation cell. The southern hemisphere residual circulation trend is weaker than in the northern hemisphere, consistent with the asymmetry in wave breaking trends, 20 which are stronger in the northern hemisphere. A detailed analysis of the BDC trends will be presented in further studies using the stratospheric age of air and age spectrum. Note that the ERA-Interim still shows a negative trend in the residual stream-function (not shown), consistent with the weakening upwelling trend estimated from w * (Seviour et al., 2011).
In addition, the time series of deseasonalized monthly mean tropical upwelling mass flux averaged between the turnaround latitudes show a comparable inter-annual variability between the ERA5 and ERA-Interim reanalyses at the pressure 25 level of 70 hPa (Fig. 12). The estimated annual mean upwelling mass flux is about 6.17 × 10 9 kg · s −1 for ERA5 and about 6.86 × 10 9 kg · s −1 for ERA-Interim for the 1979-2018 period, consistent with the multi-model mean value of 5.9 ×10 9 kg · s −1 from previous studies (e.g, Butchart et al., 2010;Hardiman et al., 2014;Linz et al., 2017). The ERA5 tropical upwelling mass flux trend at 70 hPa calculated from w * is about 1.5 % per decade, consistent with the observed BDC changes of about 1.7 % per decade (Fu et al., 2019) and climate model predictions of an increase in the strength of the BDC of about 30 2 % per decade due to the increasing GHGs in the atmosphere (Butchart et al., 2010;Garny et al., 2011;Hardiman et al., 2014).
While the climate model projected trend is statistically significant, th mass flux trend in ERA5 is not significant. Note that the reported negative long-term trend by Seviour et al. (2011) in the tropical upwelling mass flux estimated from the ERA-Interim w * disappears in ERA5. This negative long-term trend in standard w * calculations is found to be intricately related to ERA-Interim and absent in other reanalyses (Abalos et al., 2015;Miyazaki et al., 2016). As a basis for verification, we also calculated the long-term trend in w * in ERA5. This estimated mean upwelling mass flux trend from the ERA5 reanalysis is about 0.1 × 10 9 kg · s −1 · dec −1 instead of -0.47 × 10 9 kg · s −1 · dec −1 for the ERA-Interim. The relative difference between the two reanalyses in the mass flux is not large, but significant and is between 10 % and 20 % below 70 hPa, indicating that the ERA5 upwelling change is more consistent with observed and climate model predicted strengthening BDC than ERA-Interim.

Summary and Conclusions
In this study, we assess the inter-annual variability and trends of the advective stratospheric BDC in the ERA5 and ERA-Interim reanalyses, using different circulation metrics, including the w * , residual stream-function and RCTT. We also evaluate the impact of the natural variability, including the QBO and the ENSO, on the BDC as well as the impact of planetary and gravity wave drag on the BDC and its changes.

5
Based on our comparisons of the circulation from the ERA5 with the ERA-Interim, we found a remarkably good agreement in the morphology of the advective BDC, including the tropical upwelling and extratropical downwelling as well as in the QBOand ENSO-induced impact on the BDC. Despite the good agreement in the spatial structure, there are significant differences between the ERA5 and the ERA-Interim in the strength of the advective BDC and in its modulation by the natural variability in the UTLS, and in the upper stratosphere. The advective BDC and its induced-modulations by the natural variability are 10 stronger in ERA5 due to stronger planetary and gravity wave drags. However, in the tropical pipe region below 20 km, the tropical upwelling is about 40 % weaker in ERA5 than in ERA-Interim due to a weaker gravity wave breaking at the equatorial flank of the subtropical jet. In the extra-tropics, the large-scale downwelling is stronger in ERA5 than in ERA-Interim linked to significant differences in planetary and gravity wave drag. These differences in ERA5 vary seasonally and are consistent within all metrics used for evaluating the BDC, including the w * , ψ * , RCTT and modulations by the natural variability. In addition, 15 the differences between the ERA5 and ERA-Interim show a hemispheric asymmetry and a seasonal variation with a maximum difference of about 60 % stronger downwelling in ERA5 arising in the southern hemisphere.
Regarding the QBO signal in the zonal mean wind and temperatures, we found a good agreement in the morphology of the secondary meridional circulation between the two reanalyses. The analysis of the QBO amplitude variation revealed anomalously warm tropical tropopause temperatures associated with the QBO westerly shear modulation, which decreases the tropical upwelling between 15 and 20 km, and also between 30 and 50 km in both reanalyses. Conversely, the QBO easterly shear modulation between 20 and 30 km induced anomalously cold tropical temperatures between 20 and 30 km, leading to enhanced 5 tropical upwelling by the QBO secondary circulation. Besides the good structural agreement, noticeable differences are also found in the strength of the QBO easterly and westerly shear modulations. The QBO modulation is stronger in ERA5 than in ERA-Interim. In the UTLS, the QBO westerly shear in ERA5 is stronger, therefore, consistent with a weaker tropical upwelling in ERA5 compared to ERA-Interim. In addition, the QBO westerly phase in ERA5 extends further down in the UTLS region below the tropical tropopause level. 10 The regression analysis of the QBO-and ENSO-induced structural changes in w * and ψ * shows a very good agreement in the modulations of the BDC between the two reanalyses. Similarly to the QBO-induced secondary meridional circulation in the temperatures, both reanalyses agree remarkably well on a weaker upwelling in the UTLS region below 20 km and in the upper stratosphere above 30 km due to the QBO westerly shear and the stronger upwelling in the middle stratosphere induced by the QBO easterly shear. The variations in the residual circulation induced by the ENSO are also very consistent between the 15 two reanalyses. However, significant differences have also been found in the strength of the BDC modulations by the QBO and ENSO. The impact of QBO easterly and westerly shear on the BDC is stronger in ERA5 than in ERA-Interim. Analogously, the ENSO impact on the BDC is globally stronger in ERA5 than in ERA-Interim, except in the UTLS below 20 km where ERA5 seems to be weaker than ERA-Interim, consistent with the weaker ERA5 tropical upwelling.
The analysis of the planetary and gravity wave drag clearly shows that the main discrepancies in the BDC in the UTLS 20 region between the ERA5 and ERA-Interim, including the tropical upwelling and mid-latitude downwelling, are mainly due to a weaker gravity wave breaking at the equatorial flank of the subtropical jets in ERA5. The discrepancy in the planetary wave drag between the ERA5 and ERA-Interim is less evident in the UTLS region. In the upper stratosphere above 30 km, the differences in gravity and planetary wave drag between the two reanalyses are significantly large, consistent with a weaker BDC in the UTLS and stronger BDC in the middle and upper stratosphere in ERA5. These differences in wave forcings, 25 therefore, in the BDC strength and its modulations, are very likely related to a combination of several factors, including a higher resolution, higher model top, new non-orographic gravity wave parameterization scheme, and the progress made in the more realistic representation of the sponge layer in ERA5, which is presumably designed to damp upward propagating waves near the model top to avoid reflections from the upper boundary.
Finally, the estimates of the advective BDC trend in ERA5 show a global acceleration of the annual mean residual circulation the trends of the net wave forcing of the BDC. This hemispheric asymmetry in the wave breaking trends is linked to more disturbed winter polar vortex in the northern hemisphere than in the southern hemisphere as well as to ozone recovery in the southern polar vortex (Polvani et al., 2018). Data availability. The ERA5 (https://apps.ecmwf.int/data-catalogues/era5/?class=ea last access: 20 August 2020) and ERA-Interim data (https://apps.ecmwf.int/data-catalogues/era-interim/?class=ei, last access: 20 August 2020) are available. The advective BDC and ERA5 data used for this paper may be requested from the corresponding author (m.diallo@fz-juelich.de). Author contributions. MD designed the study, conducted research, performed the calculation and the complete analysis of the advective BDC diagnostics as well as drafted the first manuscript. FP calculated the RCTTs. ME calculated the wave decomposition. FP and ME