The influence of mixing on stratospheric circulation 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 1 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2018-1110 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 25 October 2018 c © Author(s) 2018. CC BY 4.0 License.

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 aging 20 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 tropical 25 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 climatological AoA model 30 spread. In the present study, we quantify the impact of mixing efficiency (relative mixing strength) differences in 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.
2 Data and methods

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 via vertical motion and corresponds to transport by v * (Garny et al., 2014). can be derived by means of the tropical 10 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 wellmixed 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 w * T (z) can 15 thus be described as Here, H denotes the scale height (7 km), α the ratio of tropical to extratropical mass approximated by α = sin(lat(w * (z) = 0)) − sin(1 − lat(w * (z) = 0)) 2 · sin(90) − sin(lat(w * (z) = 0)) , 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 20 of the vertical residual velocity w * . 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 yields ( Eq. 3 shows that is approximately proportional to the relative increase in AoA due to mixing. Note that according to the 25 concept of the TLP model, has to be viewed as a parameter that counts 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), however, the authors nevertheless found a constant for three different climate states in one model. The tropical profiles provided for the TLP model are averaged over the latitudinal band of the models' individual vertically averaged turnaround latitudes (which are also timedependent) and are interpolated to vertical coordinates relative to the tropopause height of each model. The mixing efficiency is then obtained by a best fit 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 therefore a function of the mixing efficiency and of the residual circulation strength: A_mix is proportional to the mixing efficiency, but indirectly proportional to the vertical velocity. A higher mixing efficiency 5 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) 15 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 20 across the 21st century is investigated in detail by Šácha et al. (2018). They found out 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 these models are not part of our analysis. However, NIWA-UKCA and ACCESS are the direct successors to the UMUKCA models 25 and these range rather in the lower field 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 model was not included in the SPARC report. With higher vertical resolution (47 to 90 levels in the vertical), the EMAC model tends to simulate larger ∆AoA and thereby detaches from the bulk of the other model simulations.
In consistence with that, Revell et al. (2015) show that also in the SOCOLv3 model, AoA gets on average one year older 30 when the model is run with 90 layers in the vertical (instead of 39) because of less vertical diffusion. All the statements above also count for panel b) with the tropical to middle latitude AoA differences with altitude ( Fig. 1b). The spread is somewhat   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).
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 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. 5 The models may agree in predicting a decrease in stratospheric AoA, but they do show large differences in the strength of this trend. Also, the trends between the two periods (1960-2000 and 2000-2100)  analyses of mixing trends and in the second period, two mechanisms work against each 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. Tab. 2 provides an overview over mean AoA, and tropical upwelling (at 70 hPa) averaged between 1970 and 1990 (denoted as 1970) and between 2080 and 2100 (denoted as 2100) and their differences (∆) for the ten CCMI-1 REF- and over 10-100 hPa, which means that these AoA values are not comparable with those in Fig. 2. These values will be used to 10 explain or confirm some of the results throughout the paper.  Hence, values between 0.5 and 1 signify a domination of residual transport changes (reddish), while values between 0 and 0.5 mean that the ∆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).  (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 the vstar integral at 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 the inter-model correlations 25 between ∆AoA and ∆RCTT as well as between the ∆AoA and the differences in tropical upwelling in 70 hPa in Fig. 4. These correlations are calculated across the 10 model simulations for each grid point separately. Note that for this, the data of each model was interpolated to the resolution of the grid of the model with the highest horizontal resolution.

Coherences between the components
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  that the link between the residual circulation and AoA can not be expected to be local, but is rather of remote nature. It is mainly the ascent of air in the tropics that determines the influence that changes of the residual circulation have on ∆AoA.

5
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). Dietmüller et al. (2018) have found out 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 the residual circulation. Hence, we now lay our focus on the model trends (and differences between the two periods) in mixing 10 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).  In consequence, we next investigate if, apart from their connection to residual circulation changes, ∆AoA is also linked to 15 variations in mixing efficiency. 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 can not be reflected in the (sub)tropical measure . 5 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 also controls the transit times of the air parcels that recirculate. Overall, we can again conclude that ∆AoA in climate models is connected with both, the changes in mixing and in 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 10 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. Therefore, we further 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, 15 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  EMAC-L90 this 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 aging by mixing. 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 impact of the mixing efficiency on ∆AoA can be seen through AoA*. Or the other way round, the difference between 5 AoA*(2100) and AoA(2100) reflects the impact of ∆ on AoA. In the given example of the EMAC-L90 model, the fractional impact of the mixing efficiency change on the AoA difference is 29±4% calculated between 2.5 and 25 km above the tropopause. This value varies considerably among the models. Tab. 3 provides an overview of the quantity for the ten models.
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. But in the case of NIWA-UKCA, ACCESS and ULAQ, the mixing 10 efficiency increases (see Tab. 2), which means it leads to an increase in AoA, and not to a decrease as in the other models.
The negative ∆RCTT therefore accounts for more than the entire negative ∆AoA to compensate the effect of the change.
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 15 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 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 , α, T corr , and RCTT T . These quantities are taken to calculate tropical AoA by 20 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)    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 15 dynamical reasons for the changes in the 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 CCMI data are available). We present and analyze multi-model means (MMMs) 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 5 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 * |.
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.

10
The mechanism of the increase in residual circulation in the lower stratosphere is well understood (e.g. Shepherd and McLandress, 2011). It follows the increase in upper tropospheric temperatures that leads to an upward shift and increase in the zonal 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 in the MMM between about 100 15 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 * (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 20 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 means v and q, respectively. Note that due to 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., 25 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 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 Jr. et al., 1980;Cohen et al., 2014) with Similarly as in Abalos et al. (2017), we calculate the diffusivity coefficient K yy as This relation states that horizontal eddy mixing is proportional to the EP flux divergence, and indirectly proportional to the 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.

5
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 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 10 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 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 15 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 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 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 25 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 30 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).
However, as also 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 5 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 epsilon change also show the largest |∆K yy /|v * ||, in both hemispheres. Although we can not find a clear inter-model correlation between ∆ and in ∆(K yy /|v * |) here, these are strong indications of a link between the two quantities.

10
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) 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 in the region of the subtropical barrier, mixing increases more strongly (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.
4 Summary and conclusions 5 In the present study, we analyze the AoA differences of 10 CCMI-1 (Morgenstern et al., 2017) climate prediction simulations between the two periods 1970-1990 and 2080-2100. In agreement with previous model studies (Butchart et al., 2010;SPARC, 2010), we find a reduction of AoA over time in all model simulations. The smallest differences are consistently found in the tropics and the largest in the extratropical lower stratosphere, but the magnitude of the changes varies vastly among the models.
Our investigation focuses on the reasons for this negative ∆AoA (AoA(2080-2100)-AoA ) in the analyzed model 10 simulations, as well as on the causes for the large model differences in ∆AoA.
When we linearly separate ∆AoA into the differences that are driven by aging by mixing (A_mix) changes and by residual circulation transport time (RCTTs) changes Garny et al., 2014), we find that generally ∆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 the tropical pipe and in the downwelling branches of the extratropics, but only 15 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 (i.a.) 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) and AoA, in contrast, are spread throughout the stratosphere, 20 which also 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). Recently, this was linked with the upward shift of the tropopause e.g. by Oberländer-Hayn et al. (2016) and Abalos et al. (2017). The variations of mixing over time and their impact on stratospheric 25 circulation changes and thus AoA, however, are widely uncharted up to now.
In order to gain a deeper insight into the changes of mixing, we further investigate the mixing efficiency (Garny et al., 2014). This is a diagnostic quantity that is calculated by fitting the equations of a simplified model (the TLP model; Neu and Plumb, 1999) to the GCM data and provides information about the ratio of the mixing mass flux to the net mass flux. The mixing efficiency controls the strength of additional aging by mixing. Here, we reveal that the mixing efficiency decreases 30 over time in most of the model simulations, but two models show a weak and one shows a strong positive ∆ . A decrease in the mixing efficiency indicates that mixing does not increase as strongly as the residual circulation mass flux does. We find that in models with a stronger decrease in mixing efficiency, AoA increases more strongly, 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 changes in mixing efficiency over time found here in the CCMI simulations 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 appears to be model dependent, a zero change in mixing efficiency lies within the range of ∆ found here. Moreover, note that the analysis of PV gradient changes presented in Garny et al. (2014) 5 also suggest a decrease of the relative mixing strength (see their Fig. 11).
Subsequently, we quantify the influence of changes in mixing efficiency for ∆AoA as well as the impact of variations on the ∆AoA model spread for all model simulations, individually. The influence of ∆ on ∆AoA is determined by means of calculating AoA for 2100 but with the of 1970 (AoA * ) with the TLP model and then comparing the hypothetical ∆AoA * with ∆AoA. We obtain a multi-model mean of 10.4% of the influence of mixing efficiency changes on the differences in AoA.

10
However, the models show a large spread in this quantity (σ=18.3%) and as three models possess a positive ∆ , also the sign is not consistent among the models. In a similar manner, we assess the impact of mixing efficiency variations on the model spread in ∆AoA. This reveals that first, the AoA changes are generally smaller when is kept constant (at the 1970-1990 mean values) and second, that some of the model spread is caused by variations in . This means that model differences in changes lead to a considerable enhancement of the model inconsistencies in future projections of the BDC. 15 To study the reasons for the changes in the 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.
However, changes in wave dissipation also lead to variations in the properties of mixing. To see if this is reflected in the model data, we calculate a diffusivity coefficient that bases on the ratio of wave driving to the meridional potential vorticity gradient 20 (under the premise of QG theory). This 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 stronger than mixing and thus the negative ∆ in most of the models. In the models with an increase in , the relative mixing strength is consistently found to increase, in particular in the NH, because the zonal wind 25 changes, and thus the PV gradient changes, take place at higher latitudes in these models. However, detailed process analysis experiments are required to test how robust this connection is.
Overall 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 a further reduction of about 10% in AoA in the future in most models. We further could show that the inter-model differences in simulated changes in the mixing efficiency contribute to the 30 inter-model spread in the simulated AoA changes. The decrease in the 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 anthropogenic and biomass burning emissions of air pollutants at global and regional scales during the 1980-2010 period, Climate Change,109,doi:10.1007/s10584-011-0154-1, 2011.