The influence of mixing on stratospheric age of air changes in the 21st century

Climate models consistently predict an acceleration of the Brewer-Dobson circulation (BDC) due to climate change in the 21st century. However, the strength of this acceleration varies considerably among individual models, which constitutes a notable source of uncertainty for future climate projections. To shed more light upon the magnitude of this uncertainty and on its causes, we analyze the stratospheric mean age of air (AoA) of 10 climate projection simulations from the Chemistry Climate 5 Model Initiative phase 1 (CCMI-I), covering the period between 1960 and 2100. In agreement with previous multi-model studies, we find a large model spread in the magnitude of the AoA trend over the simulation period. Differences between future and past AoA are found to be predominantly due to differences in mixing (reduced aging by mixing and recirculation) rather

same story. The time series of the observations presented in the studies by Engel et al. (2009);Ray et al. (2014); Engel et al. (2017) show a (non-significant) positive trend in the Northern Hemisphere (NH) across the recent decades, but most climate models show an AoA decrease over time. Also the trends in the (much shorter) satellite time series mostly do not coincide with the model results (Ploeger et al., 2015a). This discrepancy is still an ongoing debate and it has been addressed in numerous studies (e.g. Garcia et al., 2011;Bönisch et al., 2011;Stiller et al., 2017;Lossow et al., 2018). 5 In the present study, we focus on the AoA differences between future and past simulated by 10 chemistry-climate models that participated in the Chemistry Climate Model Initiative phase 1 (CCMI-1; Morgenstern et al., 2017). We analyze the changes of stratospheric AoA in the 1960-2100 climate projection simulations between the two periods 1970-1990and 2080-2100 is well established that in climate change simulations models predict an acceleration of the BDC, which consequently leads to younger stratospheric air (e.g. Rind et al., 1990;Butchart and Scaife, 2001;Garcia and Randel, 2008). Stratospheric transport is 10 therefore sensitive to varying greenhouse gas concentrations, but also other constituents like ozone depleting substances (ODS) have been shown to play a considerable role in modulation of stratospheric transport (e.g. Oman et al., 2009;Oberländer-Hayn et al., 2015;Polvani et al., 2017Polvani et al., , 2018. ODSs on the one hand, act as greenhouse gases themselves, and on the other hand, lead to the chemical destruction of stratospheric ozone. However, multi-model intercomparison studies have revealed that the magnitude of the BDC acceleration until the end of the century varies strongly among the various state-of-the-art climate 15 models (Butchart et al., 2006(Butchart et al., , 2010Hardiman et al., 2014). Moreover, the mechanisms driving these changes are still not entirely clear.
To analyze the reasons for the AoA changes and their differences between various chemistry climate models (CCMs), we follow the approach of Birner and Bönisch (2011) who calculated the residual circulation transit times (RCTTs) by means of backward trajectories. Garny et al. (2014) have then separated AoA into the two contributions of residual transport and 20 aging by mixing. In several previous studies (e.g. Ploeger et al., 2015a;Dietmüller et al., 2017), it has been shown that this concept is well-suited for process-based model analyses Here, the method allows us to conclude that a large fraction of the change in the stratospheric circulation is due to aging by mixing, and residual transport only regionally plays the primary role.
But these conclusions can prove fallacious, because also interactions between these two processes (Ploeger et al., 2015a) have to be considered. While the changes that originate from residual circulation changes mainly depend on the strengthening of 25 tropical upwelling (see e.g. Randel et al., 2008;Butchart et al., 2011;Butchart, 2014), the origin of changes in mixing are widely uncharted. We therefore calculate the mixing efficiency, an independent measure for the relative strength of mixing for given residual mean transport changes (ratio of mixing mass flux to net mass flux) (Garny et al., 2014), across the 21st century by means of a one-dimensional transport model of the stratosphere in the CCMI-1 model simulations. In the companion paper, Dietmüller et al. (2018) have already shown that the mixing efficiency can explain most of the AoA model spread 30 in the climatologies from 1960 to 2010. In the present study, we quantify the impact of mixing efficiency (relative mixing strength) differences between two climate states in the individual model simulations. Moreover, we show the influence of mixing efficiency variations on the model spread in AoA changes. To conclude, we also provide a theoretical explanation for the reasons of the relative mixing changes, based on the role of the ratio of wave dissipation to potential vorticity gradients in controlling the mixing properties.

CCM simulations
In this study, we analyze the model output of 10 state-of-the-art CCM simulations. All these simulations were conducted in the framework of the Chemistry Climate Model Initiative phase 1 (CCMI-1, Morgenstern et al., 2017). An overview on the model simulations that are used for analysis in this study and some aspects of their model setup is given in Tab. 1. Additionally, the 5 name of the atmospheric model component is provided to demonstrate similarities between some of the models. This sub-set of models has been chosen on the basis of availability of the required data for the analyses in this study (AoA and residual circulation velocities).

Analysis methods
The methodology of this study mostly follows the companion paper Dietmüller et al. (2018). A short description of the most important concepts used here is given in the following. 15 Stratospheric mean AoA is defined as the mean residence time of an air parcel in the stratosphere (Hall and Plumb, 1994;Waugh and Hall, 2002). In the CCMs, the AoA tracer is implemented as an inert tracer with a mixing ratio that linearly increases over time as lower boundary condition ("clock tracer" Hall and Plumb, 1994). In some models, this lower boundary condition is global, in others only in the tropics, in NIWA-UKCA and ACCESS for example, AoA is kept at 0 in the boundary layer. AoA is then calculated as the time lag between the local mixing ratio at a certain grid point and the current mixing ratio at a reference 20 point. This reference point, however, varies between the models (e.g. boundary layer, tropopause, 100 hPa), which could lead to inconsistencies in the AoA calculation. To avoid this, we subtract the mean AoA value of the respective model's tropical tropopause from the actual AoA value at all grid points. This means that in our analyses, AoA=0 at the tropical tropopause for all models. We perform this calculation for each time step (with monthly values) separately for the entire time period 1960-2100. Therefore, the AoA trend exludes any changes in tropospheric transport times due to the fact that the tropopause rises 25 over time (see Vallis et al., 2015;Oberländer-Hayn et al., 2016;Abalos et al., 2017).
The residual circulation transit time (RCTT) is the hypothetical age that air would have if it only followed the residual circulation, meaning that no processes such as eddy mixing or diffusion would come into play. These RCTTs are calculated using a concept described by Birner and Bönisch (2011), by calculating backward trajectories on the basis of the Transformed Eulerian Mean (TEM) meridional (v * ) and vertical (w * ) velocities (referred to as residual velocities, available in the CCMI-1 30 data base for the chosen models) with a standard fourth-order Runge-Kutta integration. The RCTT is then the time that these backward trajectories require to reach the tropopause from the respective starting point in the stratosphere. For more details see also Birner and Bönisch (2011) and Garny et al. (2014).
The RCTT differs from AoA because of (resolved and unresolved) mixing. In the stratosphere, this is due to the mixing of air between branches and the in-mixing of air from the mid-latitudes into the tropical pipe, which leads to recirculation of old air around the BDC branches. In global model studies, this effect has been named aging by mixing (A_mix) and is interpreted as the difference between AoA and RCTT (Garny et al., 2014;Ploeger et al., 2015a, b). However, it has to be kept in mind that the residual of AoA and RCTT does not only reflect this process alone, but actually includes resolved mixing as well as 5 parameterized and numerical diffusion.
As a measure of the relative strength of mixing (independent of the residual circulation strength), we use the so-called mixing efficiency for analysis. is defined as the ratio of the mixing mass flux to the net residual mass flux between the tropics and the extratropics across the subtropical barrier. The net mass flux is the horizontal motion that is determined by mass continuity from vertical motion and corresponds to transport by v * (Garny et al., 2014). can be derived by means of 10 the tropical leaky pipe (TLP) model (Neu and Plumb, 1999). The TLP model is a one-dimensional transport model of the stratosphere that includes advection and horizontal mixing of air between the tropics and the extratropics. It assumes two columns of well-mixed air (a tropical and an extratropical column) and can be used to quantify the strength of mixing across the subtropical barrier. If we neglect vertical diffusion (see Dietmüller et al., 2017), we can formulate an analytical solution for tropical and mid-latitude AoA. According to the TLP model, tropical AoA (AoA T ) with altitude-dependent vertical velocity (1) (Neu and Plumb, 1999;Garny et al., 2014). Here, H denotes the scale height (7 km), α(z) the ratio of tropical to extratropical mass approximated by 20 z t the height of the tropical tropopause and T corr (z) = H 1 w * T (z) − 1 w * T (zt) the correction term for the altitude-dependency of the vertical residual velocity w * (additional analytical solution term from horizontal advection, for details see Neu and Plumb (1999); Garny et al. (2014)). AoA thus depends on the advective vertical velocity w * T (i.e. the residual velocity) and on the mixing efficiency (i.e. the amount of mixing between the tropics and the extratropics). Solving Eq. 1 for the mixing efficiency Eq. 3 shows that is approximately proportional to the relative increase in AoA due to mixing. For analysis, we calculate the ten year running averages of to obtain climatologically representative values. Note that according to the concept of the TLP model, has to be viewed as a parameter valid for a given climate state in a certain model. Hence, it can vary between models and/or for different climate conditions. In the study by Garny et al. (2014), the authors nevertheless found a constant for three different climate states in one model. The tropical profiles of w * , tropopause height and AoA provided for the TLP model are averaged over the latitudinal bands of the models' individual vertically averaged turnaround latitudes (which are also time-dependent; calculated from v * for consistency between the models, refer to Supplement of Dietmüller et al. (2018)) and are interpolated to vertical coordinates relative to the tropopause height of each model. The mixing efficiency is then obtained 5 by the TLP model's best fit to the CCM AoA profile. Here, the best fit is done for the altitude range from the tropopause to 32 km (details for the calculation of the mixing efficiency are given in Garny et al., 2014). According to the TLP formulation, aging by mixing (A_mix) is a function of the mixing efficiency and of the residual circulation strength: A_mix is proportional to the mixing efficiency, but inversely proportional to the vertical velocity. A higher mixing efficiency 10 is e.g. connected with more air parcels recirculating (see Garny et al., 2014), thereby increasing A_mix. But also the vertical velocity controls the speed of the air parcels that recirculate. Thus, the mixing efficiency has been shown to be a useful diagnostic tool, as it does not depend on the speed of recirculation (e.g. Garny et al., 2014;Dietmüller et al., 2017).

Changes in AoA and in its components
A multi-model comparison of stratospheric transport changes in the 21st century has been conducted before for the Stratosphere-Troposphere Processes and their Role in Climate (SPARC) report (SPARC, 2010). In that study, the authors showed stratospheric mean AoA of 10 chemistry-climate model simulations that took part in the CCMVal-2 (Chemistry-Climate Model Validation, Eyring et al., 2013) project. To allow a direct comparison of the simulations that were analyzed in SPARC (2010) 20 with the simulations we use here, Fig. 1 presents the same depiction of AoA differences between the 2090s and the 1990s (∆AoA) a) at 50 hPa and b) as difference between tropics and middle latitudes as in their figure 5.18.
The general structures of the AoA differences agree among the two model intercomparison projects. All models predict a decrease in mean AoA at 50 hPa and the smallest decrease in the tropics (Fig. 1a). A second minimum of the decrease is found at 60 • N/S and the greatest decrease at 30 • N/S and/or at the poles. The AoA behaviour in these regions of maximum change 25 across the 21st century is investigated in detail by Šácha et al. (2018). They showed that the trends there are related to the climatological AoA distribution, the upward shift of the pressure levels and the widening of the AoA isolines. In comparison with the CCMVal-2 models, the inter-model spread of the AoA difference is reduced in our results. However, it were the two UMUKCA (Unified Model/U. K. Chemistry Aerosol) models that lead to the large spread in the SPARC report and are not part of our analysis. However, NIWA-UKCA and ACCESS are the direct successors to the UMUKCA models and these range 30 rather in the lower end of AoA changes here. The ULAQ model has changed from a very large latitudinal ∆AoA amplitude to a rather small one from CCMVal-2 to CCMI-1. Other models that appear in both studies do not show large changes. The EMAC of the CCMVal-2 and the CCMI-1 projects, Dietmüller et al. (2018) showed that it is mainly the mixing rather than the residual circulation that causes the large AoA spread and that this is likely linked to the different resolutions of the model simulations.  Engel et al. (2009Engel et al. ( , 2017. Note that the observational AoA data is relative to ground level, while the model data has been processed to be relative to the tropopause (see main text).
All the model simulations show a clear negative trend over time across the first as well as across the second period. The in situ measurement-based observations (Engel et al., 2009(Engel et al., , 2017, in contrast, display a positive, but not significant, trend for the recent couple of decades. A similar behaviour as in the in situ measurements can also be found in satellite-based observations for example in Haenel et al. (2015), although for a shorter time series (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). A number of studies have investigated this discrepancy (see e.g. Engel et al., 2009;Garcia et al., 2011;Ploeger et al., 2015a) and many reasons have been discussed 5 to resolve this contradiction (e.g. effect of mixing in models, sparse sampling of observational data, differences in changes between deep and shallow BDC branches), but a viable explanation is still missing. In the present study, however, we do not discuss this issue any further, but rather focus on the analysis of the negative AoA trends in the model simulations. The models may agree in predicting a decrease in stratospheric AoA, but they do show large differences in the strength of this trend.
(σ(∆AoA) = 0.18 for µ(∆AoA) = 0.54, see below for ∆AoA) Also, the trends between the two periods (1960-2000 and of ODSs for the trends in stratospheric dynamics. ODSs act as both, radiatively active greenhouse gases and chemically active gases controlling ozone depletion and recovery. During the period between 1960 and 2000, ODSs increase and thus act in concert with GHGs to cause an AoA decrease. Thereafter, however, ODSs decrease over time, which means that with respect to AoA trends, the ODS trend works against the trend in the continuously rising GHGs. This can possibly explain the weaker trend in the second period in most models . Due to this change in dynamical properties around the year 2000, an analysis of the trends across the entire period from 1960 to 2100 cannot be conducted without mixing-up various dynamical effects. The first period is relatively short for robust analyses of mixing trends and in the second period, two mechanisms work against each 5 other, so that in most models, the trends are rather small. In the following, we therefore analyze the differences between the periods 1970-1990 and 2080-2100. Based on the sensitivity simulations in Polvani et al. (2018), which shows that the change in the slope in the year 2000 is due to ODSs, this allows to capture the GHG effect alone. At the end of the 21st century, ODS mixing ratios have declined to similar values as between 1970 and 1990 in the simulations. We do not start our investigation at 1960 due to the calculation method of RCTTs and which are available only from 1970 onwards. and over 10-100 hPa, which means that these AoA values are not comparable with those in Fig. 2. Tab. 2 will be used to explain and confirm some of the results throughout the paper. Specifically, values between 0.5 and 1 signify a domination of residual transport changes (reddish), while values between 0 and 0.5 mean that ∆AoA is controlled by A_mix variations (blueish). Values above 1 mean that the A_mix difference is positive and values below 0 mean that ∆RCTT is positive (assuming a negative ∆AoA).
In most models, a quite similar ∆AoA pattern can be seen. Relatively small ∆AoA dominates the tropical pipe where young air is prevalent and the differences increase with higher altitudes and latitudes. Several model simulations (ACCESS,  Fig. 3 suggests that most of the ∆AoA is determined by aging by mixing changes. However, A_mix itself depends on RCTT (see Eq. 4) and is therefore not an independent measure for separation of the processes. AoA changes are commonly attributed to changes in the residual circulation (e.g. Austin and Li, 2006). In the climatologies, the tropical upwelling (as calculated via 25 integration of v * between the turn-around latitudes, see e.g. Okamoto et al., 2011) in the 70 hPa pressure level has been shown to be a good measure for the strength of the residual circulation throughout the stratosphere (see Dietmüller et al., 2018). To see if this relationship holds true for ∆AoA in our set of the CCMI model simulations, we present in Fig. 4 the inter-model correlations between ∆AoA and ∆RCTT as well as between the ∆AoA and the differences in tropical upwelling in 70 hPa.

Coherences between the components
These correlations are calculated across the 10 model simulations for each grid point separately. Note that for this, data of each 30 model was interpolated to the resolution of the grid of the highest horizontal model grid resolution.
A high correlation between ∆RCTT and ∆AoA (Fig. 4a) can only be seen in the middle to high latitudes between around 70 and 10 hPa and low or no correlations in the other regions. Particularly in the upwelling region in the tropics this is rather surprising since AoA should be controlled by the upwelling speed there. However, this connection apparently is not local. The  that the link between the residual circulation and AoA cannot be expected to be local, but rather of remote nature. It is mainly the ascent of air in the tropics that determines the ∆AoA-dependency on the residual circulation. Particularly the upwelling and downwelling regions do not correlate well in Fig. 4a, although ∆AoA is clearly dominated by the changes in transport in those regions (see Fig. 3). 10 Dietmüller et al. (2018) showed that the climatological inter-model AoA spread in the CCMVal-2/CCMI-1 hindcast simulations can mostly be explained with model differences in mixing efficiency and only to some extent with differences in residual circulation. Hence, we now focus on the model trends (and differences between the two periods) in mixing efficiency and on the processes that are responsible for its changes. provides information on the models' relative strengths in mixing at the tropical barriers as a vertically integrated measure (see Sect. 2.2 for calculation method).  Sect. 3.4. Next, we investigate if ∆AoA is also linked to variations in mixing efficiency, and not only to residual circulation changes. For this, Fig. 6 shows the inter-model correlations (correlations across the 10 model simulations) between the local ∆AoA and the models' differences in mixing efficiency.  weaker. This may be due to changes in the polar vortex strength or position in the various models, which cannot be reflected in the (sub)tropical measure .
A_mix is connected with the mixing efficiency, but it also depends on the residual circulation. A higher/lower mixing efficiency causes more/less recirculation, which increases/decreases aging by mixing, and the RCTT controls the transit times 10 of the air parcels that recirculate. Overall, we can again conclude that ∆AoA in climate models is connected with changes in both, mixing and residual transport. Most of the ∆AoA that is connected with residual transport can be attributed to tropical upwelling changes, at least for the deep BDC branch. This is caused by a climate change induced strengthening and/or upward shift of the subtropical jet streams (see e.g. Randel et al., 2008;Shepherd and McLandress, 2011;Butchart, 2014;Oberländer-Hayn et al., 2016). The differences in mixing between future and past climate, however, have not been analyzed in detail before. 15 Therefore, we investigate the overall effect of mixing on the BDC changes in the remainder of the paper.

Impact of the mixing changes
To quantify the impact that changes in mixing efficiency have on ∆AoA in the simulations, we again apply the TLP model, now by using w * , the tropopause height and a given to calculate the RCTTs and AoA (see Eq. 1). We then compare the AoA and RCTT climatologies of the two periods 1970-1990 (1970) and 2080-2100 (2100). Fig. 7 shows the tropical AoA and RCTT profiles (calculated with the TLP model) of these climatologies exemplarily for the EMAC-L90 model simula- The difference (about two thirds) is caused by aging by mixing (see Fig. 7). This ratio had already been found in Dietmüller et al. (2018) for the CCMVal-2 and CCMI-1 hindcast simulations. AoA' reveals that for EMAC-L90 this 15 ratio roughly also counts for the differences between the periods. The AoA difference between 2100 and 1970 is subdivided by AoA' into about one third the share of RCTT difference and two thirds of A_mix. This ratio can also be estimated from Fig. 3 for the EMAC-L90 simulation. However, as already stated above, the A_mix change also includes some RCTT change (see Eq. 4), which implies that this is not a clear separation of the mechanisms.
The difference between AoA*(2100) and AoA(2100) reflects the impact of ∆ on AoA. In the given example of the EMAC-20 L90 model, the fractional impact of the mixing efficiency change on the AoA difference is 29±4% as calculated between 2.5 and 25 km above the tropopause. This value varies considerably among the ten models, Tab. 3 provides an overview. NIWA-UKCA, ACCESS and ULAQ show a negative fractional impact of mixing efficiency on AoA changes. These are the three models that also show a positive ∆ (see Tab. 2). In contrast to the other models, the mixing efficiency therefore leads to an AoA increase over time. The negative ∆RCTT therefore accounts for more than the entire negative ∆AoA to compensate 25 the effect of the change. With less than 5.5% and 3.5%, NIWA-UKCA and ACCESS have the lowest contribution of ∆ on the AoA change and with up to 29%, ULAQ and EMAC-L90 have the largest. The other models have a contribution of ∆ on ∆AoA between 10 and 29%, the multi-model mean is 10.4% (σ = 18.3%). It is reasonable that a large ∆ leads to a high percentage of the -share on ∆AoA. The correlation coefficient between the two values is 0.83. This makes clear that the change in mixing efficiency does have a considerable impact on ∆AoA, as it could already be assumed from the high inter-model correlations between the two quantities (Fig. 6). Due to the large model spread in ∆ , however, this impact bears 5 large uncertainties. Next, we quantify the impact of ∆ on the model spread in ∆AoA.
To analyze the share of ∆ on the AoA model spread, we again use the TLP model-calculated values of , α(z), T corr , and RCTT T . These quantities are taken to calculate tropical AoA by 10 AoA as an average between 100 and 10 hPa. As above, this calculation is performed for each model for the two periods, respectively, so that ∆AoA T = AoA T (2100)−AoA T (1970). Using (1970) to calculate AoA T (2100) in Eq. 1 provides AoA * T (2100) and thus ∆AoA * T = AoA * T (2100) − AoA T (1970) provides the AoA difference without any changes in . Fig. 8 displays the model distribution of ∆AoA with (blue) and without (red) a changing mixing efficiency.  . This reflects that the mixing efficiency changes lead to strong variations in AoA. When is held constant, for models with negative ∆ , ∆AoA is reduced (in absolute values) (see e.g. EMAC) and ∆AoA is enhanced for models with positive ∆ (see e.g. ULAQ). The model range decreases here because being the model with the largest ∆AoA, EMAC-L90 has a negative ∆ and ULAQ, the model with the lowest ∆AoA has a positive ∆ . A large/small (negative) ∆ causes a large/small ∆AoA. From this analysis, we can conclude that |∆AoA| 10 generally increases through the changes in mixing efficiency (by 10.4% as multi-model mean) and that ∆ leads to a larger model spread in the AoA changes (∆ increases the ∆AoA model range by about 37%). Now, the question remains what the reasons for the changes in relative mixing strength are, or why changes in the course of climate change simulations. In the next section, we will explore an explanation for this by analyzing various dynamical fields.

On the mechanism of mixing changes
The relation of the residual circulation and mixing determines the mixing efficiency and changes therein. To study possible dynamical reasons for the changes in mixing efficiency, the relation of wave driving, the residual circulation and mixing is analyzed in the following based on the transformed Eulerian mean (TEM) momentum equation. This analysis is similar to that presented in Garny et al. (2014), but here we use the quasi-geostrophic (QG) formulations on pressure levels (on which the 5 CCMI data are available). We present and analyze multi-model mean (MMM) diagnostics to provide a general picture for the entire subset of CCMI-1 model simulations. Unless otherwise stated in the text, the individual models qualitatively agree fairly well in these diagnostics. The fields of the individual models can be found in the Supplement. Fig. 9 shows the differences between the two periods (1970-1990 and 2080-2100) overlaid with the MMM climatologies of the first period of zonal winds, Eliassen-Palm flux divergence (EPfd), v * , the meridional potential vorticity gradient (∂PV/∂y), as well as K yy and K yy /|v * |. 10 The calculation of K yy and the meanings of these variables will be explained below. Due to data availability, not all models could be included in these MMMs. For SOCOL and ULAQ, the EPfd fields were not available, hence, for consistency, all MMMs base on the remaining eight other models only.
The mechanism of the residual circulation increase in the lower stratosphere is well understood (e.g. Shepherd and McLandress, 2011). It follows the upper tropospheric temperatures increase that leads to an upward shift and increase of the zonal 15 mean winds. This can be seen in Fig. 9a, which shows the MMM differences of the zonal mean wind between the two time periods. Subsequently, the critical layers that allow for wave propagation shift to higher levels. This becomes evident from Fig. 9b. There, a region of strongly negative ∆EPfd (enhanced wave dissipation) can be seen between about 100 and 50 hPa (with maxima around 70 hPa) and positive differences (less wave dissipation) can be found between around 200 and 150 hPa.
The enhanced wave dissipation in the lower stratosphere leads to amplified poleward residual transport, which reflects in ∆v * 20 (Fig. 9c).
However, as shown in the sections above, transport changes in a future climate are not only due to this increase in residual transport, but also due to changes in eddy mixing. The EPfd (∇ · F ) not only drives the residual circulation, but it is also related to eddy mixing of potential vorticity q. Under the quasi-geostrophic (QG) approximation and neglecting parameterized (gravity) wave forcing, this can be formulated as where v and q denote the deviation of their zonal means v and q, respectively. Note that due to data availability, we show the full EPfd fields in Fig. 9 and not the QG EPfd. Note also, that strictly this relation holds true on the beta plane (Andrews et al., 1987), or approximately on isentropic surfaces, i.e. on Θ levels, on which mixing takes place. Therefore, the quantity may be somewhat distorted on the pressure levels presented here. This may lead to differences in the quantities, however, not 30 qualitatively different conclusions in the following. Commonly, a flux-gradient relationship is assumed for the eddy PV flux, so that v q = −K yy ∂q y , with the diffusivity coefficient K yy . ∂q y is the gradient of the QG PV (see e.g. Edmon et al., 1980; Similarly as in Abalos et al. (2016), we calculate the diffusivity coefficient K yy as This relation states that horizontal eddy mixing is proportional to the EP flux divergence, and inversely proportional to the 5 meridional PV gradient. Thus, a weak PV gradient indicates strong mixing, while a strong PV gradient acts as barrier to mixing, and thus mixing is suppressed.
The MMM differences show that next to the enhanced wave dissipation in most of the stratosphere (−∇ · F increases), also the PV gradient increases in the subtropics to mid-latitudes (see Fig. 9d) due to the strengthened winds. Thus, K yy tends to be enhanced due to increased wave dissipation, but reduced due to the increased PV gradient. As seen in Fig. 9e, the 10 diffusivity changes are dominated by the increase in wave dissipation in the lower stratosphere (between around 100 to 50 hPa in both hemispheres) and dominated by the decrease in wave dissipation below. This means that the vertical shift of EP flux convergence is reflected in ∆K yy , or in other words, that the region of strong mixing is shifted upward and the absolute strength of mixing increases in the lower stratosphere. At higher altitudes, the PV gradient contributes more strongly to K yy , and in the SH mid-latitudes, the mixing strength K yy even decreases despite enhanced wave dissipation due to the strongly enhanced PV 15 gradient. In the NH, ∆K yy is smaller and mostly positive in the MMM. However, it should be noted that particularly in the lower stratospheric subtropics, the region of the climatological maxima of orographic gravity wave drag, the QG assumption can be misleading due to the neglection of parameterized wave drag.
The mixing efficiency derived from the AoA data represents the relative strength of mixing, i.e. the ratio of the mixing mass flux to the residual meridional mass flux. The mixing mass flux is proportional to the mixing velocity that can be expressed as 20 v mix = K yy /L for a given horizontal length scale L. Thus, the ratio of mixing versus residual mass flux can be approximated Given a constant length scale L, the mixing-to-mean-advection ratio is proportional to K yy /|v * |. As shown above, residual transport as well as mixing increases in most regions ( Fig. 9c and e). Fig. 9f, however, shows that the ratio K yy /|v * | decreases 25 in large parts of the stratosphere. This means that mixing increases less strongly (or even decreases) than residual transport.
Only in the (sub-)tropical lower to mid stratosphere, the relative mixing strength increases. The changes in K yy /|v * | mostly reflect the inverse change in the PV gradient, because the mixing diffusivity is inversely related to it. Note that if assuming v * = −f −1 ∇ · F , the ratio K yy /|v * | equals f /∂q y , i.e. the ratio is only related to the PV gradient. Changes in f /∂q y are overall similar to changes in K yy /|v * | (not shown), except for some details that might be related to the simplification of the equations and/or to the neglection of parametrized wave drag (i.e. small scale gravity wave drag).
Overall, the analysis presented above suggest that enhanced wave dissipation (caused by the zonal mean wind increase, see Shepherd and McLandress, 2011) amplifies the residual circulation as well as mixing, particularly in the lower stratosphere.
However, the zonal mean wind changes also cause changes in the meridional PV gradient, and the stronger PV gradients act to inhibit mixing. Therefore, mixing increases less strongly than residual transport does, which explains the decrease in the mixing efficiency in the future (see Sect. 3.2).

5
However, as shown in Sect. 3.2, the mixing efficiency does not decrease in all models. In Fig. 10, the mean ∆K yy /|v * | in the region of the subtropical barrier averaged between 20 • S-40 • S and between 20 • N-40 • N is shown with altitude for the subset of eight CCMI-1 model simulations. |∆(K yy /|v * |)| is larger in the SH as compared to the NH, in particular at higher altitudes. This is due to the acceleration of the Antarctic polar vortex (Fig. 9a) that leads to an increase in ∂q y (Fig. 9d), thereby reducing effective mixing. In most 10 models, ∆(K yy /|v * |) is negative within these latitudes throughout the stratosphere in the SH as well as in the NH. Between 20 and 50 hPa there are only two exceptions, namely NIWA-UKCA and ACCESS in the NH, those two models that also have a slight positive ∆ (see Tab. 2). Moreover, the two EMAC model versions, which have the largest ∆ also show the largest |∆K yy /|v * ||, in both hemispheres. Although we cannot find a clear inter-model correlation between ∆ and ∆(K yy /|v * |) here, these are strong indications of a link between the two quantities.
Hence, the inter-model differences in mixing efficiency changes are consistent with differences in the relative rate of change in mixing to residual transport. The differences in K yy /|v * | changes between the models appear to be related to structural differences in zonal wind changes, i.e. the models with an increasing mixing efficiency (namely NIWA-UKCA, ACCESS) 5 simulate zonal wind increases in the NH at higher latitudes than other models. In consequence, the PV gradient also increases at higher latitudes, so that mixing increases more strongly in the region of the subtropical barrier (see figures in Supplement).
Our approximations as well as the limited number of models do not warrant a robust conclusion from inter-model correlations.
However, the results shown here suggest that the mixing efficiency changes, as well as its inter-model spread are consistent with changes in the relative mixing strength due to changes in the background PV gradient.  ) in the analyzed model simulations, as well as on the causes for the large model differences in ∆AoA.
Linear separation of ∆AoA into changes by aging by mixing (A_mix) and changes by residual circulation Garny et al., 2014), shows that the contribution due to ∆A_mix dominates in all models. In particular, the influence of ∆A_mix controls almost the entire changes in AoA in the subtropical lower stratosphere. ∆RCTT is important in 20 the tropical pipe and in the downwelling branches of the extratropics, but only dominate in some models. This linear separation of ∆AoA in its components, however, is intricate in its interpretation. First, A_mix itself is a function of the residual circulation and therefore the individual processes are not independent from each other and second, this dependence between the processes is not necessarily local. By means of inter-model correlation analyses, we find that the changes in RCTTs and AoA are locally correlated only in the extratropical middle stratosphere. The connection between the changes in tropical upwelling (at 70 hPa) 25 and AoA, in contrast, are spread throughout the stratosphere, which points towards a non-local dependence of AoA from residual transport. The changes of the residual circulation as a consequence of tropical upwelling changes had been associated mainly with a subtropical jet stream acceleration due to a thermal wind response to upper tropospheric warming and lower stratospheric cooling previously (see e.g. Lorenz and DeWeaver, 2007;Randel et al., 2008;Butchart, 2014 We have shown that the mixing efficiency (Neu and Plumb, 1999;Garny et al., 2014) controls the strength of additional aging by mixing. Here, we reveal that the mixing efficiency decreases over time in most models, but two of them show a weak and one shows a strong positive ∆ . A decrease in mixing efficiency indicates that mixing does not increase as strongly as the residual circulation mass flux does. We find that, AoA increases more strongly in models with a stronger decrease in mixing efficiency, as less relative mixing reduces the fraction of air that recirculates around the BDC branches (Garny et al., 2014). Hence, stronger tropical upwelling as well as a reduced mixing efficiency both lead to a decrease of AoA. The temporal 5 changes in mixing efficiency contrast the results presented in Garny et al. (2014), who found a constant mixing efficiency in different climate states. However, as the changes in mixing efficiency appear to be model dependent, no change in mixing efficiency lies within the range of ∆ found here. Moreover, the analysis of PV gradient changes presented in Garny et al.
(2014) also suggests a decrease of relative mixing strength (see their Fig. 11).
Subsequently, we quantify the influence of mixing efficiency changes for ∆AoA as well as the impact of variations on To study the reasons for the changes in mixing efficiency, we analyze several dynamical fields as multi-model means of the differences between the two periods. These show the well-known coherence that the increase and upward shift of the zonal mean winds lead to enhanced wave dissipation in the lower stratosphere and thereby an amplified poleward residual transport. 20 However, changes in wave dissipation also lead to variations in the properties of mixing. A diffusivity coefficient that bases on the ratio of wave driving to the meridional potential vorticity gradient (under the premise of QG theory) reveals that the increase in wave dissipation shifts the region of strong mixing upward, thereby increasing the absolute strength of mixing in the lower stratosphere. Above, however, an enhanced PV gradient (which is due to the zonal wind increase) leads to a decrease in mixing strength. This counteraction of the two effects can explain why residual transport increases faster and that ∆ is negative in 25 most models. In the models with positive ∆ , the relative mixing strength is consistently found to increase, in particular in the NH, because the zonal wind changes, and thus the PV gradient changes, take place at higher latitudes. However, detailed process analysis experiments are required to test how robust this connection is.
In summary, we separated the effects of mixing changes from changes in residual circulation in causing the decrease in AoA.
We found that a decrease in the relative mixing strength leads to an AoA reduction of about 10% in the future in most models. 30 We further showed that the inter-model differences in simulated changes in mixing efficiency contribute to the inter-model spread in the simulated AoA changes. The decrease in relative mixing strength appears to be related to changes in background PV gradients. However, clear causalities can only be determined with model experiments that are specifically designed for that purpose, e.g. by varying certain parameters or model characteristics. The influence of mixing on the BDC and its changes can be important to project future climate conditions. We therefore suggest to conduct more in-depth analyses with the aim to study 35 the changes in residual circulation and mixing as well as their uncertainties and possible connections in more detail.
Author contributions. RE, SD and HG designed and conducted the analysis and wrote the paper. PŠ, TB and HB helped with discussions on content and structure of the study. In their role as CCMI model PIs, the other authors contributed information concerning their individual models and helped revise the article.

5
Data availability. All CCMI-1 data used in this study can be obtained through the British Atmospheric Data Centre (BADC) archive Moreover, we thank Andreas Engel for providing the in-situ AoA data. Andrews, D. G., Holton, J. R., and Leovy, C. B.: Middle atmosphere dynamics, vol. 40, Academic press, 1987.