Interactive comment on “ Time varying changes in the simulated structure of the Brewer Dobson Circulation ”

This paper presents GEOSCCM model simulations of the Brewer-Dobson circulation (BDC) and mean age over the past decades. The authors consider different ensemble members from their simulation which show a very different evolution of the BDC. One member even shows increasing mean age in the NH since 1988, similar to observed mean age trends. In the lower stratosphere, on the contrary, the BDC continues to accelerate. Hence, structural changes evolve in the BDC after 1988. ODS and volcanic eruptions are identified as the main forcing agents of these structural changes.


Introduction
The global circulation in the stratosphere -the Brewer-Dobson circulation (BDC) -consists of air masses rising across the tropical tropopause, moving poleward, and sinking into the extratropical troposphere (Holton et al., 1995;Waugh and Hall, 2002;Butchart, 2014).Since the BDC and its changes have important implications on both stratospheric and tropospheric climate as well as stratospheric ozone chemistry (SPARC-CCMVal, 2010;World Meteorological Organization, 2011, 2014;Manzini et al., 2014), it is important to assess the factors that lead to simulated BDC changes and whether historical changes in the BDC as simulated by models are inconsistent with available observational constraints.
The BDC has historically been deduced either from the residual circulation or from the average time for an air parcel to travel from the tropical troposphere to a given stratospheric sampling region (i.e., the mean age of air or mean age).While these two diagnostics are clearly related, differences can arise due to isentropic mixing and recirculation (Waugh and Hall, 2002;Strahan et al., 2009;Li et al., 2012).Specifically, the vertical component of the residual circulation (w * ) measures the instantaneous advection, whereas mean age is an integrated measure of the total transport.It would not be a surprise if the two metrics had a different evolution.
Chemistry-climate models robustly predict a strengthened BDC under climate change in the middle and lower stratosphere of approximately 1-5 % per decade (the precise rate depends on the level considered and varies among models; C. I. Garfinkel et al.: Changes in BDC structure Butchart and Scaife, 2001;Butchart et al., 2006;Garcia and Randel, 2008;Li et al., 2008;Waugh, 2009;Shepherd and McLandress, 2011;Garcia et al., 2011;Lin and Fu, 2013;Butchart, 2014;Oberländer-Hayn et al., 2015).The model used in this study, the Goddard Earth Observing System Chemistry Climate Model, Version 2 (GEOSCCM), predicts a trend quantitatively similar to those in other models both for the historical period and for the future (Oman et al., 2009;Waugh, 2009;Butchart et al., 2010;Li et al., 2012).
It has been argued that observational estimates of historical changes do not agree with the simulated acceleration trend.Specifically, the analysis of historical tracer data does not provide evidence for an acceleration trend in the midstratosphere Northern Hemisphere (NH), where mean age actually appears to have increased or remained unchanged (Engel et al., 2009;Bönisch et al., 2011;Stiller et al., 2012;Hegglin et al., 2014;Ray et al., 2014).In particular, the mean age evolution in the figures of Engel et al. (2009) and Ray et al. (2014) indicates aging since the late 1980s, with earlier changes less clear.ERA-Interim reanalysis data also suggest aging of NH mid-stratosphere air since the late 1980s (Diallo et al., 2012;Monge-Sanz et al., 2013;Ploeger et al., 2015).While these observational and reanalysis-based studies disagree about the sign of changes in other regions of the stratosphere, they all indicate aging of air in the NH midlatitude mid-stratosphere.Ray et al. (2014) argue that the large uncertainty estimates on the trends presented by Engel et al. (2009) are overly conservative and that this aging trend is statistically significant even after taking into account the nonlinear growth rates of these trace gases (i.e., the arguments of Garcia et al., 2011).Several of the aforementioned studies suggest that these changes in the mean age imply a redistribution of the BDC, and specifically a slowdown of the deep (i.e., mid-stratospheric) NH branch of the BDC and/or less mixing of fresh tropical air into this region.However, it is not clear what forcings (if any) could be responsible for this redistribution and also whether models of the kind used in WMO ozone assessments (World Meteorological Organization, 2011, 2014) can capture such a slowdown given known forcings.Oman et al. (2009) conducted time-slice and transient simulations using a previous version of the model used in this study, and they found a trend towards younger mean age everywhere in the stratosphere between 1960 and 2004, though the decrease in mean age in the Southern Hemisphere (SH) is larger due to Antarctic ozone depletion.Oman et al. (2009) also found that ozone recovery would lead to a slowdown of the BDC if not for warming sea surface temperatures (SSTs) due to increasing greenhouse gas (GHG) emissions.Oberländer-Hayn et al. (2015) recently presented differences in mean age and the residual circulation in time-slice simulations using the ECHAM/MESSy Atmospheric Chemistry Model (EMAC).The changes in tropical upward mass flux indicate a strengthening of the BDC between 1960 and 2000 in the NH winter season in the lower stratosphere and a weakening in the upper stratosphere with a change in sign at 10 hPa.Changes in mean age show a decrease of about 0.13 year decade −1 in the lower and middle stratosphere and a slight increase in the Arctic upper stratosphere and lower mesosphere.While there is some hint of a structural change in the properties of the BDC, the changes occur higher in the stratosphere and at more polar latitudes than is suggested by available observations.
The motivation for studying historical changes in the BDC in free-running climate simulations is not to form a best estimate of the actual historical evolution; for that purpose, reanalyses and/or nudged experiments are far better (though care must be exercised when interpreting trends).Rather, the motivation is twofold: one, future projections of the BDC can only be produced by free-running climate simulations, and these projections are of limited value if a model's simulation of the past is inconsistent with observational constraints; two, assuming the model is capable of following the observed evolution, the forcings that caused these historical changes can be systematically diagnosed by sequentially adding these forcings.Combined, the goal of this study is to understand whether historical forcings could have led to the structural changes in the BDC that have been inferred by observational studies.
We show that over the full duration of the experiments (i.e., for a start date in 1960) we recover the result from previous modeling studies: anthropogenic climate change leads to acceleration of the BDC throughout the stratosphere.However, our model can simulate statistically significant aging of air in the midlatitude NH near 20hPa between the early 1990s and the present as suggested by available observations.This suggests that historical forcings have caused structural changes of the BDC since the late 1980s, whereby the BDC accelerated in the lower stratosphere but decelerated in the mid-stratosphere, in both the tropics and in the NH.Mean age and the residual circulation (as measured by tropical w * ) change in unison.The cause of this mid-stratospheric deceleration trend is a combination of forcings -ozone-depleting substance (ODS) recovery and the timing of volcanic eruptions -that together has outweighed greenhouse-gas-induced acceleration since the late 1980s.We therefore emphasize that, if one wishes to capture observed historical changes in a model simulation, careful attention must be paid to the forcings included and the start and end dates used for trend calculation.

Methods
The model focused on for this study is GEOSCCM, an aerosol-and chemistry-focused version of the GEOS-5 Earth system model.GEOSCCM couples the GEOS-5 (Rienecker et al., 2008;Molod et al., 2012) atmospheric general circulation model to the comprehensive stratospheric chemistry module StratChem described in Pawson et al. (2008) and the Goddard Chemistry, Aerosol, Radiation, and Transport Model (GOCART) described in Colarco et al. (2010).The model has 72 vertical layers, with a model top at 0.01 hPa, and all simulations discussed here were performed at 2 • latitude × 2.5 • longitude horizontal resolution.Previous versions of GEOSCCM have been graded highly in the two phases of the Chemistry-Climate Model Validation Activity (Eyring et al., 2006;SPARC-CCMVal, 2010).Improvements to the model since then are described in Oman and Douglass (2014) and Aquila et al. (2016).Note that this model version includes a spontaneous quasi-biennial oscillation (QBO) and that the QBO phase differs among the three ensemble members.
A series of simulations of the period from January 1960 to December 2014 have been performed in order to understand the past evolution of the stratosphere.These simulations were presented in Aquila et al. (2016), where the focus was on changes in temperatures.Here we examine changes in the BDC.We start with an ensemble in which the only time-varying forcing is changing SSTs and sea ice, and then sequentially add the following forcings: greenhouse gases, ozone-depleting substances, volcanic eruptions, and solar variability.More specifically, these simulations are grouped into the following five experiments: 1. SST, which uses time-varying observed SSTs and sea ice up to November 2006 from the Met Office Hadley Centre observational dataset (Rayner et al., 2006) and since then from Reynolds et al. (2002) and subsequent updates (Fig. 1a).GHGs and ODS concentrations are fixed at 1960 levels.Volcanic eruptions are not included in this experiment, and the solar forcing is held constant.
All simulations used emissions of tropospheric aerosol and aerosol precursors following Granier et al. (2011).We focus on transient simulations as it is difficult to compare historical tracer observations to time-slice simulations.
Each experiment is composed of three ensemble members initialized with different initial conditions from a 1960 timeslice simulation.Because we have three members for each forcing combination, we can also assess at least partially the range of internal atmospheric variability.This internal variability is not a source of model uncertainty; rather it is an inherent part of the climate system.If the BDC in one ensemble member, but not in the other two, evolves consistently with observational constraints, one can reasonably conclude that models can capture the observed trends if part of the observed trend was due to internal variability and was not forced (cf.Santer et al., 2008).
We assume that BDC perturbations induced by each forcing agent add linearly to the others, as previous work focusing on forcing agents for lower-stratospheric mass flux and mid-stratospheric mean age suggests that nonlinearities are small (e.g., Oman et al., 2009;McLandress et al., 2010).The model version used to perform the integrations described in this paper is no longer supported, and hence we cannot explicitly test this assumption.
The use of observed SSTs in our simulations, rather than internally calculated by the model, produces a climate state closer to the observed one.However, partitioning trends into an SST-driven component and a component from other radiative or chemical forcings is somewhat artificial, as the prescribed SST changes occur in response to and in tandem with the changing direct atmospheric forcing; however such a partitioning is an effective tool for disentangling the physical mechanisms leading to changes in the atmospheric circulation (Deser and Phillips, 2009).Specifically, Oman et al. (2009) found that the mean tropospheric warming associated with rising SSTs had a bigger impact on mean age than the direct radiative impacts of CO 2 .Note that interannual and decadal variability in SSTs drives changes in the BDC that likely have nothing to do with climate change, and hence we include a smoothed version of the SST variations in Fig. 1a (all experiments were forced with the full time-evolving SST fields).
A passive tracer is used to derive the mean age.The mixing ratio of this tracer increases linearly with time, and the time lag in tracer concentrations between a certain grid point in the stratosphere and a reference point in the troposphere provides an estimation of mean age at this stratospheric grid point (Waugh and Hall, 2002).We adopt a reference point of 200 hPa at the Equator.Mean age was first initialized in a 10-year time-slice control run with 1960 conditions before we branched off for each experiment, and so mean age can be defined at the beginning of the experiments.Note that diagnostic output necessary to compute the full age spectrum was not saved for these model experiments, and hence we are limited in our ability to quantify mixing changes.The aging of NH mid-stratospheric air in observations is pronounced mainly after the late 1980s (see figures in Engel et al., 2009;Ray et al., 2014), and several of the reanalysisbased studies begin their analysis in the late 1980s as they need to initialize trajectories for 10 years before computing mean age.Hence, in our presentation and discussion of the results, we consider trends and variability both over the full period of integration and also since 1988.Motivated by our results, we also discuss trends for a start date of 1992, as the recent slowdown of the deep branch of the BDC is most pronounced (and locally statistically significant) for this start date due to the eruption of Mt.Pinatubo.Finally, global mean age profiles as retrieved by satellites are only available from 2002 onward (Stiller et al., 2012;Haenel et al., 2015), and hence we show trends since 2002 separately as well.The trends are calculated with a linear least-squares fit.Statistical significance of the trends in individual ensemble members of GEOSCCM is computed using a two-tailed Student's t test, and the reduction in degrees of freedom due to autocorrelation of the residuals is taken into account with the formula N (1−r 1 )(1+r 1 ) −1 , where N is the number of years and r 1 is the lag-1 autocorrelation (Eq.6 of Santer et al., 2008).In computing the ensemble mean response, we first averaged the mean age among the three ensemble members for each season/year (in order to damp the internal, unforced, atmospheric variability) and then computed the trend based on the ensemble mean mean age.

Results
We now examine how the BDC as simulated by GEOSCCM has changed in structure over the past 55 years.We show that BDC changes (or "trends") vary with period considered and location.These structural changes are associated with several distinct forcings, and these forcings transiently drive changes in the BDC.When combined with internal variability, it is possible that these forcings can drive, over the 27year period of 1988-2014, a deceleration trend in the NH mid-stratosphere.

Changes in the all-forcing ensemble
We begin with changes in the BDC for the all-forcing ensemble.We first consider changes in the residual circulation and in mean age as a function of latitude and pressure over the full duration of the experiment in the ensemble mean (Fig. 2).The BDC accelerates throughout the stratosphere as mean age decreases and the residual circulation accelerates.Hence, changes over the full integration period in our experiments are consistent with previous work (e.g., Oman et al., 2009).Statistical significance of the trends are considered in Figs. 3 and 4, which show the 95 % confidence bounds on the trends for a range of start dates for mean age in the selected regions and for tropical mean upwelling.The changes in both metrics of the BDC are statistically significant at the 95 % level throughout the stratosphere for a trend start date of 1960, and all three ensemble members indicate quantitatively similar acceleration trends.
However, the ensemble mean acceleration trend weakens (and even reverses locally) and its robustness goes away as we consider more recent periods.To demonstrate this, we start by showing trends since 1988 in Fig. 5.In the lower stratosphere (i.e., the shallow branch of the BDC), the BDC continues to accelerate in all three ensemble members, and this change is statistically significant at the 95 % level in each ensemble member individually and in the ensemble mean (e.g., Figs.3d, e, f and 4b).In the mid-and upper stratosphere, however, trends are not robust across the various ensemble members.One of the three ensemble members shows decreasing mean age and an accelerated residual circulation (Fig. 5a and gray lines in Figs.3a, b, c and 4a), while another shows the opposite (Fig. 5e and red lines in Figs.3a,  b, c and 4a).In this ensemble member with aging air, upwelling decreases throughout the tropics, such that both the residual circulation and mean age diagnostics suggest deceleration of the BDC.None of the mean age trends in the mid-and upper stratosphere in Fig. 5a, c, and e are statistically significant at the 95 % level (e.g., Fig. 3a, b, c).Note that if the start date for the trend is advanced to 1992 then none of the three ensemble members indicates acceleration of the BDC in the NH midlatitude mid-stratosphere (Fig. 3a), but we prefer to demonstrate that, even for a start date of 1988, aging can be simulated given the large amount of internal atmospheric variability.
These differences in trends since 1988 among ensemble members can be reconciled with the changes in the wave forcing of the BDC.Previous modeling and theoretical work have demonstrated that changes in the wave forcing directly force changes in the residual circulation (Rosenlof, 1995;Butchart, 2014).Specifically, enhanced wave convergence (i.e., deceleration of the mean flow) leads to enhanced upwelling on the tropical flank of the enhanced wave convergence and downwelling on the poleward flank (Haynes et al., 1991).As both resolved waves and gravity waves are important for the total wave driving, we evaluate their combined impact in common units of acceleration of the mean flow (m s −1 day −1 decade −1 ); positive values indicate less wave convergence and deceleration of the residual circulation, while negative values indicates enhanced wave convergence and acceleration of the residual circulation.The right column of Fig. 5 demonstrates that the difference between the ensemble member with a weakened BDC and the one with an accelerated BDC is related to differences in the SH wave driving and subsequent changes in the SH residual circulation.For the ensemble member with a weakened residual circulation in the SH (and thus aging of air in the mid- stratosphere), there is less wave flux converging in the SH (Fig. 5f).In contrast, for the member with an accelerated residual circulation in the SH (and thus decreasing mean age in the NH mid-stratosphere), there is enhanced wave flux converging in the SH (Fig. 5b).Note that a difference in wave fluxes in the SH can influence mean age in the NH because mean age is an integral measure of transport, and thus changes in the residual circulation in the tropics due to wave flux changes in the SH can impact the transport pathway into the NH.Hence, the difference between aging and freshening of mid-stratospheric air is associated with the internal atmospheric variability associated with wave fluxes.The Supplement discusses implications of this intra-ensemble difference for trends in temperature and zonal wind.There is reduced wave convergence in the NH upper stratosphere in all three ensemble members, and this effect is most pronounced in boreal winter (not shown).The cause of the decrease in NH stratospheric wave convergence is discussed in Garfinkel et al. (2015) and in the next section.Overall, our model simulations indicate that known forcings could have led to a slowdown of the deep branch of BDC since 1988 given the large amount of internal variability in wave fluxes.
We now consider the time evolution of changes in mean age in various regions in Fig. 6.As the eruption of Mt.Pinatubo had a pronounced impact on the BDC in these simulations, we separately consider the evolution before and after its eruption.Changes at NH midlatitudes in the midstratosphere are considered in Fig. 6a.Note that this is the region where available observations suggest that air has aged over the past several decades.In the all-forcing ensemble mean (blue line), mean age decreases by 0.45 year between 1960 and 1992 (i.e., 0.14 year per decade, blue in Fig. 7a) but then ages by 0.12 year after 1992 (i.e., 0.05 year per decade, blue in Fig. 7b).The aging trend since 1992 is statistically significant at the 95 % level in two of the three ensemble members and in the ensemble mean (Fig. 3a). Figure 8 considers the evolution of each of the three ensemble members individually.The ensemble spread in mean age in any given year can exceed 5 %.One ensemble member simulates anomalously old air after 2010, and for this member aging trends are simulated for a start date of the trend calculation of either 1992 or 1988.Note that a different ensemble member simulates anomalously younger air relative to the other two after El Chichón, highlighting the large role of internal atmospheric variability in masking the response to climate forcings.
In the tropical and SH mid-stratosphere, mean age also increases from 1992 through the end of the simulation (blue curves in Figs.6b, c and 3b, c).The large uncertainty in the trends in the SH arise due to the correction of Santer et al. (2008) for the effective number of degrees of freedom associated with residuals with large autocorrelation.In contrast, in the lower stratosphere, mean age continues to decline in the tropics and in the NH, though not in the SH (blue curves in Fig. 6d, e, f).The decline in mean age is statistically significant at the 95 % level in two of three ensemble members and in the ensemble mean in the NH lower stratosphere (Fig. 3d).Hence, there are time-and space-varying variations in recent mean age trends, and only in some regions of the stratosphere has mean age continued to decline.
Similar structural changes are evident for the tropical upwelling mass flux between the turnaround latitudes ( w * =0 SH ρw * dA, where ρ is the density and A is the area of a given zonal band; Figs. 4 and 9).In the all-forcing experiment, tropical upwelling accelerated until 1992 in both the mid-stratosphere and lower stratosphere (i.e., blue line rises in Fig. 9, blue in Fig. 10a, c, e, g), but since 1992 it has decreased at all levels above 70 hPa (e.g., Figs.9a and 10b, d, f,  h).These changes are reminiscent of those proposed by Ray et al. (2010Ray et al. ( , 2014) ) in order to explain how NH midlatitude air in the mid-stratosphere anomalously ages.Specifically, Ray et al. (2010) concluded that for a tropical leaky-pipe model, "the best quantitative agreement with the observed mean age and ozone trends over the past three decades is found assuming a small strengthening of the mean circulation in the lower stratosphere, [and] a moderate weakening of the mean circulation in the middle and upper stratosphere", as simulated by GEOSCCM.Overall, it is clear that GEOSCCM can simulate structural changes in the BDC that resemble those inferred from observations.

Forcing of the trends
We now consider the forcing mechanisms behind these structural changes in the BDC.

Mid-stratosphere
As shown above, mean age at NH midlatitudes in the midstratosphere decreases by 0.45 year between 1960 and 1992 in the all-forcing ensemble but then ages by 0.12 year after 1992 (blue line; Figs.6a and 7).This evolution can be broken down into its various forcing components.
1. SSTs.The red line in Fig. 6a shows that SSTs lead to a decrease in mean age of 0.1 year over the course of these 55 years, albeit with substantial interannual and decadal variability.Specifically, over the last 30 years of the integrations (from 1985 to the end), there is a weak and insignificant aging trend (e.g., red in Fig. 7b).The likely cause of this is a tug-of-war of opposing effects.On the one hand, gradual warming of the oceans in isolation leads to an acceleration of the BDC (Fig. 1a and Oman et al., 2009).On the other hand, the spatial pattern of recent changes in SSTs has led to a decline in planetary wave flux (especially wave 1) entering the NH stratosphere at midlatitudes (Garfinkel et al., 2015): the vertical component of the Eliassen-Palm flux at 100 hPa area averaged between 40 and 80 • N declines in all three all forcing and in all three SST-only experiments, and the decrease is statistically significant at the 95 % level in the ensemble mean both in January through March (the focus of Garfinkel et al., 2015) and in the annual mean.This decline in upward-propagating midlatitude planetary waves at 100 hPa impacts the deep branch more strongly (Plumb, 2002;Ueyama et al., 2013).Hence, it is not surprising that little change has occurred over the last 30 years.2. Greenhouse gases.The influence of GHG can be deduced from the difference between the red and magenta curves in Fig. 6a, as the SST-only ensemble is conducted with radiative forcings fixed at 1960 levels.The difference between the curves exceeds 0.1 year towards the end of the integrations.
3. Declining ODS concentrations.The effect of increasing ODS concentrations can be deduced from the difference between the magenta and green curves in Fig. 6a, and it suggests that increasing ODS concentrations led to a decrease in mean age of 0.25 year by the late 1990s (Fig. 6a) when the ODS burden was peaking.More recently, declining ODS concentrations (Fig. 1c) lead to a recovery towards older air; that is, the gap between the magenta and green curves decreases between the late 1990s and the present in caveat that, while declining ODS concentrations clearly impact the BDC in these integrations, a statistically significant recovery of ozone has been detected in observations only in the upper stratosphere at midlatitudes and in the tropics (World Meteorological Organization, 2014).
3. Volcanoes and solar.The influence of volcanic eruptions can be deduced from the difference between the green and cyan curves in Fig. 6a, and in our model simulations the eruption of Mt.Pinatubo and El Chichón led to a decrease in mean age of 0.2 year, which gradually decayed over 4 to 6 years.Minor volcanic eruptions in the past 10 years may have led to an additional decrease in mean age of 0.05 to 0.1 year.The net effect is that large eruptions (or lack thereof) can influence decadal variability in mean age in GEOSCCM.Solar influences appear to be relatively minor (Fig. 6a), and we therefore focus our attention on the other forcings in this paper.
The net effect is that, over the second half of the experiments (when observations are more numerous), decreasing ODS concentrations and the recovery from Pinatubo overcame the influence of rising GHG concentrations and led to aging of 0.12 year.The same forcings that led to a statistically significant aging trend in the NH after 1992 in the ensemble mean and in two integrations also led to similar significant aging in the tropics (blue curve in Figs.3b and 6b).Specifically, changes in this region are dominated by ODS concentrations, and since ODS concentrations decrease after the late 1990s (Fig. 1c), mean age increases despite rising GHG concentrations.As for the NH, the earliest start date of significant aging trends is in the early 1990s due to the eruption of Mt.Pinatubo.ODS concentrations also dominate the SH midstratospheric mean age evolution (Fig. 6c), and since ODS concentrations decrease after the late 1990s, mean age does not change significantly after 1992 in the all-forcing ensemble as well.Note that the influence of the eruption of Mt.Pinatubo is weaker in the SH than in the NH; Mt.Pinatubo is located at 15 • N, and in our experiments the majority of the aerosols stay in the NH.

Lower stratosphere
In the tropical and NH lower stratosphere (Fig. 6d-e), increasing GHG concentrations and warming SSTs drive a decrease in mean age throughout the period of the all-forcing ensemble integrations.In agreement with the modeling results of Lin et al. (2015), warming SSTs (cf.Fig. 1a) impact the shallow branch more strongly than the deep branch.As midlatitude planetary waves are less important for the shallow branch than for the deep branch (Plumb, 2002;Ueyama et al., 2013;Abalos et al., 2014), it is reasonable to expect that a reduction in their strength has a relatively smaller impact on the shallow branch.Volcanic eruptions have a weaker effect in the lower stratosphere as compared to the middle stratosphere (compare the difference between the blue and green curves for the years 1963/1964, 1983/1984, and 1991/1992 between Fig. 6a and d).A possible explanation is that longwave and near-infrared heating due to volcanic aerosols occurs at the level of the aerosols (not shown), and in our integrations the aerosols are quickly lofted higher in the stratosphere (compare the volcanic influence on temperature as a function of time in the various levels of Aquila et al., 2016).This allows for stronger and longer-lasting changes mainly in the "deep" branch of the BDC.As for volcanoes, changing ODS concentrations impact the deep branch of the BDC more strongly (except in the SH where ozone depletion is strongest), and this effect is consistent with the idealized modeling results of Gerber (2012).The colder vortex that follows ozone depletion creates a waveguide higher into the stratosphere, raising the breaking level of Rossby waves and deepening the BDC.Hence, it is to be expected that ozone depletion and recovery have a disproportionate impact on the deep branch.Overall, the aging in the deep branch since 1992 does not extend to the shallow branch because the two factors that led to aging (recovery from the Pinatubo eruption and declining ODS concentrations) preferentially impact the deep branch, while the two factors that led to freshening (GHG increases and SST warming) preferentially impact the shallow branch.
In the SH lower stratosphere (Fig. 6f), ODS concentrations are the dominant forcing, but the gradual decline in ODS concentrations is balanced out by rising GHG concentrations and mean age is flat after 1995 in the blue all-forcing curve.It is known that ozone-depletion-induced polar cooling can directly modulate extratropical wave propagation down to the troposphere (Chen and Held, 2007;Oman et al., 2009;Garfinkel et al., 2013), though future work is needed in order to understand how this influence led to an accelerated shallow branch.

Residual circulation
The same forcings led to structural changes in tropical upwelling (Figs. 4 and 9).In the all-forcing experiment, tropical upwelling accelerated until 1992 in both the mid-stratosphere and lower stratosphere (i.e., blue line rises in both panels of Fig. 9), and this acceleration is driven largely by rising greenhouse gas concentrations, warming SSTs, and ozone depletion (Fig. 10a, c, e, g).Since 1992, however, the residual circulation has decelerated at 50 and 70 hPa (Figs. 4a,9a,and 10b,d).The recent deceleration of w * at 50 hPa comes about due to competition between changing SSTs and declining ODS concentrations: changing SSTs lead to continual though insignificant acceleration (i.e., the red curve continues to rise in Fig. 9a though trends in Fig. 10b, d are not significant), but ODS recovery leads to a slight deceleration (i.e., the gap between the red curve and green curve gradually decreases after 2000 in Fig. 9a).Because of the eruption of Pinatubo, the deceleration trend starts in 1992 in the all-forcing experiment (blue curve in Fig. 9a) rather than in the late 1990s, when ODS concentrations began to decrease (Fig. 1a).For a start date of 1991, the deceleration trend at 50hPa is statistically significant at the 95 % level in two ensemble members and in the ensemble mean (Fig. 4a).At 100 hPa (Fig. 9b), on the other hand, the dominant forcing is SSTs, and w * continues to increase throughout the experiment (Fig. 4b).The implications of these changes for ozone and temperature in the lower stratosphere are discussed in Polvani et al. (2017).Overall, the same forcings that control the age of mid-stratospheric air also control the residual circulation, and these forcings can explain the lack of acceleration of the deep branch of the BDC since 1992.

Summary of key forcing agents
In summary, from 1960 to the late 1980s (the first half of the experiments), ozone depletion, rising GHG, and warming SSTs all led to a decrease in mean age in all regions of the stratosphere.Over the second half of the experiments (since the early 1990s), rising GHGs have continued to lead to decreasing mean age, though the decrease is more prominent in the lower stratosphere.However, declining ODS concentrations and the proximity of the start date to the eruption of Pinatubo lead to an aging trend that is most prominent in the mid-stratosphere.The degree of compensation between these forcings is region-specific, and for the NH midlatitude midstratosphere the volcanic effects and declining ODS concentrations dominate, while in the lower stratosphere the SSTs and GHGs dominate.Hence, structural changes occurred in the BDC in our simulations.

Discussion of observed changes
There are no direct measurements of historical changes in the BDC.However, its past evolution can be deduced from trace gas measurements or from satellite data, and here we consider whether the modeled evolution of the BDC in GEOSCCM is consistent with these constraints.

Comparison with BDC changes inferred from in
situ CO 2 and SF 6 concentrations since 1975 Balloon measurements of CO 2 and SF 6 concentrations are available from 1975, and these data do not provide evidence for an acceleration trend in the mid-stratosphere NH, where mean age actually appears to have increased (Engel et al., 2009;Ray et al., 2014).In particular, the mean age evolution in the figures of Engel et al. (2009) and Ray et al. (2014) indicates aging since the late 1980s, with earlier changes less clear.As discussed in Sect.3.1, a similar evolution is present in our simulations.In order to make the comparison more precise, we subsample the simulated mean age on the day of each flight analyzed by Ray et al. (2014).As daily threedimensional fields of mean age were not archived from the simulations, we cannot map the simulated age into equivalent latitude space (as done by Ray et al., 2014).We show all three GEOSCCM members in order to estimate the internal variability in the model-simulated mean age (see Fig. 11).GEOSCCM captures the climatological value of mean age averaged over this period accurately: the difference between the observations and model for these data points is 3 months.A similar 3-month offset is evident when comparing GEOSCCM to the mean ages reported by Engel et al. (2009), which falls within the 6-month uncertainty in the observations (Engel et al., 2009).The climatological value of mean age in other regions also agrees well with satellite-based estimates presented in Stiller et al. (2008).
GEOSCCM mean age lies within the error bar for most measurements.While the weak (non-significant) aging trend noted in observations since 1975 is not present in GEOSCCM, observed and modeled trends agree within the 95 % uncertainty level.Note that, if we use the wider uncertainties reported by Engel et al. (2009), mean age in GEOSCCM agrees with all balloon flights and trends agree within the 90 % uncertainty level (Fig. 12).Over the recent period (since 1992), one of the three ensemble members simulates a trend in close agreement with the observed trend (0.12 ± 0.1 year decade −1 for GEOSCCM and 0.14±0.14year decade −1 for observations).That being said, there is apparently less subseasonal and QBO variability in GEOSCCM than in the observations (and also the tropical leaky-pipe model of Ray et al., 2014), and the trend towards older air is weaker in the GEOSCCM ensemble mean than in the observations.(Changes in mixing cannot be diagnosed from these simulations, but it is conceivable that mixing could account for some of the difference).For future work, we will consider whether GEOSCCM is consistent with more recent tracer measurements.

Comparison with BDC changes since 2002
While extreme caution must be exercised in interpreting a trend over such a short period due to the large stochastic variability in the atmospheric circulation, we now assess whether BDC changes since 2002 in GEOSCCM are consistent with observational constraints.
Vertically and latitudinally resolved changes in satellitemeasured SF 6 are available from 2002 onward, and Haenel et al. (2015) infer mean age trends from these data (their Fig. 6).They find that mean age declines in the tropical lower and mid-stratosphere south of the Equator and increases at NH midlatitudes and in the SH polar stratosphere.We show changes in annual averaged mean age from January 2002 to December 2011 in Fig. 13.The model simulates younger mean age in the lower stratosphere in all three ensemble members, but changes higher in the stratosphere are not robust among the various ensemble members.These intraensemble differences highlight the fact that one should not base any conclusions on the long-term behavior of the BDC on 10-year trends, as trends over 1 decade are strongly influenced by unforced (internal) variability.While not one of the ensemble members captures the interhemispheric dipole in the trends as suggested by satellite data (though individual integrations separately capture half of the dipole), we suggest that such a correspondence should not necessarily be expected as the wave forcing of the BDC differs in any realization of the atmospheric state (e.g., see the discussion in Santer et al., 2008).Aschmann et al. (2014) deduce changes in the tropical mass upwelling from changes in ozone, and they infer a lack of acceleration in mass upwelling above 70 hPa since 2002 (no significant changes) and an acceleration below 70 hPa of 11 % per decade.Quantitatively similar behavior is evident in Fig. 9 -upwelling increases at 100 hPa by approximately 10 % between 2000 and 2014 but decreases at 50 hPa.Hence, the trend towards younger air in the NH lower stratosphere in GEOSCCM is consistent with observational constraints.Aschmann et al. (2014) further speculate that this slowdown  2009), and the uncertainty of the model simulated mean age can be deduced from the intra-ensemble spread.The GEOSCCM mean age is offset by 3 months, i.e., the bias in the mean age (which is less than the 6month uncertainty in the observed mean age as evaluated by Engel et al., 2009).
of the upwelling above 70hPa is associated with the La Niñalike sea surface temperature trends over this period.This conjecture is supported by our modeling results: in the SST-only experiment, w * at 100 hPa continues to increase over this period but is largely flat at 50 hPa (and also at 70 and 30 hPa, not shown).Garfinkel et al. (2015) also noted that recent changes in SSTs (including the La Niña-like sea surface temperature trends) lead to less planetary wave flux entering the stratosphere at midlatitudes, and changes in midlatitude planetary waves will impact the deep branch more strongly.

Response to the eruption of Mt. Pinatubo
The eruption of Mt.Pinatubo led to younger mean age throughout the stratosphere and enhanced tropical upwelling in our model.Figure 2 of Garcia et al. (2011) suggests that similar behavior is present in the Whole Atmosphere Community Climate Model (WACCM).Similar behavior is evident in the mid-stratosphere in SOCOL (SOlar Climate Ozone Links), though not in the lower stratosphere (Muthers et al., 2016).Diallo et al. (2012) also infer older mean age in the lower stratosphere following Pinatubo using ERA-Interim data.However, changes in the residual vertical velocity following Pinatubo differ among reanalysis products and for varying methodologies used for computing the residual vertical velocity (Abalos et al., 2015), and hence the actual response of the BDC to Pinatubo is poorly constrained by existing reanalysis data.Specifically, Diallo et al. (2012) use diabatic heating rates in ERA-Interim to define w * , and tropical diabatic heating rates show cooling after Pinatubo in the reanalyses (as the reanalyses do not assimilate aerosol burden; www.atmos-chem-phys.net/17/1313/2017/cf.Fig. 1 of Abalos et al., 2015) but warm in GEOSCCM due to increased shortwave heating from the aerosol plume (not shown).Future work is needed in order to better constrain the response of the BDC to volcanic eruptions using observations.It is worth noting that the response to the eruption of Mt.Pinatubo is stronger in the NH mid-stratosphere than in the SH (cf.Fig. 6).The likely cause of this is as follows: Mt.Pinatubo is located at 15 • N, and in our experiments the majority of the aerosols stay in the NH.Therefore, the increased shortwave heating is stronger in the NH.Future work is needed in order to explore sensitivity to the details of the prescribed volcanic forcing in the model.

Summary of the comparison of GEOSCCM to observations
In conclusion, the evolution of the BDC in GEOSCCM is qualitatively, and by most measures considered here quantitatively, consistent with observational constraints.There is a transition between declining mean age throughout the stratosphere before the late 1980s and regionally-specific changes in mean age afterwards (including the possibility of aging in the midlatitude mid-stratosphere in the NH).The statement that is often made that climate models simulate a decreasing age throughout the stratosphere only applies over long time periods and is not the case for the past 25 years, for which we have most tracer measurements.

Conclusions
The Brewer-Dobson Circulation and its changes have important implications for both stratospheric and tropospheric climate as well as stratospheric ozone chemistry (SPARC-CCMVal, 2010;World Meteorological Organization, 2011, 2014;Manzini et al., 2014).Hence, it is crucial to understand (1) the structure of historical changes in the BDC and (2) the factors that lead to these changes.It is also important, for predicting future changes, to know how well models can simulate historical changes of the BDC as given by available observational constraints.Analysis of a series of chemistry-climate model experiments of the period January 1960 through December 2014 yielded the following conclusions: 1.Over the full duration of the experiments (i.e., for a start date in 1960), we recover the result from previous modeling studies: anthropogenic climate change leads to acceleration of the BDC throughout the stratosphere.Ozone depletion, rising GHG concentrations, and warming SSTs all led to declining mean age in all regions of the stratosphere.
2. Since the late 1980s, structural changes have occurred in the BDC.The BDC accelerated in the lower stratosphere in the NH and tropics but not in the midstratosphere.Specifically, since 1992 in the ensemble mean, mean age has increased by 0.12±0.09year (95 % confidence intervals) in the mid-stratosphere of the midlatitude NH, and tropical mass upwelling has slowed down by 2 % (statistically significant in two ensemble members).The trend in one ensemble member is quantitatively similar to that in available observations.Hence, there is no inconsistency in trends between our model and available observations.
3. The source of this structural change is the time-varying evolution of the forcing factors.While warming SSTs and rising greenhouse gas concentrations both lead to acceleration of the BDC (consistent with previous work), their influence is stronger in the lower stratosphere.In contrast, volcanic eruptions and ODS concentrations generally impact the deep branch more strongly.
Declining ODS concentrations and the proximity of the start of declining ODS concentrations to the eruption of Pinatubo led to an aging trend beginning in the early 1990s in the midlatitude NH mid-stratosphere.If internal atmospheric variability is taken into consideration, then the start date of an aging trend in the midlatitude NH mid-stratosphere can be pushed back to 1988.
In light of these results, we wish to emphasize that, if one wishes to understand the causes of observed historical changes in the BDC, careful attention must be paid to the start and end dates used for trend calculation and the forcings included in a model simulation.In addition, it should not be expected that any single free-running model experiment should capture the precise observed trend or necessarily resemble a second integration using that same model, as the wave forcing of the BDC differs in any realization of the atmospheric state.
Many questions as to the historical changes in the BDC are left unanswered by this study.Diagnostic output necessary to compute the full age spectrum was not saved for these model experiments, and hence we are limited in our ability to quantify mixing changes, but it is conceivable that mixing changes contributed to recent observed mean age trends (Ray et al., 2014;Ploeger et al., 2015).Future work is needed in order to better constrain the response of the BDC to volcanic eruptions using observations and to explore sensitivity to the details of the prescribed volcanic forcing in the model.Finally, a more quantitative comparison of model mean age from a range of modeling centers to recent in situ observations is needed in order to more firmly diagnose areas of agreement and disagreement between models and observations.

Data availability
Data related to this article can be found in the Supplement.Additional data requests should be addressed to Chaim Garfinkel (chaim.garfinkel@mail.huji.ac.ikl).
The Supplement related to this article is available online at doi:10.5194/acp-17-1313-2017-supplement.

Figure 1 .
Figure 1.Forcing applied in the simulations.(a) 10 • S-10 • N average sea surface temperatures, with a smoothed version of the curve bolded; (b) atmospheric concentrations of CO 2 following RCP4.5;(c) equivalent effective stratospheric chlorine (EESC; Newman et al., 2007, Eq. 1 with a = 1 and α for Br y = 60, and using 3year mean age); (d) ensemble mean of the aerosol optical thickness from explosive volcanic eruptions, resulting from prescribed injections of volcanic SO 2 ; (e) total solar irradiance.

Figure 2 .
Figure 2. Trends in annual averaged BDC from 1960 to 2014 in the ensemble mean of the three all-forcing integrations.Mean age trends are indicated by contours with a contour interval of 6 days decade −1 , and the residual circulation trends are indicated with streamlines.The thickness of the streamline is proportional to the magnitude of the wind speed.

Figure 3 .
Figure 3. Trend through 2014 in annual averaged mean age in the NH (left), tropics (middle), and SH (right) for start dates starting in 1960 and ending in 2002.Each line is for one member of the all-forcing ensemble or for the all-forcing ensemble mean.Changes in the midstratosphere (deep branch) are shown in (a-c), and changes in the lower stratosphere (shallow branch) are shown in(d-f).95 % confidence intervals on the calculated trends are indicated; note that in the SH mid-stratosphere the correction ofSanter et al. (2008) for the number of degrees of freedom leads to very few degrees of freedom for the trend and therefore to high uncertainty.

Figure 4 .
Figure 4.As in Fig. 3 but for the upwelling mass flux in the tropics.

Figure 5 .
Figure 5. (left) As in Fig. 2 but for trends in mean age from 1988 to 2014 in the three ensemble members.Note that (a) corresponds to the ensemble member with decreasing mean age and (e) to the ensemble member with aging air in the NH midlatitude stratosphere.The right column shows the changes in total wave forcing of the BDC (EP flux divergence plus gravity wave drag).

Figure 6 .
Figure 6.Time series of annual averaged mean age in the NH (left), tropics (middle), and SH (right).Each line is for one ensemble mean.Blue is all forcing, red is SST only, and intermediate colors span the other three experiments performed.Changes in the mid-stratosphere (deep branch) are shown in (a-c), and changes in the lower stratosphere (shallow branch) are shown in (d-f).

Figure 7 .
Figure 7. Trends in midlatitude NH mean age in the ensemble mean of each of the five experiments (a) from 1960 to 1992 and (b) from 1992 to 2014.The trend over the full duration of the experiments is shown in panel (c).The horizontal line indicates the ensemble mean trend, and the vertical line indicates the 95 % uncertainty bounds.

Figure 9 .
Figure9.Time series of annual averaged mass flux between the turnaround latitudes at 50 hPa (top) and 100 hPa (bottom).Each line is for one ensemble mean.Blue is all forcing, red is SST only, and green is for SSTGHGODS.For clarity, we suppress two of the intermediate experiments.

Figure 10 .
Figure10.As in Fig.7but for tropical mean upwelling at a variety of pressure levels.

Figure 11 .
Figure 11.Mean age estimates in the data from Fig. 7 of Ray et al. (2014) in the NH mid-stratosphere between 20 and 25 km and in the three all-forcing integrations for the same days.The trend from 1992 to 2012 is included.Note that in computing the trends we do not include any information as to the difference in uncertainty among the measurements for the Ray et al. (2014) data because all samples in GEOSCCM have the same uncertainty, and we prefer to use the same statistical test for both data sources for consistency.This leads to overly conservative estimates on the uncertainty for theRay et al. (2014) data.

Figure 12 .
Figure12.Mean age estimates in the data fromEngel et al. (2009) in the NH mid-stratosphere between 30 and 5 hPa and in the three all-forcing integrations for the same locations and months (results are similar if we use 30 to 10 hPa).The uncertainty for the observational estimates is taken fromEngel et al. (2009), and the uncertainty of the model simulated mean age can be deduced from the intra-ensemble spread.The GEOSCCM mean age is offset by 3 months, i.e., the bias in the mean age (which is less than the 6month uncertainty in the observed mean age as evaluated byEngel et al., 2009).

Figure 13 .
Figure 13.As in Fig. 2 but for changes in the BDC from January 2002 to December 2011 in the annual mean of the three allforcing ensembles.Note that the contour interval differs from Fig. 2.