The Brewer-Dobson circulation in CMIP6

The Brewer-Dobson circulation (BDC) is a key feature of the stratosphere that models need to accurately represent in order to improve the representation of surface climate variability. For the first time, the Climate Model Intercomparison Project includes in its phase 6 (CMIP6) a set of diagnostics that allow for careful evaluation of the BDC. Here, the BDC is evaluated against observations and reanalyses using historical simulations. CMIP6 results confirm the well-known inconsistency in BDC trends between observations and models in the middle and upper stratosphere. The increasing CO2 simulations feature a robust 5 acceleration of the BDC but also reveal large uncertainties in the deep branch trends. The very close connection between the shallow branch and surface temperature is highlighted, which is absent in the deep branch. The trends in mean age of air are shown to be more robust throughout the stratosphere than those in the residual circulation. The paper reflects the current knowledge and main uncertainties regarding the BDC.

. List of CMIP6 models with TEM and/or AoA diagnostics used in this study. AoA is only provided in the historical runs for MRI.

Model
Model

Data and Methods
CMIP6 models providing the necessary TEM output, as described in Gerber and Manzini (2016), have been used. We refer 50 to that paper for the specific diagnostic description. In addition to the TEM diagnostics, some models additionally provided mean age of air (AoA). Table 1 shows the models used, their number of levels, the model top and the corresponding available variables. While each model has a different horizontal resolution (not shown), they are all in the range between 1 and 2 degrees in longitude and latitude. Note that two more models output TEM variables, CESM2 and Can-ESM5, but we did not include them in the analyses because they have low model tops (2.25 hPa and 1 hPa, respectively), and did not represent the residual 55 circulation structure adequately (not shown).
Model results have been compared to reanalysis over the historical period. Because there is a large spread in the BDC values obtained from reanalyses , we used three reanalyses: ERA-Interim (Interim European Centre for Medium-Range Weather forecasts Reanalysis, Dee et al. (2011)), JRA-55 (Japanese 55 year Reanalysis, Kobayashi et al. (2015)), and 60 MERRA (Modern-Era Retrospective Analysis for Research and Applications, Rienecker et al. (2011)). We have used the residual circulation for these three reanalyses from Abalos et al. (2015), which covers the period 1979-2012. In order to compare AoA with observational estimates, we used AoA derived from the Michelson Interferometer MIPAS (Stiller et al., 2012;Haenel et al., 2015). Here, we use a new version of AoA derived from an updated retrieval of SF 6 (Stiller et al., 2020). Furthermore, AoA derived from GOZCARDS N 2 O data is used (Linz et al., 2016), as well as AoA data derived from in-situ gases and other external forcings, and will be used to examine past climatology and trends and compare to observations or reanalysis when possible. For comparison purposes we have focused on the period 1975-2014. Note that this period encompassess the reanalysis period considered , but it is slightly longer. This is done to enhance statistical significance trend calculations for model output. The small difference in the period considered is assumed to have a negligible impact on the climatological BDC. The 1pctCO2 simulations are initialized from preindustrial (1850) conditions and are 150 years long, with 75 CO 2 concentrations increasing gradually at a 1%/year rate (Eyring et al., 2016). In order to establish statistical significance we have used a two-tailed Student t test at the 95% confidence level. We consider there is intermodel agreement if 66% of the models are consistent.

Representation of the BDC in CMIP6 historical simulations
This section aims to assess the degree of agreement in the BDC climatology and trends between CMIP6 historical simulations 80 and observations or reanalysis data. previous model intercomparison studies such as CMIP5 .

Climatology and seasonality
In order to examine the quantitative differences in more detail, the tropical upwelling mass flux is examined. This is computed as the net upwelling between the annual mean turnaround latitudes (i.e., the latitudes separating the upwelling and downwelling regions). The calculation is based on the streamfunction, which in turn is computed from the meridional component of the 90 residual circulation provided as model output,v * . The streamfunction is obtained as where p is pressure, φ is latitude, g the gravitational constant on Earth, and it is assumed thatv * tends to zero as p → 0. The upwelling mass flux is then computed at each level as 95 whereΨ * max andΨ * min are the maximum and minimum values of the residual streamfunction at each pressure level, which correspond to the northern and southern turnaround latitudes, respectively (Rosenlof, 1995).   Figure 1. Annual mean climatology for the multi reanalysis mean (a) and multi model mean for the historical simulations (b) of the vertical component of the residual circulation (w * , in mm · s −1 , shading) and residual streamfunction (Ψ * , in kg · m −1 · s −1 , contours). Black thick contours indicate the location of the turnaround latitudes. The red contour in the right panel shows the turnaround latitudes for reanalyses.
Black dots represent regions where there is disagreement in the sign ofw * for more than 66% (2/3) of the individual reanalyses or models (which happens only around turnaround latitudes). Figure 2 shows the seasonality in the tropical upwelling mass flux for the lower (70 hPa) and upper (1.5 hPa) stratosphere, representative of the shallow and deep branches, respectively. Note that, while the level of 70 hPa is commonly used to represent the shallow branch, 1.5 hPa is higher than usually considered for the deep branch. We argue that this level is optimal for the 100 characterization of the deep branch, since tropical upwelling maximizes at this level in the upper stratosphere ( Fig. 1). All models show a generally consistent seasonality, with an annual cycle peaking in November-December in the shallow branch, and an amplitude of about 50% of the climatological mean, and a semi-annual cycle peaking in June and December for the deep branch, with an amplitude of about 80%. The seasonality is consistent with that of reanalyses, and the intermodel spread is of similar magnitude to the reanalysis spread. In particular, the intermodel spread is over 40% of the climatological mean 105 for the lower stratosphere and over 30% for the upper stratosphere. The annual cycle in the lower stratosphere has been linked to seasonality of wave forcing in the extratropics, subtropics and tropics (e.g. Randel et al., 2008;Ueyama et al., 2013;Ortland and Alexander, 2014;Kim et al., 2016). The semi-annual cycle in the upper stratosphere has been less studied, and is probably linked to the combined annual cycles of the downwelling in each hemisphere and to the secondary circulation associated with the Semi-Annual Oscillation (e.g. Garcia et al., 1997;Young et al., 2011).

110
As mentioned in the Introduction, the mean age of air provides an estimate of the net transport circulation strength that can be compared to observational estimates. Figure 3 shows the AoA climatology at 50 hPa for the models that provide this quantity (see Table 1), together with the observational estimates described in Section 2. The simulated AoA values show considerable spread across models, as previously shown for Chemistry-Climate Model Intercomparison project (CCMI) simulations (e.g., Dietmüller et al., 2017). The global mean age values vary by a factor of 2, between 2.5 and 5 years approximately. Nevertheless, J a n AoA (e.g. Linz et al., 2016).

Past trends
In this section we examine the BDC trends over the historical period, in particular over the last four decades. Figure 4 shows the (SH) than in the Northern Hemisphere (NH). Consistently, the residual circulation accelerates throughout the stratosphere, with enhanced tropical upwelling and polar downwelling, strongest in the SH. Note that there is also reduced downwelling in midlatitudes in both hemispheres. The larger BDC polar downwelling trends in the SH are consistent with recent results using CCMI models, and reflect the contribution of ozone depletion in the Antarctic lower stratosphere to the BDC trends (Polvani et al., 2018Abalos et al., 2019). In order to better capture this signal, the trends are shown separately for the AoA 50 hPa 1991-1997 from the timeseries in order to avoid the influence of Pinatubo volcanic eruption on the trends. We caution that the 135 trends are computed over short periods and the residual circulation presents high interannual variability, such that there is no statistical significance of the trends (panels d-f). Nevertheless, the influence of the ozone hole is clearly seen in the different trend magnitudes.
While Figure 4 shows the general trend behavior for the multi-model mean, it is important to assess the robustness of the trends across different members and models, especially given the relatively short periods under consideration. Figure 5 140 examines the inter-model spread as well as the inter-member spread for each model, which informs on the influence of internal climate variability on the trends. To do so, a histogram of trends for AoA and tropical upwelling is shown in Fig. 5 using all the members of each model. For upwelling, the same levels as in Fig. 2 are used in order to separate the shallow and deep branches. For AoA, the trends are shown over the regions with available trend estimates from observations, that is, the NH mid-latitude lower stratosphere (80 hPa) and mid-stratosphere (30 hPa). Observational trend estimates of AoA are shown for 145 80 hPa, and displayed in the legend for 30 hPa, as they are outside the range of the abscissa. In addition, the MMM trends from CCMI have been included in all panels. We do not include reanalysis trend estimates for upwelling because the reanalyses do not provide robust trend estimates . In the lower stratosphere there is very good agreement between the CMIP6 MMM AoA trends and the observed trends (Fig. 5a). The trend in the CCMI MMM is also negative but not as strong. Note, however, the reduced number of models with AoA output in CMIP6. In the middle stratosphere ( Fig. 5b) the MMM of 150 both intercomparison projects produce a negative mean age trend between -2 and -3 %/decade, stronger for CMIP6. These values disagree strongly with the observed estimate of +3 %/decade, even when taking the large uncertainty into account.
Even if one considers the updated estimate of AoA trends by Fritsch et al. (2020) of +1.5 ±3%/decade, the range of model estimates barely overlaps with the observational uncertainty range. The residual circulation trends in the shallow branch ( Fig.   5c) range from 0.5 to 3.5 %/decade, with a maximum in the distribution slightly below 2 %/decade, consistent with previous 155 climate model simulations. The CMIP6 MMM trend is in excellent agreement with that from the CCMI MMM (Fig. 5c). In the middle stratosphere, the value of the CMIP6 MMM trend is similar to that in the shallow branch (slightly below 2%/decade), while the trend in the CCMI MMM is weaker (Fig. 5d).
When looking at individual simulations, the AoA trends show a similar spread in the trends across simulations of less than 3 %/decade at the two levels. In contrast, the residual circulation trends show a larger spread in the deep branch than in the 160 shallow branch. In fact, some members feature slightly negative trends (ranging from -0.5 to over 5 %/decade) in the deep branch upwelling. Therefore, a deceleration of the deep branch over 1975-2014 is compatible with the internal variability in some of the CMIP6 models (although in the tail of the distribution), which could be consistent with observational AoA estimates. Nevertheless, we note that negative upwelling trends do not necessarily imply positive AoA trends, because the latter is an integrated quantity, affected non-locally by both advection and mixing (e.g. Garny et al., 2014). Indeed, the GISS and 165  To explore further the role of internal variability, we analyze how the trends depend on the length of the period considered.  (1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985) to 40 years . Trends of individual ensemble members (crosses) are displayed as relative deviations from the ensemble mean trend. In this section we examine the BDC trends and wave forcing in the 1pctCO2 simulations. Figure 7 shows the trend inw * for the 1pctCO2 simulations in the different models and for the MMM. This figure clearly demonstrates the increasing strength of the residual circulation due to CO 2 increase, in both the deep and shallow branches.

BDC trends
The trend is particularly strong in the lower stratosphere and near the stratopause, mirroring the climatological structure (Fig. 1).

190
Changes in the turnaround latitudes indicate that the upwelling region narrows in the lower stratosphere in almost all models.
These features are consistent with previous results (Palmeiro et al., 2014;Hardiman et al., 2014). In the upper stratosphere,  negative trend periods (Fig. 9d). These large oscillations in the deep branch upwelling trends are reflected in the MMM, which shows a quasi-periodicity of about 30 years. In contrast, the upwelling trend in the shallow branch is more consistently positive throughout the period, and shows an increase in the trend magnitude over time, from 2 to 4 %/decade in the MMM (Fig. 9c).

220
The AoA trends are more similar at the two levels, with consistent negative values throughout the period, despite the large oscillations ( Figs. 9a and b). This contrasts with more stable and less uncertain trends in the shallow branch. They also reveal that AoA is a less noisy 225 variable, featuring consistently negative trends in the deep branch across models and members, in contrast to the residual circulation.  Figure 10 shows the contribution of the forcing from different waves to the vertical distribution of the net tropical upwelling for the MMM. To do so, the downward control principle has been applied between the annual-mean climatological turnaround 230 latitudes (Palmeiro et al., 2014). As a novelty from previous multimodel assessments, here we examine the entire stratosphere, including the forcing of the deep branch. The CMIP6 MMM results in Fig. 10a show that the climatological behavior of the tropical upwelling is driven fundamentally by resolved waves throughout the stratosphere. Their contribution is about 70% in the shallow branch, peaks in the middle stratosphere reaching 89% near 7 hPa, and is about 80% in the upper stratosphere.

Drivers of the BDC trends
Parameterized orographic gravity waves contribute 18% to the upwelling at 70hPa, in good agreement with previous multi-235 model studies (i.e. 20% in CCMVal2 models in Butchart et al. (2010)), and their contribution decreases at higher altitudes (6% at 1.5hPa). Parameterized non-orographic gravity waves (NOGW) are not negligible for the shallow branch and account for about 11% at 70hPa, while they become the second contributor in the upper stratosphere (15% at 1.5hPa). Note that model output for NOGW drag was available only for a small number of models in previous multi-model assessments, hampering a direct comparison.

240
As noted above, the vertical structure of the trends in upwelling (Fig. 10b) approximately mirrors that of the climatology (Fig. 10a). As shown in previous assessments (Butchart et al., 2010;WMO, 2014), resolved waves play the primary role in driving trends in the shallow branch. This is due to the intensification and upward displacement of the subtropical jets (not shown) and the upward displacement of the critical lines as discussed in other studies (i.e. Garcia and Randel, 2008;McLandress and Shepherd, 2009;Shepherd and McLandress, 2011;Hardiman et al., 2014). In particular, at 70hPa, the con-245 tribution to the total trend is 63% resolved waves, 25% OGWs and 11% NOGWs. In the upper stratosphere (above 10 hPa), resolved waves and NOGW are equally important to the MMM trends, while the contribution from orographic gravity waves is much smaller. This is in agreement with the results of Palmeiro et al. (2014) for the previous version of WACCM, who explained the key role of NOGW due to changes in the filtering associated with changes in the background winds with increasing greenhouse gases. The resolved waves contribution peaks approximately at 7 hPa with a 59% for the MMM. At that level,

250
NOGW contribute a 32% and OGW a 9%. At 1.5 hPa the percentages are 48%, 41% and 11%, respectively. red: upwelling due to resolved waves, computed from the Eliassen-Palm flux divergence (DELF); green: upwelling due to orographic gravity waves (OGW); blue: upwelling due to non-orographic gravity waves (NOGW). Note that the GISS model is not included because the NOGW output was not available. Figure 11 shows the intermodel spread in wave forcing. The forcing is more consistent across models for the climatology (panels a and b) than for the trends (panels c and d), in agreement with Butchart et al. (2010) results in the lower stratosphere (see also WMO (2014) and references therein). In the shallow branch, most models (except one) attribute the trends primarily to resolved wave drag. In addition, all the models show a relatively small contribution from NOGW (less than 20%). Finally, 4 out 255 of 7 models feature a larger contribution from OGW than NOGW (with OGW being the main forcing for GFDL-ESM4). Note that a present-day climatological source of NOGW is launched at 70 hPa in the extratropics in MIROC6 (Watanabe, 2008), which may explain the small contribution to the trends at this level and the negative contribution in the upper stratosphere.
The deep branch trends feature much larger spread, with a factor of 4 difference between the smallest trends (MIROC6) and the largest trends (UKESM) at 1.5 hPa. The larger role of NOGW compared to the shallow branch is clear, although there is 260 quite a range in the contributions of resolved and NOGW. Four out of 7 models show comparable contributions, resolved wave forcing dominates in GISS-E2 and MIROC6 while in CESM2-WACCM there is a comparable contribution from NOGW and OGW, and a negligible contribution from resolved waves at 1.5 hPa. These results highlight a wide diversity of forcings of the CO 2 -driven trends among models in the deep branch.
Note that there is no proportionality between climatologicalw * and its trends across models, in contrast with results from Note that for the GISS model NOGW output is not available, although these waves are present in the model. in the range from 5 to 13% per degree of surface warming. The strong connection between the shallow branch of the BDC and tropical surface temperature is consistent with the findings from previous works Chrysanthou et al., 2020;Orbe et al., 2020).

280
In contrast, the deep branch sensitivity to surface temperature is near zero in the historical runs, while in the 1pctCO2 runs it is small but consistently positive across the models. This difference between the two simulation types is likely due to the contribution of ozone depletion to the deep branch trends in the historical run, shown in Fig. 4. The impact of ozone depletion on the BDC (present in the historical but not in the 1pctCO2 simulations) is independent from surface temperature, and therefore it reduces further the -already weak-connection between deep branch trends and surface warming. Note that the 285 behavior is very similar when the global or the tropical mean surface temperature trends are considered. The disconnection between the deep branch of the BDC and surface temperature is consistent with the study of Chrysanthou et al. (2020), which finds that the acceleration of the deep branch is largely due to the direct CO 2 radiative effect.
The connection between surface temperature and BDC is explored at interannual and decadal timescales in Figure 13. The decadal regression coefficients (panel b) show a very high consistency with the trend behavior in Fig. 12, both qualitative and 290 quantitative. The corresponding coefficients of determination (R 2 , panel d) reveal that the fraction of variance of the shallow branch tropical upwelling explained by the surface temperature is around 35-60%, with a wide inter-model spread (ranging from 15 to 85%). These values are about a 15% higher for the 1pctCO2 than for the historical simulations, and also for the tropical than for the global surface temperature.
The results for the interannual variations are different in that the regression values for the shallow branch are notably higher 295 (approximately twice as large) for global than for tropical surface temperatures (Fig. 13a). This likely reflects the fact that interannual temperature variations often have an antisymmetric pattern between tropics and extratropics. In particular, the connection between the BDC and surface temperature on interannual timescales is dominated by ENSO (e.g., Calvo et al., 2010), which features strong sea surface temperature anomalies (SST) confined to the tropics, and opposite sign anomalies in the extratropics. The globally averaged surface temperatures attenuate the signal (associated with ENSO) that is mainly G lo b a l, 7 0 h P a T ro p ic a l, 7 0 h P a G lo b a l, 1 .5 h P a T ro p ic a l, 1 .5 h P a 20 In contrast, long-term trends and low-frequency variability have a more uniform latitudinal pattern. On the other hand, the fact that explained variances (R 2 ) are comparable for the tropical as for the global surface temperature, for both interannual and decadal timescales, implies that the surface warming signal on the lower branch of the BDC is controlled primarily by tropical temperature.

Conclusions and outlook
The present paper examines the BDC in CMIP6 models, focusing on the residual circulation and the mean age of air.
First, historical simulations are used to compare the climatology and past trends to observations and reanalyses. The climatological average and seasonality of the BDC in CMIP6 models lie within the spread of observational (and reanalysis) estimates.
However, there is a large spread in the magnitude among models, which is as large as that across observational estimates. The 310 BDC trends in historical simulations are stronger during the ozone depletion period than after, which reflects the important contribution of ozone depletion to BDC acceleration shown in previous chemistry-climate model studies. In agreement with previous multi-model studies, there remains a clear disagreement in AoA trends between models and observations in the middle and upper stratosphere. On the other hand, there is very good agreement in the shallow branch trends. The trends in the deep branch of the residual circulation reveal a large spread among models and across different members of the same model (Figs. 5,6). In particular, the inter-member spread in deep branch trends is about ±200% of the ensemble mean trends, even when considering long (40 year) periods. In contrast, for the shallow branch the spread is 4 times smaller. This reveals a notably stronger influence of internal variability on the deep branch than the shallow branch trends. In addition, the trends are more robust for AoA than for upwelling, with inter-member spread in the trends below ±30% for periods slightly longer than 20 years, both in the lower and middle stratosphere.

320
The sensitivity of the BDC to CO 2 increase is examined, and the robustness of the trends and their wave forcing are explored.
In contrast with previous BDC analyses based on multi-model assessments with CCMI and CCMVal, we focus here on the response to CO 2 alone, using the 1pctCO2 simulations. All models produce stronger BDC acceleration in the NH than in the SH, and the BDC actually decelerates in the SH polar lower stratosphere, possibly due to the CO 2 effects on the polar vortex discussed in recent works (e.g. Ceppi and Shepherd, 2019). An analysis of the wave forcing of the residual circulation 325 shows that shallow branch forcing of climatology and trends is mainly due to resolved waves with a contribution from OGW, consistent with previous studies. For the deep branch, the main drivers of climatology and trends are resolved waves and NOGW, but there is a wide spread across models, especially for the trends. There is a very large uncertainty in deep branch trends (factor of 4 inter-model spread), which could be linked to the spread in forcing. In contrast, the spread in the shallow branch trends is less than 30%. The uncertainty in deep branch residual circulation trends is emphasized in the large multi-330 decadal fluctuations found over the 150 simulation years. On the other hand, the shallow branch trends are found to increase over time with CO 2 increase, by approximately a factor of 2 for the MMM. In contrast, the AoA trends are more robust over time, consistent with the results for the historical simulations.
Finally, the connection between surface temperature and the BDC is investigated. We find a strong connection between the shallow branch and the tropical and global surface temperature. Long-term trends in lower stratospheric upwelling feature a 335 sensitivity of 7-10% per degree of surface warming in the models (Fig. 12). On interannual and decadal timescales, surface temperature explains 35-60% of the shallow branch variance on average (Fig. 13). Note that the strong connection of shallow branch acceleration with surface warming is consistent with the correlation with the upward shift of the tropopause pointed out by (Oberländer-Hayn et al., 2016). In contrast, the deep branch variability is not correlated with surface temperature on any timescale.

340
One of the key results of the present paper is the difference between shallow and deep branches of the residual circulation.
The CMIP6 models confirm that, while trends in the shallow branch can now be reconciled with observations (Fu et al., 2015; WMO, 2018), a clear inconsistency remains for the trends in the deep branch. Our analyses reveal much larger uncertainties in the deep than in the shallow branch trends, associated with larger internal variability as seen by the inter-member spread. We note that, while a robust mechanism for the acceleration for the shallow branch has been described (Shepherd and McLandress,345 2011), the drivers of deep branch acceleration remain largely unexplored. Previous studies point to the effects of stratospheric zonal wind trends on the filtering of NOGW. The CMIP6 model results in the present paper confirm the important role of NOGW for the deep branch trends. The zonal mean wind trends show acceleration of the polar jet in both hemispheres for the MMM, but there is a large spread in the trends in the NH, with some models featuring deceleration (not shown). This is con-sistent with large differences in resolved versus NOGW forcing contributions to the deep branch among models. However, the compensation mechanism (Cohen et al., 2013;Sigmond and Shepherd, 2014) implies that the different relative contributions from resolved versus parameterized gravity waves does not necessarily lead to differences in the net BDC trends. Overall, open questions remain regarding the deep branch trends and their forcing mechanisms. Reducing uncertainties in deep branch trends is particularly relevant to better constrain the future distribution of ozone in the polar stratosphere, affected not only by direct transport but also by the descent of ozone-depleting chemical compounds from the mesosphere (Maliniemi et al., 2020). 355 The results show that AoA is a much less noisy variable thanw * , implying that robust trends could be extracted from relatively short periods (20 years). The advective circulation can be approximated from AoA using the leaky pipe model, in order to compare with observations, as done in Linz et al. (2016). Unfortunately, the global AoA observational estimates available are not long enough to evaluate trends. In addition, it is crucial to account for the large uncertainties in deriving AoA trends from realistic tracers (Fritsch et al., 2020). We note that, for the few models providing both AoA andw * (3 models for 1pctCO2 360 and 4 for historical simulations), it is not possible to extract robust conclusions on the relationship between the strength of the climatological values or the trends of the two variables. Establishing relations between these two magnitudes would help evaluate the spread in the magnitude and variability of mixing across models (e.g., Dietmüller et al., 2017;Eichinger et al., 2018). Based on the results of the present study that highlight the robustness of AoA trends, we suggest that AoA should be a first-priority consistently defined diagnostic for the next CMIP project (Gerber and Manzini, 2016).  Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer,