The advective Brewer-Dobson circulation in the ERA5 reanalysis: climatology, 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, climatology and trends of the BDC in the ERA5 reanalysis and inter-compare with its predecessor, the ERA-Interim reanalysis, for the 1979–2018 period. We also assess the modulation of the circulation 5 by the Quasi-Biennial Oscillation (QBO) and the El Niño-Southern Oscillation (ENSO), and the forcings of the circulation by the planetary and gravity wave drag. The 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 and signiﬁcant differences in the strength of the BDC and natural variability impacts on the BDC between the two reanalyses, particularly in the upper troposphere and lower stratosphere 10 (UTLS), and in the upper 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 up to 40 % weaker mainly due to a signiﬁcantly weaker gravity wave forcing at the equatorial-ward of the upper ﬂank of the subtropical jet. In the extra-tropics, the large-scale downwelling is stronger in ERA5 than in ERA-Interim pressure 100-70 hPa) and the shallow branch (at pressure 70-30 hPa). The two-way mixing is deﬁned as a quasi-horizontal stirring and irreversible dis- 20 placement of air masses induced by breaking of large- and small-scale waves in the surf zone, and turbulent mixing (McIntyre and Palmer, 1984; Randel, 1993; Shuckburgh and Haynes, 2003).

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-5 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.
2 Data and Methodology 10

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-hourly) and better horizontal (31 km versus 80 km) and vertical (60 levels versus 137 levels) resolutions and extends higher 15 into the middle atmosphere (0.1 hPa versus 0.01 hPa or 65 km versus 80 km). The ERA5 reanalysis has replaced the ERA-Interim reanalysis from 31 August 2019. Due to stratospheric temperature biases for the 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 above 40 km persists for much of the period from 1979 (Hoffmann et al., 2019). In addition to the higher spatial and temporal resolution, other key improvements of ERA5 compared 20 to ERA-Interim are a better ability to resolve synoptic-scale features like hurricanes and tropical cyclones as well as a better representation of the tropospheric circulation. Moreover, data from many recent satellite instruments are now additionally assimilated (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. 25 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 on their original model levels interpolated to log-pressure levels for the residual circulation calculations. The difference between ERA-Interim and ERA5 is that the assimilation cycles start at 06:00 UT and 18:00 UTC for ERA5 and at 00:00 UTC and 12:00 UTC for ERA-Interim. In addition, the dynamical fields in ERA5 are archived hourly, and no longer as accumulation over every 3 hours as in ERA-Interim. Note that both re- 30 analyses are part of the SPARC Reanalysis Intercomparison Project (S-RIP) (Fujiwara et al., 2017), which is a coordinated inter-comparison of modern global atmospheric reanalyses. In particular, the Chapter 5 of the S-RIP report evaluates several commonly used metrics of the BDC calculated from reanalysis fields. This work contributes to this assessment of the BDC metrics in the ERA5 reanalysis.

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 5 given by Andrews et al. (1987) ∂ u ∂t where a is the Earth's radius, φ is latitude. z = −H · ln( p p s ) is 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 10 respectively the zonal mean, meridional and vertical wind velocity. ℑ corresponds to the total zonal momentum forcing. The Coriolis frequency is f = 2 · Ω · sin φ and Ω is the Earth's rotation rate. Here, and in the following, the zonal mean and the deviations from the zonal mean are respectively indicated by an over-bar and a prime. The partial derivatives with respect to z (vertical) and φ (meridional) direction are indicated by ∂ ∂ z (.) and ∂ ∂ φ (.) respectively. As the advective BDC is a Lagrangian mean circulation, the Eulerian mean velocity is not a sufficient diagnostic metric. 15 However, a useful proxy for the Lagrangian mean circulation under time-averaged conditions is provided by the residual mean meridional circulation (v * , w * ) in the TEM framework McIntyre, 1976, 1978;Dunkerton et al., 1978;Holton, 1990). The latitudinal and vertical components of the residual mean meridional circulation (v * , w * ) and the mean mass streamfunction 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, 1978) where X urGW D is the residual mean nonconservative forcing unresolved by the model and DF is the normalised EP-flux 25 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 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 the 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 10 and gravity waves. DF can be decomposed into wave numbers ranging from 1 to 180 following previous studies (Ern et al., 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 15 equal to the unresolved wave drag by the ECMWF model grid, therefore, contributing to the gravity wave forcings in the 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 20 waves, then, to maintain the steady conditions, the BDC has to stay non-zero McIntyre, 1976, 1978;Eliassen 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). Using the downward control principle (e.g. Haynes et al., 1991;Rosenlof and Holton, 1993;Garcia and Boville, 1994), it can be written as 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 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 5 circulation transit time (RCTT) using a two dimensional backward trajectory model driven by the residual mean meridional circulation (v * , w * ). Used as 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 estimates see Ploeger et al. (2015b).

10
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 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 15 residual (ε). For a given metric of the advective BDC strength, BDC metric (herein U, T, w * , ψ * and wave drags), at any given latitude (φ ) and altitude (z) position in the stratosphere, the regression model can be written simply as following The estimated QBO and ENSO coefficients with the regression fit as a function of latitude and altitude are normalized by the standard deviation (SD) of the QBO and ENSO predictors, which are the QBO index at 50 hPa and the Multivariate 20 ENSO Index (Wolter and Timlin, 2011). We named these normalized coefficients as the QBO and ENSO amplitude variations.
Because of the presence of lags (τ qbo and τ enso ,) in the QBO and ENSO terms, the problem is nonlinear and the residual may have multiple minima as a function of the parameters. In order to determine the optimal values of the lags the residual is first minimized at fixed lag and then selected from a range of possible lags. 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 25 solar cycle are intentionally omitted for simplifying the description as these terms are not diagnosed in this study. For more details, please see our previous studies (Diallo et al., 2012(Diallo et al., , 2017(Diallo et al., , 2018Tao et al., 2019).

Inter-annual variability of the upwelling
Evaluating the large-scale ascent and descent of air masses 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 mass flux from the ERA5 reanalysis to the ERA-Interim reanalysis using the residual vertical velocity (w * ). Figure 1a-i show 5 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 and statistically significant differences at 95% confidence interval in the strength of the tropical upwelling and extratropical downwelling arise in three distinct regions of the stratosphere (tropical pipe, mid-latitude surf zone 10 and polar regions) (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 statistically significant weaker annual mean tropical upwelling in the UTLS, consistent with the well-known too fast ERA-Interim tropical diabatic upwelling compared to observations (Dee et al., 2011;Seviour et al., 2011;Ploeger et al., 2012). Using trace gas reconstruction and decomposition with a Lagrangian model driven by ERA-Interim, Ploeger et al. (2012) found tropical diabatic upwelling to be about 40% too fast in the tropical tropopause layer (Fueglistaler et al., 2009). The reduced tropical upwelling in ERA5 vertical residual circulation velocities, as compared to ERA-Interim, is also consistent with the diabatic ERA5 heating rates, which are 30-40% weaker in ERA5, therefore, correcting the known ERA-Interim bias (Ploeger et al., 2021). 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 significant weaker downwelling at 95% confidence interval. Conversely in the polar vortex region, the ERA5 reanalysis shows a 20 stronger large-scale downwelling of the air mass than the ERA-Interim, suggesting potential differences in the polar vortex strength, but further studies are needed to prove that speculation. 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 a significant seasonal variation 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, 25 these large and significant differences in the w * are displaced toward the winter hemisphere in the upper stratosphere with a maximum of about 1.5 mm · s −1 . The reasons of these discrepancies in tropical upwelling and extratropical downwelling are likely associated with the differences in wave forcings resulting from the improvements in the ERA5, including higher vertical and horizontal resolutions, higher model top, and new gravity wave parameterization scheme (Orr et al., 2010;Hersbach et al., 2020). 30 To quantify the differences in the circulation circulation differences between the two reanalyses, we average the annual mean w * into vertical profiles based on the turnaround latitudes of tropical upwelling and extratropical downwelling together with the associated differences and their statistical significance at 2σ level for the 1979-2018 period (Fig. 2 a-f). The vertical profiles of the w * The w * vertical profiles show a good agreement in the w * structure between the two reanalyses for the 1979-2018 time the ERA-Interim reanalysis (e) and its associated differences (f). The panels (g-i) show the JJA mean from the ERA5 reanalysis (g), the ERA-Interim reanalysis (h) and its associated differences (i). Grey line indicates the zero w * contours. Grey dots in panels (c, f, i) indicate regions where the differences between the ERA5 and ERA-Interim reanalyses are statistically significant at 95% estimated using student t-test. The climatological tropopause from the ERA-Interim is indicated as the black dashed horizontal line.  Differences between the ERA5 and ERA-Interim reanalyses are significant where deviations from zero exceed the 2σ range. period ( Fig. 2a-c). Despite the similarities in the vertical structure, the relative differences in the w * the relative w * differences show a significantly stronger large-scale downwelling in ERA5 than in ERA-Interim with a hemispheric asymmetry ( Fig. 2 d-f). In the southern hemisphere, the differences in the large-scale downwelling the large-scale downwelling differences are about 20 % stronger in ERA5 than in ERA-Interim below 30 km, and are even larger above with a maximum of about 60 % at 35 km (Fig. 2 c, f). These differences are statistically significant at 2σ level as indicated by the errorbars. In the tropics, the ERA-Interim reanalyses are statistically significant at 95% confident level estimated using student t-test. 35 km, ERA5 reveals a statistically significant and stronger tropical upwelling than ERA-Interim. In the northern hemisphere, the downwelling differences between the two reanalyses are negligible below 30 km and reach 30 % above (Fig. 2c, f). The errorbars, which correspond to the statistically significant area at 2σ level, indicate that the large differences in upwelling and downwelling correspond to the regions where the variability is large. 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-5 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).
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 intercomparisons of the upwelling strength of the air mass entering the stratosphere (e.g., Butchart et al., 2010). The tropical 10 upwelling and extratropical downwelling patterns also agree 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-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 significant seasonal variation of the differences in the w * w * differences between the two reanalyses is also shown in Figure 3c. The significantly weaker tropical upwelling in ERA5 varies seasonally between 0 • and 20 • N with 15 a maximum of about 0.2 mm · s −1 . In the mid-latitude, the large-scale downwelling is also significantly weaker and varies seasonally with a maximum occurring during the boreal summer (June-July-August). In the polar vortex region, ERA5 shows a significant large-scale downwelling, which is stronger than in ERA-Interim, as already evident from Fig. 1-2. The mass flux differences maximize in the polar vortex region of the southern hemisphere. These significant seasonal variations in the tropical upwelling and the extratropical downwelling result 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 between the two reanalyses for the 1979-2018 time period. Figure 4a-i show a remarkably good agreement between the two 5 reanalyses in the morphology of the advective BDC, regarding an ascent of air mass in the tropics, a motion in the stratosphere toward the mid-and high altitudes and latitudes, and descent into the mid-and high latitude regions. Note that clear differences exist in the annual mean upwelling and downwelling circulation cells of the mass stream-function in the tropical pipe, midlatitude surf zone, and polar vortex regions. Similarly to the w * , also the residual circulation mass stream-function shows weaker upwelling in the tropical pipe and the mid-latitude surf zone than in ERA-Interim, but stronger in the polar regions 10 ( Fig. 4a-c). Similarly to the w * , the residual circulation mass stream-function upwelling in the tropical pipe and the midlatitude surf zone is also weaker in ERA5 than in ERA-Interim, but stronger in the polar regions ( Fig. 4a-c). The maximum and significant 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 15 stream-function between the two reanalyses are negative between 0 and 20 • N, which extends from the UTLS into the upper stratosphere. This negative difference indicates a weaker circulation in ERA5 and is consistent with the weaker upwelling found in Fig. 3. However, the differences in the mass stream-function between the two reanalyses are significantly positive between 50 • 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 significantly slower in ERA5 than in ERA-Interim. 20 During the boreal summer (June-July-August), these discrepancies in the residual circulation are reversed ( Fig. 4g-i)). These dominating seasonal features of the differences in the residual circulation suggest that the significant improvements in ERA5 likely induce a southward shift of upwelling cells in the UTLS during winter and northward shift during summer compared to ERA-Interim.
For further insights into the circulation differences between the two reanalyses, we also evaluate the annual and seasonal 25 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 and calculated using v * and w * in the equation (2) (e.g., Birner and Bönisch, 2011;Garny et al., 2014;Ploeger et al., 2015a, b). Note 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 reanalyses, e.g. the shorter transit time in the tropics (faster BDC) and longer transit time in the extra-30 tropics (slower BDC) for the 2010-2018 period (Fig. 5a, b). However, the RCTT also shows clear and significant differences in the residual circulation, consistent with the discrepancies in residual mass stream-function and residual vertical velocity previously discussed. The annual mean RCTT in ERA5 exhibits longer transit time (below about 0.5 years) associated with the significantly slower tropical ascent of the transition and shallow branch of the BDC (Lin and Fu, 2013; indicate regions where the differences between the ERA5 and ERA-Interim reanalyses are statistically significant at 95% estimated using student t-test. The climatological tropopause from the ERA-Interim is indicated as the black dashed horizontal line. slow integrated residual circulation in the ERA5 reanalysis is consistent with the differences in the w * w * differences as well as with the diabatic RCTT in ERA5 (Ploeger et al., 2021). In the northern hemisphere, the RCTT from the ERA5 reanalysis shows a longer transit time (differences below about 0.5 years). In the southern hemisphere, the residence time in ERA5 tends to be regions where the differences between the ERA5 and ERA-Interim reanalyses are statistically significant at 95% estimated using student t-test. The climatological tropopause from the ERA-Interim is indicated as the black dashed horizontal line.
longer than in ERA-Interim. The analysis of the seasonal variations also shows significant patterns of differences consistent with those shown in the annual mean RCTT (Fig. 5d-i). The differences vary seasonally and the maximum difference is found in the polar northern 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 modes of variability in the ascending branch of the BDC on seasonal to intrannual timescales is the QBO (Lindzen and Holton, 1968;Plumb and Bell, 1982). Composed of alternating westerly and easterly zonal wind shears, the QBO propagates downward from the tropical 5 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 phase (Osprey et al., 2016;Newman et al., 2016), is clearly visible in both reanalyses. Apparent differences are also observed 10 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 reanalysis exhibits stronger QBO westerly and easterly phases compared to the ERA-Interim reanalysis in the tropical stratosphere above 20 km (Fig. 6c). However, the QBO phases in the ERA5 reanalysis are weaker than in the ERA-Interim reanalysis between the 15 and 20 km. According to previous findings, the QBO westerly phase is sensitive to many model details, including the parametrized non-orographic gravity wave 15 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 enhanced non-orographic gravity waves. In the ERA5, the use of a new non-orographic gravity wave drag parameterization, which is different from the ERA-Interim, is likely the causes of the observed differences Therefore, the use of a new non-orographic gravity wave drag parameterization in ERA5 likely is the cause of the observed differences (Orr et al., 2010;Scaife et al., 20 2002; Lott et al., 2012;Richter et al., 2014). In ERA-Interim, the Rayleigh drag used to be a substitute for the non-orographic gravity wave drag. In ERA-Interim, Rayleigh drag was applied as a substitute for the non-orographic gravity wave drag. When a Warner and McIntyre (2001) type scheme, a non-orographic spectral gravity wave parameterization model, was introduced then this Rayleigh drag could be switched off. For ERA5, a Warner and McIntyre (2001) type non-orographic spectral gravity wave scheme was introduced as parameterization model and hence the Rayleigh drag could be switched off. In addition, the 25 finer vertical and horizontal resolution in ERA5 likely leads to better representation of the synoptic and small-scale waves as sources of gravity waves are better represented, such as convection, which is important for gravity wave driving of the QBO (Lindzen and Fox-Rabinovitz, 1989).  Both reanalyses agree very well in the phase and periodicity of the westerlies and easterlies. In addition to gravity waves, the QBO 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., a) b) c) 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 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). Westerly shear reduces the tropical upwelling 5 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;Randel and Wu, 1996;Randel et al., 1999;Randel, 1987;Hamilton, 1998). The magnitude 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 10 temperatures between 20 and 30 km, leading to enhanced tropical upwelling by the secondary meridional circulation associated with the QBO (Randel et al., 1999;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  Diallo et al., 2018;Ray et al., 2020). The magnitude 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 than in ERA-Interim likely related to the Holton-Tan effect (i.e. 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 5 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 the secondary meridional circulation, we analyzed the QBO-and ENSO-induced variability in the stratospheric residual circulation from variations in w * and ψ * . Considered as one of the major modes of BDC modulation , the coupled atmosphere-ocean phenomenon, ENSO, is defined as extreme sea surface temperature (SSTs) Zonal mean climatology of w * is overplotted as dashed grey lines.
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, 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). 5 To quantify the QBO-and ENSO-induced variability in w * and ψ * , the statistical analysis (eq. 10) is performed by explicitly including ENSO and QBO predictors to isolate their modulations of the BDC. The QBO and ENSO coefficients in the regression fit are normalized by the standard deviation of the predictors' proxies for the 1979-2018 period, except for the temperature and zonal mean wind figures, and will be called here the QBO and ENSO amplitude variations. This direct approach gives similar results as the differentiating approach of the residuals (Diallo et al., 2017(Diallo et al., , 2018. The results of the QBO-induced variability in 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 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 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 5 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 Our analysis of the QBO-induced variability in the residual stream-function circulation is also consistent with the QBO impact on w * (Fig. 8a, b). ERA5 exhibits weaker anomalies than ERA-Interim. This specific difference in the UTLS 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 5 30 km where ERA-Interim shows larger negative anomalies than the ERA5. In the northern hemisphere upper stratosphere, ERA5 also reveals larger negative w * anomalies than the ERA-Interim, which is likely due to the differences in wave activity (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 agree-ENSO strengthens the residual circulation of the BDC, particularly the shallow branch (Fig. 9c, d) Yang et al., 2014). The ENSO-induced changes in the deep branch are less evident compared to the shallow and transition branches, nevertheless they are stronger in ERA5 than in ERA-Interim. These 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 trop-5 ical 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 less evident ENSO-induced changes in the deep BDC branch . As all differences in the circulation and in the natural variability disclose possible discrepancies in the BDC forcings, we analyze the contribution of the planetary 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 reanalyses, we evaluate the net forcing, the planetary and the gravity wave driving of the BDC (Haynes et al., 1991;Rosenlof and Holton, 1993;Newman and Nash, 2000;Plumb, 2002;Shepherd, 2007). Figures 10a-i show the 15 annual mean climatology of the net forcings of the BDC and the contributions due to the planetary and gravity wave drag together with the associated differences between the ERA5 and ERA-Interim reanalyses for the 1979-2018 time period. Note that the zonal wave forcing ℑ can induce meridional wind contributions δ v * (Holton, 1986;Holton et al., 1995). Assuming in the zonal momentum balance (eq. 4) steady state 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 20 to ℑ = − f * v * (Holton, 1986;Holton et al., 1995). 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 indicating eastward forcing would weakens the BDC (δ v * is equatorward), and negative wave drag indicating westward forcing would weakens 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). 25 Overall, the climatological structure of the mean wave drag agrees well between the ERA5 and ERA-Interim reanalyses, including the positive wave drag in the easterly zonal mean wind regions and negative wave drag in the westerly zonal mean wind regions (Fig. 10a-f). However, significant differences occur in the amplitude of the wave drag ( Fig. 10g-i). The net forcings of the BDC shows that the gravity wave drag contributes the most to the significant lower stratospheric circulation differences between the ERA5 and ERA-Interim. The breaking of synoptic and small-scale waves near the tropical tropopause 30 layer and the equatorward 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;. In the UTLS below 20 km, the previously reported weaker tropical upwelling in ERA5 compared to ERA-Interim is due to a significantly weaker gravity wave breaking at the equatorial-ward upper flank of the subtropical jets in the ERA5 (Fig. 10g, i). Thus, this difference in the strength of gravity wave breaking in the lower stratosphere UTLS explains the weaker ERA5 tropical upwelling as observed in the residual vertical velocity, residual mass stream-function and the RCTTs. The significant contribution of planetary wave is mainly located either below the tropopause or far from the tropical upwelling driving regions, therefore, has a minor contribute to the differences in shallow branch of the BDC between the ERA5 and ERA-Interim reanalyses. In the extratropical upper stratosphere above 30 km, both planetary and gravity wave contribute to the circulation differences between the ERA5 and 5 ERA-Interim reanalyses with a stronger wave forcings in ERA5 than in ERA-Interim (see blue area in 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 significantly weaker gravity wave breaking at the equatorial flank of the subtropical jets in the ERA5 (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 layer and the upper flank of the subtropical jets are the primary forcing sources of the transition and 10 shallow branches of the BDC (Plumb, 2002;Plumb, 2011;Plumb, 2019). 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 governed by differences in the contribution of both the planetary and gravity wave breaking ( Fig. 10g-i).
While the planetary wave drag is stronger in the extra-tropics between 30 and 40 km in ERA5 than in ERA-Interim, it is weaker 15 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 the non-orographic gravity wave parametrization from Warner and McIntyre (2001). This means that ERA5 does not use 20 this simple form of Rayleigh drag 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 significant differences in wave drag between the ERA5 and ERA-Interim reanalyses reveal seasonal vari- . DJF (a, c) and JJA (b, d) seasonal mean variations of the mass stream-function function (ψ * in kg · m −1 · s −1 ) from the ERA5 (a, b) and the associated differences between ERA5 and ERA-Interim (c, d) for the 1979-2018 time period. The ψ * is calculated from the w * , resolved wave drag (DF), unresolved wave drag (X urGW ) and wind tendency using downward control (DC) principle (Haynes et al., 1991).
the southern hemisphere during the boreal winter, which is, in turn, consistent with the reported stronger large-scale downwelling in polar regions (Fig. A1h, i). During boreal summer, the differences in the tropical lower and extratropical upper stratosphere induced by mainly planetary and gravity wave dissipation are also consistent with the discrepancies in the BDC between the ERA5 and ERA-Interim reanalyses. In the lower stratosphere, the gravity wave breaking at the equatorial flank of the subtropical jet causes the significant differences between the ERA5 and ERA-Interim tropical upwelling (Fig. A2a-i). The 5 observed maximum of the seasonal variations in the tropical upwelling at 70 hPa between 0 and 20 • N is due to a significantly weaker gravity wave breaking in the ERA5 lower stratosphere 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, consistent 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 To further link the differences in the upwelling and wave drag, we calculate the DJF and JJA mean stream-function at 70 hPa (∼18.5 km) using the downward control principle (Haynes et al., 1991), described in equation (8). In addition to being commonly used as the reference level for model inter-comparisons of the upwelling strength , the 70 hPa pressure level works best for RCTT and age of air in inter-model correlations for almost all the stratosphere (Dietmüller et al., 2017). Figure 11 shows the DJF and JJA mass stream-function seasonal variation at 70 hPa calculated from the ERA5 5 and ERA-Interim w * , resolved wave drag (DF) and unresolved wave drag (X urGW ) for the ERA5 reanalysis together with the differences between ERA5 and ERA-Interim for the 1979-2018 time period. Clearly, the gravity wave drag is the main driver causing the weaker tropical upwelling in ERA5 compared to ERA-Interim with largest effects occurring within the tropics ( Fig. 11 c, d). Due to small differences between ERA5 and ERA-Interim in the contribution of the resolved wave drag (DF), including the resolved planetary and gravity wave drag, we can also conclude further that the upwelling differences between the 10 two reanalyses in term of the gravity wave drag results from the unresolved gravity wave contribution, i.e. the parameterized and imbalance gravity wave drag (Alexander and Rosenlof, 1996;McLandress et al., 2012;Ern et al., 2014). This suggests that the improvement in parameterized non-orographic gravity wave drag, which impacts tropical regions where sub-grid-scale gravity waves from convective sources are dominant, may certainly be the reason of these differences (Chun et al., 2004), consistent with the upwelling and wave drag differences. According to Butchart et al. (2010), the state-of-the-art chemistry-15 climate models show a good representation of the resolved wave contribution to driving the annual mean tropical upwelling at 70 hPa with a contribution of 70.7%from resolved waves, 21.1% from orographic gravity wave drag (parametrized in ERA5 and ERA-Interim) and 7.1% from non-orographic gravity waved drag (not included in ERA-Interim but include in ERA5). With a contribution of about 30% according to climate model mean, the differences in the unresolved gravity wave drag between the two reanalyses is consistent with the contribution range to the upwelling. Nevertheless, one should keep in mind that the 20 representation of gravity waves in reanalyses and climate models remains not well known (Seviour et al., 2011;Shepherd, 2014). The uncertain aspect of the reanalyses related to the representation of unresolved (sub-grid scale) processes such as gravity-wave drag, convection, and boundary-layer physics have been improved compared to the ERA-Interim (Hersbach et al., 2020).
4 Advective BDC trends 25 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. 12a-d). The ERA5 reanalysis shows a long-term acceleration of the BDC, consistent with all reanalysis data sets, except the CFSR reanalysis (Diallo et al., 2012;Abalos et al., 2015;Miyazaki et al., 2016;Chabrillat et al., 2018;Ploeger et al., 2019), and with long-term climate model simulations projections (e.g., Butchart et al., 2010;Hardiman et al., 2014). The annual mean trend of the residual circulation 30 mass stream-function shows a positive trend in the northern hemisphere and weakly negative trend in parts of the southern hemisphere, indicating a not significant strengthening of the shallow and deep branches of the BDC in both hemispheres (Fig. 12a). The negative BDC trend in the inner tropical UTLS suggests a weakening of the transition branch probably linked to the tropopause rise induced by the combined effect of tropospheric warming by greenhouse gases and stratospheric cooling by ozone depletion (Randel et al., 2000;Santer et al., 2003;Seidel and Randel, 2006;Son et al., 2009;Vallis et al., 2015;Oberländer-Hayn et al., 2016;Šácha et al., 2019;Eichinger and Sacha, 2020). The decomposition of the wave drag trend into planetary and gravity wave forcing provides an estimate of the impact of these two wave forcings on the structural changes in the transition, shallow and deep branches of the BDC. The BDC trend in ERA5 is induced by the combined contribution of 5 the gravity and planetary wave breaking (Fig. 12c, d). The negative trends in the net wave forcings of the BDC between 35 and 45 km indicate an intensified planetary and gravity wave breaking in both hemispheres (Fig. 12b-c). While the shallow branch is accelerated mainly by enhanced gravity wave breaking in the tropical and subtropical UTLS region lower stratosphere, the deep branch is driven by a sum of the contribution from the planetary and gravity wave breaking at high altitudes. The planetary wave breaking seems to contribute the most to the acceleration of the deep branch of the BDC, consistent with previous studies 10 (Sigmond and Shepherd, 2014;Abalos et al., 2015;Šácha et al., 2019). 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. 12b). Also note that the southern hemispheric circulation cell exhibits an apparent structural difference compared to the northern hemispheric circulation cell. A detailed analysis of the BDC trends is presented in further studies using the stratospheric age of air and its spectrum (Ploeger et al., 2021). Note that the ERA-Interim still shows a 15 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 shows a comparable inter-annual variability between the ERA5 and ERA-Interim reanalyses at the pressure level of 70 hPa (Fig. 13). The estimated annual mean upwelling mass flux is about 6.17 × 10 9 kg · s −1 for ERA5 with a standard 20 deviation of 0.71 and about 6.86 × 10 9 kg · s −1 for ERA-Interim with a standard deviation of 0.99 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 * for the 1979-2018 time period is about 1.5 % per decade, consistent with the observed BDC changes of about 1.7 % per decade (Linz et al., 2017;Fu et al., 2019) and climate model predictions of an increase in the strength of the BDC of about 2 % per decade due to the 25 increasing GHGs in the atmosphere Garny et al., 2011;Hardiman et al., 2014;Eichinger and Sacha, 2020). While the climate model projected trend is statistically significant, the mass flux trend in ERA5 is not statistically significant at 95% using student t-test with two tail distribution. 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 * is not consistent with the ERA5 trend.
Therefore, the ERA5 w * trend is similar to the trend from climate models, which implies that the ERA5 w * estimated with 30 the standard TEM formula can be used as a proxy for climate model validation, which was not the case for the ERA-Interim because of a questionable w * (Seviour et al., 2011). This negative long-term trend is only present in the ERA-Interim reanalysis when using the standard TEM formula to calculate the w * and is 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± 0.053 × 10 9 kg · s −1 · dec −1 with a sigma of 0.053 instead of -0.47± 0.53 × 10 9 kg · s −1 · dec −1 for the ERA-Interim with a sigma of 0.53. 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 climatology, 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, in particular the QBO and the ENSO, on the advective BDC as well as the impact of planetary and gravity wave drag on the advective BDC and its trend. 5 In our comparisons of the circulation from the ERA5 reanalysis with the ERA-Interim reanalysis, 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 QBO-and ENSO-induced impact. Despite the good agreement in the spatial structure, there are significant differences between the ERA5 and the ERA-Interim reanalyses in the strength of the advective BDC and in its modulation by the natural variability in the UTLS and upper stratosphere regions. The slower advective BDC and its strong modulations by the 10 natural variability in the ERA5 reanalysis is due to weaker planetary and gravity wave drags. In the tropical pipe region below 20 km, the tropical upwelling is up to 40 % weaker in ERA5 than in ERA-Interim mainly due to a significantly weaker gravity wave breaking at the equatorial flank of the subtropical jet. The differences in planetary wave drag are significant but weaker than the differences in the gravity wave drag at key regions of upwelling forcings, therefore, they have minor contribution to the slower upwelling differences. In the extra-tropics, the large-scale downwelling near the polar vortex is stronger in ERA5 than 15 in ERA-Interim linked to significant differences in planetary and gravity wave drag, which might impact on the representation and strength of the polar vortex. 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, the differences between the ERA5 and ERA-Interim reanalyses show a hemispheric asymmetry and a seasonal variation with a maximum difference of about 60 % in downelling (stronger in ERA5) 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 5 upwelling in both reanalyses. Conversely, the QBO easterly shear modulation induced anomalously cold tropical temperatures, leading to enhanced tropical upwelling by the QBO secondary circulation. Similarly to the QBO-induced secondary meridional circulation in the temperatures, the regression analysis of the QBO-induced structural changes in w * and ψ * shows a very good agreement in the modulations of the BDC between the two reanalyses. 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 easterly and westerly 10 shear modulations. Besides the good structural agreement, noticeable differences are also found in the strength of the QBO easterly and westerly shear modulations. Between 15 and 20 km, the QBO phases in ERA-Interim are stronger than in ERA5, while above 20 km, the QBO modulation is stronger in ERA5 than in ERA-Interim.
The regression analysis of the ENSO-induced structural changes in w * and ψ * shows a very good agreement in the modulations of the BDC between the two reanalyses. However, significant differences have also been found in the strength of the 15 advective BDC modulations by ENSO. Analogously to the QBO effects, 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 differences in the advective BDC in the UTLS region between the ERA5 and ERA-Interim reanalyses, including the tropical upwelling and mid-latitude downwelling, 20 are mainly due to a weaker gravity wave breaking at the equatorial flank of the subtropical jets in ERA5. The differences in the planetary wave drag between the ERA5 and ERA-Interim reanalyses are significant, but weaker than the differences in gravity wave drag in the lower stratosphere, therefore, not dominant. 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, in particular, therefore in the 25 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. In addition, progress has also been 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 (Hersbach et al., 2020).
Finally, the estimates of the advective BDC trend in ERA5 show a global acceleration of the annual mean residual circulation 30 stream-function of about 1.5 % per decade for the 1979-2018 period. Althrough Even not statistically significant at 95%, this acceleration rate is consistent with observed acceleration of the BDC of about 1.7 % per decade (Thompson and Solomon, 2005;Linz et al., 2017;Fu et al., 2019) and with future climate model predicted rate of the strength of the BDC of about 2 % per decade due to the increasing level of atmospheric GHGs Garny et al., 2011;Hardiman et al., 2014). This implies a recommended use of ERA5 w * as proxy for climate model validations instead of the ERA-Interim. The 35 strengthening of the BDC is induced by the increasing planetary and gravity wave breaking in both hemispheres above 30 km.
The trends in the residual circulation stream-function are stronger in the northern hemisphere than in the southern hemisphere, consistent with the trends of the net wave forcing of the BDC. This hemispheric asymmetry in the wave breaking trends is likely 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). 5 Our findings suggest that the advective BDC from the kinematic ERA5 reanalysis is well suited for climate model validation in the UTLS and mid-stratosphere when using the standard formula of zonally-averaged zonal momentum equation. The reported differences between the two reanalyses may also affect the nudged climate model simulations. According to recent findings, nudged climate models show slightly stronger upwelling compared to the free-running simulations and exhibit marked differences compared to the direct estimates of the upwelling from the reanalysis (Chrysanthou et al., 2019;Davis et al., 10 2020). Therefore, additional studies are needed to investigate whether or not nudging climate models toward ERA5 reanalysis will accurately reproduce the tropical upwelling trends comparable to the free-running versions and to directly estimated residual circulation from the ERA5 reanalysis. Finally, studies seeking for a better understanding of the impact of the new non-orographic gravity wave parameterization scheme, higher model top, and the representation of the sponge layer in ERA5 on the differences in the upper stratosphere and polar regions are also needed. Data availability. The ERA5 reanalysis (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