The global diabatic circulation of the stratosphere as a metric for the Brewer–Dobson circulation

The circulation of the stratosphere, also known as the Brewer–Dobson circulation, transports water vapor and ozone, with implications for radiative forcing and climate. This circulation is typically quantified from model output by calculating the tropical upwelling vertical velocity in the residual circulation framework, and it is estimated from observations by using time series of tropical water vapor to infer a vertical velocity. Recent theory has introduced a method to calculate the strength of the global mean diabatic circulation through isentropes from satellite measurements of long-lived tracers. In this paper, we explore this global diabatic circulation as it relates to the residual circulation vertical velocity, stratospheric water vapor, and ozone at interannual timescales. We use a comprehensive climate model, three reanalysis data products, and satellite ozone data. The different metrics for the circulation have different properties, especially with regards to the vertical autocorrelation. In the model, the different residual circulation metrics agree closely and are well correlated with the global diabatic circulation, except in the lowermost stratosphere. In the reanalysis products, however, there are more differences throughout, indicating the dynamical inconsistencies of these products. The vertical velocity derived from the time series of water vapor in the tropics is significantly correlated with the global diabatic circulation, but this relationship is not as strong as that between the global diabatic circulation and the residual circulation vertical velocity. We find that the global diabatic circulation in the lower to middle stratosphere (up to 500 K) is correlated with the total column ozone in the high latitudes and in the tropics. The upper-level circulation is also correlated with the total column ozone, primarily in the subtropics, and we show that this is due to the correlation of both the circulation and the ozone with upper-level temperatures.

1 Introduction

The Brewer–Dobson circulation (BDC) is important for the distribution of trace gases in the stratosphere (Butchart, 2014) including water vapor, the radiative effects of which have been shown to impact surface climate (Dessler et al., 2013), and ozone, which impacts tropospheric circulation (e.g., Polvani et al., 2011) and human health (e.g., Abarca and Casiccia, 2002). In models and reanalysis products, the BDC is frequently quantified by the vertical velocity in the transformed Eulerian mean (TEM) framework (Andrews et al., 1987) averaged over the well-mixed tropics (e.g., Butchart et al., 2006; Li et al., 2008; Seviour et al., 2012; Hardiman et al., 2017). In steady state, the total upwelling and downwelling mass fluxes must be equal, so characterizing the tropics alone is considered sufficient. The TEM framework provides formalism that approximates the Lagrangian-mean mass transport, and in the limit of adiabatic, small-amplitude eddies, the TEM residual mean circulation is equivalent to the density-weighted isentropic mean circulation. Multimodel comparisons (Butchart et al., 2010) and inter-reanalysis comparisons (Abalos et al., 2015; Kobayashi and Iwasaki, 2016) have used the residual mean circulation at 70 hPa, averaged in the tropics, as a metric to evaluate the mean and trends of the BDC.

The 70 hPa level is consistently within the stratosphere even in climate models that do not accurately simulate tropopause height. As it is in the lower stratosphere, it approximates the mass flux between the troposphere and stratosphere and, as such, is related to water vapor flux and ozone transport.

Models predict that the residual mean circulation through a given pressure surface will increase in the future by about 2 % decade⁻¹ in the lower stratosphere and about 1 % decade⁻¹ in the middle and upper stratosphere (Butchart et al., 2010). This is a natural consequence of the lifting of the atmospheric circulation (e.g., Singh and O'Gorman, 2012; Oberländer-Hayn et al., 2016), and there are also dynamical reasons why one might expect a true acceleration of the BDC (e.g., McLandress and Shepherd, 2009; Garny et al., 2011; Shepherd and McLandress, 2011). However, observations have not shown such a robust trend (e.g., Stiller et al., 2012; Haenel et al., 2015; Engel et al., 2017). This disagreement can be attributed partially to the large internal variability in the system that prevents a 2 % decade⁻¹ trend from being detected without almost 30 years of data (Hardiman et al., 2017) and partially to the fact that there is no truly “like-to-like” comparison; a modeled tracer that is sampled like the observations can also fail to show a trend even when such a trend exists in the model (Garcia et al., 2011). Models also show that polar ozone loss has dampened the acceleration of the circulation, with an asymmetric effect on the different hemispheres (Polvani et al., 2018).

The TEM vertical velocity, which shows a robust trend in models, is a useful metric for understanding stratospheric dynamics. However, apart from its theoretical relationship with the Lagrangian-mean mass transport, it is not straightforward to relate the TEM vertical velocity to the tracer transport that is so important to climate due to the presence of other transport processes such as mixing (Ray et al., 2010, 2016; Dietmüller et al., 2017, 2018). In contrast, the global average diabatic overturning circulation through isentropes can be theoretically related to observed tracer distributions through the idealized tracer “age of air” (Neu and Plumb, 1999; Linz et al., 2016). This global diabatic circulation has been calculated from two different satellite data products (Linz et al., 2017), thus motivating the use of the global diabatic circulation as a metric for the BDC strength in addition to the TEM vertical velocity. In this paper, we explore differences between the global diabatic circulation and other calculations for the strength of the circulation in order to understand the relationship of this new constraint to more common metrics.

The calculation of the global diabatic circulation in Linz et al. (2017) is the first of its kind, but not the first observational estimate of the stratospheric circulation strength. Water vapor is transported into the stratosphere through the cold tropical tropopause, which has a strong seasonal cycle in temperature. The resulting time series of water vapor at the cold-point tropopause has a similarly strong seasonal cycle. By tracking the upward movement of the dry and wet phases over time, the water vapor signal, which is nearly conserved above the cold-point tropopause, can be used to calculate an effective velocity (w_TR). “Effective” refers to the aggregated transport, which includes the effects of advection and mixing. As a result, this “water vapor tape recorder” (Mote et al., 1996) method must be used with caution when studying the tropical tropopause layer (Podglajen et al., 2017) and with even more caution when comparing models (Dietmüller et al., 2018). This study minimizes such issues by focusing on the region above the tropical tropopause layer and by using a zonal mean between 10^∘ S and 10^∘ N to reduce the influence of horizontal mixing between the subtropics and the midlatitudes. We will explore the relationship between the vertical velocity derived from water vapor in the deep tropics and the global diabatic circulation.

One of the primary motivations for studying the BDC and its variability is its influence on stratospheric ozone. The circulation is known to transport ozone – this is why Dobson proposed it in the first place (Dobson et al., 1929), even if he concluded that this circulation was far-fetched. While the qualitative description of the influence of the stratospheric circulation on ozone variability is well established – variations in the transport of ozone from its maximum production location in the middle stratosphere in the tropics to the midlatitudes and poles – quantifying this effect is not simple. Furthermore, the interplay between the temperature, ozone, and circulation can lead to complex feedbacks. We know from observational studies that changes to dynamical quantities impact polar ozone (Hassler et al., 2011) and that the ozone hole recovery is currently being modulated by the dynamics (Solomon et al., 2016). In turn, variability and trends in the ozone affect the circulation (e.g., Polvani et al., 2011; Bandoro et al., 2014). In the Northern Hemisphere, the variability in hemispherically averaged upward Eliassen–Palm (EP) flux at 100 hPa from the early NCEP reanalysis data product has been shown to explain about 50 % of the interannual variability of total column ozone in wintertime (Fusco and Salby, 1999), with the influence of the wave driving dependent on the latitude (Reinsel et al., 2005). These strong relationships are a motivating factor in using the TEM residual mean vertical velocity, which is directly related to the EP flux, as a metric for the BDC strength. The global diabatic overturning circulation on isentropes has been demonstrated to be related to tracer transport more directly, but its relationship with ozone is unknown.

This paper serves to explore the global diabatic circulation as a metric for the stratospheric circulation strength. Section 2 describes the model runs, reanalysis products, satellite data, and regression methods. In Sect. 3, we provide an explanation of the steps for calculating the global diabatic circulation on isentropes, the necessary model output, and its advantages and disadvantages. In Sect. 4, we examine three different calculations for the TEM vertical velocity, including the underlying assumptions, with different tropical averaging. Then we compare the global diabatic circulation to the more traditionally used TEM vertical velocity calculated in these three different ways (Abalos et al., 2015) for three different reanalysis products and for the Whole Atmosphere Community Climate Model (WACCM). Thus, we determine how the information provided by this new metric compares to the information more typically used. We find close agreement between the global diabatic circulation strength and one of the three calculation methods for the TEM vertical velocity for the reanalyses (regardless of averaging choice); we also find close agreement between the global diabatic circulation strength and all three calculations for the TEM vertical velocity in the model (though only with fixed-latitude tropics). In Sect. 5, we compare the tropical vertical velocity calculated from the water vapor tape recorder (Niwano et al., 2003) from the WACCM model to the total overturning circulation. Similar to the good agreement found for the modeled residual circulation and global diabatic circulation, the global diabatic circulation strength and the water vapor tape recorder are closely linked in the model, although the correlation is weaker. In Sect. 6, we examine the relationship between the diabatic overturning circulation and stratospheric ozone using data, reanalyses, and WACCM. We find that the lower branch of the circulation is important for polar ozone, while the upper branch is the most important for subtropical ozone. The latter relationship is driven by the temperature dependence of the photochemistry and covariance of temperature and the global diabatic circulation. The ozone results are consistent with known relationships between TEM vertical velocity and ozone, demonstrating that the global diabatic circulation is just as good a metric for ozone variability. Section 7 summarizes the results and discusses implications and future work.

2 Model, reanalysis products, satellite data, and methods

A summary of the products, their resolutions, and associated references is given in Table 1.

For the model in this study, we use the Whole Atmosphere Community Climate Model (WACCM), a state-of-the-art chemistry–climate model. This model uses the physical parameterizations and finite-volume dynamical core (Lin, 2004) from the Community Atmosphere Model version 4 (Neale et al., 2013), with a vertical extent from the surface to the lower thermosphere and 31 pressure levels from 193 to 0.3 hPa. The WACCM simulation is the first member of an ensemble run based on the Chemistry Climate Model Initiative REF-C1 scenario (Morgenstern et al., 2017) and is forced with prescribed observed sea surface temperatures. This model simulation was shown to have a global diabatic circulation strength that closely agrees with the satellite tracer observations at 460 K (Linz et al., 2017). This study covers the time period 1980–2014.

Three different reanalysis data products are used in this study, building upon the work by Abalos et al. (2015) and Linz et al. (2017). These are the ECMWF Reanalysis Interim (ERA-Interim; Dee et al., 2011), the Modern Era Retrospective analysis for Research and Applications (MERRA; Rienecker et al., 2011), and the Japanese 55-year Reanalysis (JRA55; Kobayashi et al., 2015). Reanalyses are used for the same time period as WACCM for consistency in the comparisons. These reanalyses are based on the assimilation of different sets of data into three different models and using different assimilation schemes, leading to some significant differences in their output, especially above 10 hPa. Beneath 10 hPa, Abalos et al. (2015) showed that more uncertainty arose from the choice of calculation method for the vertical velocity than from the choice of reanalysis, concluding that differences between reanalyses were relatively small (except for trends). We build upon that result here and suggest that because of uncertainties in radiative heating rates, the reanalyses are not as robust in certain contexts.

Finally, we consider the total column ozone measurements from the Solar Backscatter Ultraviolet Instrument (SBUV) for 1980–2013 from the version 8.6 SBUV ozone data record (McPeters et al., 2013). This data product is based on nine recalibrated SBUV instruments with total column ozone measurements that are consistent with ground-based observations of total column ozone to within 1 %. We use the total column ozone as it has the least uncertainty for use in long-term correlation calculations.

Table 1 Model output, reanalysis products, and ozone data used in this study.

As the primary motivation of this paper is to evaluate relationships between dynamical and tracer quantities, it utilizes correlations extensively. The standard Pearson correlation coefficient is reported for each relationship. Only results significant at the 95 % level or greater are reported. Time series are deseasonalized by removing the climatology of each variable. Cross-correlations are used to examine the differences in the vertical structures of the different quantities. When these cross-correlations are between transport metrics that have different vertical coordinates, a climatological relationship between tropical (20^∘ S–20^∘ N) potential temperature and pressure is shown.

Figure 1 The seasonal cycle (a) and interannual variability (b) for the 450 K global diabatic circulation and the 750 K global diabatic circulation in the WACCM model from 1980 to 2014. Note that in (b) the sign of the anomalies has been reversed for the 450 K level in order to see the agreement. Panel (c) shows the coherence between these two time series, demonstrating that the visual correlation evident in (b) is related both to the quasi-biennial oscillation and to higher frequencies.

In the stratosphere, correlations of circulation metrics might be expected to be driven by the quasi-biennial oscillation (QBO) in addition to the seasonal cycle. Rather than explicitly removing this, we account for it by examining filtered time series and coherence (e.g., Fig. 1) and highlight the cases in which this is important. Many of the relationships examined are coherent at timescales shorter than the 2–3-year QBO period, though coherence is particularly high at that period. The relationships of dynamical variables with trace gases have less high-frequency variability and therefore tend to be dominated by the QBO.

3 Calculating the global diabatic circulation on isentropes

Why would we need a different metric for the BDC? The residual mean tropical upwelling at 70 hPa has been used for at least a decade (Butchart et al., 2006). However, it is neither directly observable nor easily relatable to observations. A metric for models and reanalyses should ideally be able to be constrained by data. The tropical leaky pipe model (Neu and Plumb, 1999), a three-box model of the stratospheric circulation that treats the tropics as largely isolated from the extratropics, results in the conclusion that the difference between midlatitude age and tropical age is related to the circulation. Linz et al. (2016) translated this model into isentropic coordinates to show a direct relationship between the idealized age of air (Hall and Plumb, 1994) and the diabatic circulation through an isentropic surface, demonstrating that the difference between the age of air that is downwelling and the age of air that is upwelling through each isentrope is inversely proportional to the diabatic mass flux through that surface, in statistically steady state and neglecting diabatic diffusion. Thus, the global diabatic circulation through an isentrope reflects the total tracer flux and should be considered an alternative, or at least additional, metric. This global diabatic circulation can also be calculated from satellite data.

3.1 Definition of the global diabatic circulation

The time-dependent, global, diabatic overturning mass flux through an isentrope is defined to be the average of the upwelling and downwelling mass fluxes, as follows.

As in Linz et al. (2016), we define the total upwelling mass flux, ℳ_u, and the total downwelling mass flux, ℳ_d, through an isentropic surface:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathcal{M}}_{\mathrm{u}}=\underset{\mathrm{up}}{\int }\mathit{\sigma }\dot{\mathit{\theta }}\mathrm{d}A,\text{and}\text{(2)}& {\displaystyle }{\mathcal{M}}_{\mathrm{d}}=-\underset{\mathrm{down}}{\int }\mathit{\sigma }\dot{\mathit{\theta }}\mathrm{d}A.\end{array}$$

$\dot{\mathit{\theta }}$ is the total diabatic heating rate, and $\mathit{\sigma }=-g^{-\mathrm{1}}\partial p/\partial \mathit{\theta }$ is the isentropic density. The limits of integration are the regions of the isentrope through which air is upwelling ($\dot{\mathit{\theta }}\mathit{>=}\mathrm{0}$) and downwelling ($\dot{\mathit{\theta }}<\mathrm{0}$) instantaneously. Since the monthly mean is not in equilibrium, some amount of storage may take place, and these two will not necessarily be equal. We therefore define the total overturning circulation as the average:

$$\begin{array}{}\text{(3)}& \mathcal{M}\left(\mathit{\theta }\right)=\left({\mathcal{M}}_{\mathrm{u}}-{\mathcal{M}}_{\mathrm{d}}\right)/\mathrm{2}.\end{array}$$

This is an arbitrary but sensible definition. Although one could certainly consider the extratropics or tropics alone, the treatment in Eq. (3) accounts for simultaneous variability in the extratropics and in the tropics, thus providing an instantaneous global average overturning circulation strength. Furthermore, it is this quantity that is directly related to the age of air distribution (Linz et al., 2016).

A note about the use of isentropic coordinates: the isentropic framework makes the separation of the diabatic and adiabatic components completely natural – they are simply the vertical and horizontal motions, respectively. In the annual mean and in steady state, the circulation on isentropes is much the same as the circulation on pressure surfaces. The seasonal variability of circulation on pressure surfaces and on isentropes differs, however. For example, the isentropes descend at the poles during springtime, which leads to upward motion of the air relative to the isentropes. Such springtime polar upwelling is not visible if pressure surfaces are used instead. Seasonal variability is removed from all time series in this study, and thus we attempt to minimize the effect of this difference on the comparisons. For trends, however, the longer-term motion of the isentropes may well be different from the motion of the pressure surfaces, which will be moving up as the tropopause lifts (e.g., Singh and O'Gorman, 2012). Thus, we might expect trends to have significantly different results depending on the choice of coordinate system, perhaps as different from trends through pressure surfaces as those calculated relative to the tropopause height (Oberländer-Hayn et al., 2016).

To calculate the global diabatic circulation from model output or reanalysis, one needs the total diabatic heating rate, the temperature, and the pressure. The diabatic heating rate is output differently in different models, but it is straightforward. The diabatic heating rate consists predominantly of two terms, the latent heat flux from phase changes of water and the radiative heating and cooling (Fueglistaler et al., 2009; Wright and Fueglistaler, 2013). For levels wholly within the stratosphere, water vapor concentrations are so low that the former is negligible. Models may output other diabatic terms, such as mixing from parameterized gravity waves; alternatively, they may output a total temperature tendency, which contains all of the necessary information in just one term. Almost all models will output these terms on either pressure or model levels. The diabatic heating rate on those levels must then be interpolated to isentropic levels, for which the temperature and pressure fields are necessary. The isentropic density can be calculated by finite difference in pressure and then interpolated to isentropes as well.

Since eddies serve to predominantly mix adiabatically, they are (naturally) less important for the global diabatic circulation than for the residual circulation. In the conversion from the diabatic vertical velocity on pressure surfaces to the diabatic vertical velocity on isentropes, the covariance of the diabatic vertical velocity and the isentropic levels could nevertheless be important. However, this covariance is small enough that monthly mean temporal resolution is sufficient to accurately calculate the circulation; specifically, in ERA-Interim using monthly means instead of 6-hourly means results in no bias throughout most of the stratosphere and up to a 10 % bias at the poles in wintertime, which, as the pole is a small area of the globe, is a much smaller bias on the total overturning mass flux. While many models do output monthly mean eddy fluxes to calculate the residual circulation, others, especially older model runs, do not. Almost all models output shortwave and longwave radiation, and as these are by far the dominant terms in the total diabatic heating rate, this metric can be calculated using models that did not report the necessary terms or have the necessary temporal resolution for the residual circulation vertical velocity calculation. The comparatively minimal data requirements for this metric recommend it for intermodel comparisons.

Although the global diabatic circulation strength is a good indication of the integrated eddy forcing on the circulation, it does not diagnose which eddies are responsible. Thus, models could get the right circulation from the wrong waves, and, indeed, there is reason to expect compensation between resolved and parameterized wave driving (Cohen et al., 2013). Because of this compensation, the analysis of different wave forcing contributions to the BDC in the residual mean framework is also potentially problematic. Finding an appropriate way to relate any BDC metric directly to tropospheric forcing in a way that could inform model development and tuning is an interesting area of research.

3.2 Global diabatic circulation characteristics

The mean value of the global diabatic circulation at 460 K for WACCM is 7.1×10⁹ kg s⁻¹, decreasing to 1.8×10⁹ kg s⁻¹ at 1000 K (Linz et al., 2017). The seasonal cycle, which is subsequently removed, is shown in Fig. 1a for two different levels for the global diabatic circulation from WACCM. The lower stratosphere has a single peak, while the upper stratosphere has a semi-annual cycle as well. This climatology is subtracted to obtain the time series shown in Fig. 1b. Note that the negative of the anomaly is plotted for the lower level to enable a clear comparison of these two anticorrelated time series. The different timescales of variability are visible, with an obvious QBO signal and shorter timescale variability. Although the correlation between the upper and lower levels is clear and in phase at QBO timescales, the higher-frequency variation is also correlated, but with a 20–90^∘ phase lag (not shown). The coherence between these two time series is shown in Fig. 1c. There is high coherence at periods of 2–3 years, as expected with the QBO. There is also coherence for periods of shorter than about 9 months, which is unexplained by the QBO.

Figure 2 Correlation coefficient (r) for the autocorrelation of the deseasonalized time series of (a) the global diabatic circulation and of the three different TEM vertical velocities calculated from WACCM with 30^∘ tropics: (b) ${\overline{w}}^{*}$, (c) ${\overline{w}}_{\mathrm{M}}^{*}$, and (d) ${\overline{w}}_{\mathrm{Q}}^{*}$, and with the true turnaround latitudes (e) ${\overline{w}}^{*}$, (f) ${\overline{w}}_{\mathrm{M}}^{*}$, and (g) ${\overline{w}}_{\mathrm{Q}}^{*}$. As the diagonal reflection is redundant, it is not shown. Contours are spaced every 0.1.

The vertical autocorrelation coefficient (r) of the global diabatic circulation (with the seasonal variability removed) is shown for WACCM in Fig. 2a. The autocorrelation is relatively narrow in width, showing that the variability of the lower-level circulation is relatively uncorrelated with that of the upper-level circulation. An interesting feature is the weak anticorrelation between lower and upper levels, which can also be seen in the vertical autocorrelation function of the heating rates themselves (in either pressure or isentropic coordinates). Some of this anticorrelation is due to the anticorrelation at QBO timescales, but the higher-frequency variability is also anticorrelated, as can be seen from Fig. 1, and the dynamical reasons for this are the subject of ongoing investigation.

3.3 Global diabatic circulation trends

The trends in the global diabatic circulation have not previously been examined. We calculate the trends (1980–2014) in the global diabatic circulation from the three different reanalyses and also from the WACCM model run over the same time period, and the results of this are shown in Fig. 3. These results are similar to those found by Abalos et al. (2015) for the TEM vertical velocity calculated using the thermodynamic equation, ${\overline{w}}_{\mathrm{Q}}^{*}$. Because of the different coordinate system, however, some differences exist. (Note that Abalos et al., 2015, found that the removal of interannual variability does not change the long-term trends significantly.) ERA-Interim shows an acceleration of the lower branch of the circulation and a deceleration of the upper branch. MERRA shows an acceleration around the mid-stratosphere, where the upper branch begins, and in the uppermost stratosphere. JRA55, meanwhile, only has a small region in the mid-stratosphere where it shows a statistically significant trend. This is also an acceleration. The WACCM simulation for this time period has no statistically significant trend in the global diabatic circulation at any level, despite significant trends in the thermal structure. This result of no trend in the WACCM overturning is intriguing; although the isentropic levels are changing location over these decades, the total overturning through each isentrope is not significantly changing. This is perhaps related to the lifting of the circulation described by Oberländer-Hayn et al. (2016) so that the circulation strength stays the same through isentropes but moves upwards in pressure. This is consistent with the results of Abalos et al. (2017), who found that trends in the residual stream function for a WACCM model run from 1955 to 2099 were far weaker when calculated with respect to the tropopause (though the trends were still significantly positive over that long period). The differences in trends in this metric compared with the more standard TEM vertical velocity calculation (Andrews et al., 1987), which show significant positive trends at most levels for MERRA and JRA55 regardless of the definition of the tropics (Abalos et al., 2015), demonstrate that although the global diabatic circulation varies closely with the other metrics, trends will appear different because changes to the thermal structure and the circulation play a role. Note that since the time series of heating rates in reanalyses are somewhat questionable above 800 K, where they are influenced by changes in the observing systems (Simmons et al., 2014), the trends there are to be treated with caution.

Figure 3 Trends in the global diabatic circulation at each level calculated from the three reanalysis data products (JRA55: blue; ERA-Interim: yellow; MERRA: green) and from the WACCM model (black). Dashed lines show the calculated trends that are not significant at the 95 % confidence level, while bold lines are significant. There are no significant trends in the WACCM model run.

4 The global diabatic circulation and TEM vertical mass flux in three reanalyses and a model

The BDC was originally hypothesized to explain observed tracer distributions (Dobson et al., 1929; Brewer, 1949), and therefore the Lagrangian mean transport is, in some sense, the appropriate formalism to study. The TEM residual circulation is not the same as the Lagrangian mean mass transport, as noted explicitly in Andrews and McIntyre (1976). However, under certain conditions (small amplitude, adiabatic eddies), the Lagrangian mean circulation and the TEM residual circulation are identical. The TEM equations also provide unique insight into the forcing of the mean flow by eddies; when the quasi-geostrophic approximation holds, the internal forcing of the mean state by the eddies is encompassed by the divergence of the Eliassen–Palm flux (Edmon et al., 1980). Because of the ready interpretation of the wave-mean flow interactions, the TEM residual mean circulation has been the primary diagnostic of the stratosphere for models. It cannot, however, be derived from observations. Here, we try to understand differences in the common methods for calculating the TEM residual circulation vertical velocity and the relationship between it and the global diabatic circulation.

A note about the QBO: although the QBO influences both the residual circulation and the global diabatic circulation, the relationships between the metrics in this section are significantly coherent at all frequencies (see Fig. 7 for a comparison of time series of ${\overline{w}}^{*}$ and ℳ.)

4.1 Comparison of TEM vertical velocity calculation methods

Abalos et al. (2015) performed an extensive reanalysis intercomparison of the trends in the TEM vertical mass flux calculated in multiple ways from ERA-Interim, MERRA, and JRA-55. For this work, the calculations were repeated for the WACCM model output. The three different methods for calculating the mass flux are summarized below, and for more details see the original paper. The first method is the residual circulation (Andrews et al., 1987), ${\overline{w}}^{*}$, in which the residual vertical velocity is calculated based on the Eulerian mean vertical velocity and the meridional eddy heat flux. This method, which we will refer to as the “direct” method relies on performing vertical integrals of the velocity fields from reanalyses. The second calculation of the BDC, ${\overline{w}}_{\mathrm{M}}^{*}$, is based on the “downward-control” principle (Haynes et al., 1991) and is calculated using the momentum balance equation by integrating the difference of the divergence of the Eliassen–Palm flux and the zonal mean zonal wind tendencies on surfaces of constant “angular” momentum (in this case, constant latitude) to derive a stream function (Randel et al., 2002). The assumption of isolines of angular momentum being equivalent to latitude lines could lead to errors in this estimate. Both of these methods also rely on the applicability of the quasi-geostrophic approximation to interpret their results. The final estimate, ${\overline{w}}_{\mathrm{Q}}^{*}$, is calculated by iterating the thermodynamic balance and the continuity equation with no net mass flux across a pressure surface (Murgatroyd and Singleton, 1961; Rosenlof, 1995). Any errors in heating rates will be reflected in this calculation. Because this estimate is also derived from the heating rates, this should be the most closely related to the global diabatic circulation. For this study, we use the deseasonalized time series of these estimates of the BDC strength integrated over 30^∘ S–30^∘ N and integrated between the turnaround latitudes (Abalos et al., 2015, Fig. 8) at all levels throughout the depth of the stratosphere.

The first two of these methods both require at least 6-hourly data, while the thermodynamic ${\overline{w}}_{\mathrm{Q}}^{*}$ needs only monthly mean data (Lin et al., 2015). For the purposes of this study, the same frequency of data (6-hourly instantaneous values) was used for all three methods and then monthly averages were taken. The interpretations of the results in terms of eddy forcing for both the direct method and the downward-control method rely upon quasi-geostrophic balance, whereas the thermodynamic method does not. Thus, we might expect that the two quasi-geostrophic estimates derived from high-frequency data would be very similar. Butchart et al. (2006) calculated the mean and the trend in the residual vertical velocity using both methods in a variety of models and found that they were generally similar in magnitude and structure, though differences between the two calculations varied more than differences in the interannual variability of each individual calculation. Rosenlof (1995) compared the thermodynamically calculated stream function to the downward-control stream function and found them to be similar, but with the strongest differences in the lower stratosphere. Abalos et al. (2015) also performed a qualitative comparison between the mean stream function for these three estimates, noting that the thermodynamic calculation is larger in the lower stratosphere and with more differences between the downward-control calculation and the other two estimates higher up in the stratosphere at the poles.

In order to understand some of the properties of the different vertical velocity calculations, we examine the vertical autocorrelation of the tropical upwelling velocity for the WACCM model. Note that the downward-control calculation of ${\overline{w}}_{\mathrm{M}}^{*}$ includes the gravity wave drag, since this made an important contribution in the model (but not in the reanalyses, in which the gravity wave drag is much smaller).

Figure 2b–d show these autocorrelation coefficients for the three methods with the tropical averaging latitudes fixed between 30^∘ S and 30^∘ N. The direct method is shown in panel (b), and the correlation is broader than the equivalent autocorrelation for the global diabatic circulation. The downward-control method is shown in panel (c), and the vertical autocorrelation is even greater. The downward-control method means that upper levels directly impact lower levels (through integration), so it is consistent that the vertical autocorrelation of ${\overline{w}}_{\mathrm{M}}^{*}$ is the broadest of all metrics. Note then that the vertical velocity calculated in this way is essentially a single piece of information for the extent of the stratosphere. Differing variability in the upper and lower branches will be comparatively indistinguishable using such a calculation. Previous results suggest that the upper and lower branches may be distinguished by this metric for subseasonal variability in winter, however (Abalos et al., 2014). The thermodynamic vertical velocity in panel (d) demonstrates that the vertical covariance is not necessarily a result of the flow itself, since vertical correlation is much narrower in this case. Unlike the global diabatic circulation, however, there is little apparent anticorrelation between the upper and lower branches of the circulation. There is an interesting feature in the lower stratosphere for this radiatively determined vertical velocity; beneath 70 hPa, the behavior is much more weakly correlated with upper levels than for the other calculations of vertical velocity. This is consistent with the results of Rosenlof (1995), who speculated that the reason for this low-level discrepancy was the relatively simple way the radiative heating was calculated by using the radiative transfer code developed for two-dimensional models by Yang et al. (1991) and Olaguer et al. (1992). However, this different behavior in the lowermost stratosphere was also found by Abalos et al. (2015) with the same three complex reanalyses used here, and the result holds with the WACCM model here. These three calculations, often treated as the same, are actually somewhat different, especially with respect to the vertical structure of their variability.

Figure 2e–g show the autocorrelation coefficients for the three methods of calculating the vertical velocity with the tropical averaging latitudes set by the location of the turnaround latitudes determined each month from the location where ${\overline{w}}^{*}$ switches from upwelling to downwelling. These turnaround latitudes vary from narrower than 30^∘ at the lowermost levels to closer to 40^∘ at the upper levels (see Abalos et al., 2015, Fig. 5 for the mean and climatology of these in the reanalyses). Using the true turnaround latitudes instead of set latitudes for the tropics makes the vertical velocity calculated using all three methods have a narrower extent of the vertical autocorrelation. The implication of this is that a good deal of the difference in variability between levels occurs at the edges of the “pipe”, where mixing plays a role. The difference between the two different edge treatments is greatest in the direct calculation of the vertical velocity (panel e). This shows that the variability of the vertical velocity in the lower stratosphere is no longer positively correlated with the variability of the vertical velocity in the upper stratosphere. The lower stratospheric structure now resembles that of the fixed-latitude thermodynamic calculation for all three calculations, with very little relation between the variability beneath 70 hPa and above that level. This suggests that the difference between the thermodynamic calculation and the others is unrelated to the treatment of radiation. As above, however, we can conclude that the three different methods of calculation provide different vertical information.

To examine the relationship between the fixed-latitude and turnaround-latitude calculations, we show the cross-correlation between the two for each calculation method in Fig. 4. The y axis is the turnaround latitude and the x axis is the fixed-latitude calculation. It is evident that the different tropical boundaries matter most for the direct calculation method. The high degree of symmetry in panels (b) and (c) implies that, although small differences are visible in the autocorrelations in Fig. 2, the choice of boundary is far less important for the momentum and thermodynamic methods.

Figure 4 Cross-correlation of the variability in the calculations of the residual circulation with the turnaround latitudes (y axis) and the fixed tropics, 30^∘ S and 30^∘ N, (x axis). (a) Direct calculation, (b) momentum calculation, and (c) thermodynamic calculation. Contour lines are every 0.1.

In Fig. 5, we show the matrix of correlation coefficients (r) for each version of the residual circulation vertical velocity with each other version. Panels (a–c) show the WACCM turnaround-latitude calculation, panels (d–f) show the WACCM fixed-latitude calculation, and panels (g–i) show the fixed-latitude calculation from one of the reanalyses, ERA-Interim. (Behavior is similar for the other two reanalyses.) Panels (a), (d), and (g) show the correlation between the direct calculation on the y axis with the downward-control calculation on the x axis. Panels (b), (e), and (h) show the correlation of the thermodynamic TEM vertical velocity ${\overline{w}}_{\mathrm{Q}}^{*}$ with the downward-control calculation ${\overline{w}}_{\mathrm{M}}^{*}$. Panels (c), (f), and (i) show ${\overline{w}}^{*}$ on the y axis and ${\overline{w}}_{\mathrm{Q}}^{*}$ on the x axis. These correlation coefficients demonstrate the interchangeability (or lack thereof) of these different calculations for the vertical velocity. Examining the turnaround-latitude calculations (a–c), one notes that the correlation of the downward-control calculation with either of the other calculations is quite weak, never getting above r=0.9. We hypothesize that this is because the vertical integration, which smears out information in the vertical, makes the averaging using turnaround latitudes less clear, since the turnaround latitudes themselves vary with height. The comparison between the direct calculation and the thermodynamic calculation in Fig. 5c has much closer agreement than either comparison with the downward-control method. Correlations between the same vertical velocities in WACCM, but now with fixed averaging latitudes (30^∘ S–30^∘ N), are much higher.

Because the calculation for the vertical velocity averaged between turnaround latitudes is less well defined (sometimes there is more than one zero crossing, for example), and because the fixed-latitude calculation is simpler and therefore more common, we shall default to using the fixed-latitude calculation for the most of the rest of this study, though some comparisons with averaging between turnaround latitudes are included as well.

Now we focus on the lower six panels of Fig. 5 to see the differences between the methods with the fixed-latitude averages and the differences between the model and the reanalysis. In panel (d), the correlation of the two momentum-based calculations at the same level is very high, with the WACCM correlations appearing very similar to the autocorrelations in Fig. 2 and r>0.9 along the diagonal between 50 and 10 hPa for the reanalysis (g). We see evidence of the broad autocorrelation of ${\overline{w}}_{\mathrm{M}}^{*}$, as the correlations of the lower-level ${\overline{w}}^{*}$ with the upper levels of ${\overline{w}}_{\mathrm{M}}^{*}$ are much higher than the opposite. We note that when the full downward-control calculation – using contours of angular momentum rather than latitude lines – is applied to calculate the ${\overline{w}}_{\mathrm{M}}^{*}$ from ERA-Interim, the correlation with ${\overline{w}}^{*}$ is actually somewhat worse (r<0.7 along most of the diagonal; not shown) and even lower (r<0.3 along the diagonal) when the correlation is calculated with 6-hourly data rather than monthly (c.f. the impact of this calculation on the mean in Ming et al., 2016). We speculate that the worse agreement at higher frequencies is related to either small-scale torques that are not captured by the momentum budget at these high frequencies or due to the assumption of instantaneous net-zero flow through each pressure surface, which cannot account for short-term storage. In panels (e) and (h), the correlations with the downward-control calculation and the thermodynamic calculation again reach much deeper along one axis than the other, associated with the broad vertical autocorrelation of the downward-control calculation method. Interestingly, at upper levels in the model, this cross-correlation is strongest, while in the reanalysis, the upper levels are where the cross-correlation is weak. The weak correlation at upper levels in the reanalysis product could be a result of the discontinuities in the heating rates above 5 hPa noted by Abalos et al. (2015). The correlation beneath 70 hPa is weak in the model and is not significant in the reanalysis, again consistent with the substantial differences at low levels seen in the mean by both Rosenlof (1995) and Abalos et al. (2015). There are major discrepancies between the lower stratospheric heating rates in different reanalyses, which could explain this feature to some extent (Wright and Fueglistaler, 2013). Panels (f) and (i) show ${\overline{w}}^{*}$ on the y axis and ${\overline{w}}_{\mathrm{Q}}^{*}$ on the x axis. These compare more favorably than panels (b), (e), and (h), but it is important to note that even in the WACCM model with these fixed latitudes, these are not equivalent beneath 70 hPa. In the reanalysis, the correlation of these is a bit higher than for the comparison in panel (h), but again there is limited correlation in the upper stratosphere. Because of their different time evolution, it is not entirely surprising that the trends in the circulation calculated using these different methods disagree with each other for the reanalyses (Abalos et al., 2015).

Figure 5 Correlation coefficients (r) for the deseasonalized time series of the three different TEM vertical velocities calculated from WACCM (a–f) and for ERA-Interim (g–i): (a, d, g) ${\overline{w}}^{*}$ vs. ${\overline{w}}_{\mathrm{M}}^{*}$, (b, e, h) ${\overline{w}}_{\mathrm{Q}}^{*}$ vs. ${\overline{w}}_{\mathrm{M}}^{*}$ , and (c, f, i) ${\overline{w}}^{*}$ vs. ${\overline{w}}_{\mathrm{Q}}^{*}$. The quantity plotted on the y axis is listed first. (a–c) Averages between the turnaround latitudes and (d–f) averages between 30^∘ S and 30^∘ N. Contours are spaced every 0.1.

4.2 TEM vertical velocity compared to the global diabatic circulation

Next, we seek to answer the question of how the global diabatic circulation on isentropes relates to these metrics. We calculate the correlation of the three different calculations of the TEM vertical velocities averaged over 30^∘ S–30^∘ N with the deseasonalized global diabatic circulation on each isentrope (as defined above) for each of the three reanalysis data products and for the WACCM model. For WACCM and for JRA55, we also show these cross-correlations with the TEM vertical velocities averaged between the turnaround latitudes. (We show only JRA55 because its pattern is very similar to the other two reanalyses. ERA-Interim has slightly higher correlations throughout and MERRA has slightly lower correlations throughout.) These 18 correlation coefficient matrices are shown in Fig. 6.

Figure 6 The correlation of three different estimates of the TEM vertical velocity with the global diabatic circulation ℳ. (a–c) WACCM model with TEM vertical velocities averaged between the true turnaround latitudes; (d–f) WACCM model with fixed 30^∘ S–30^∘ N tropical averaging latitudes; (g–i) JRA-55 reanalysis averaged between turnaround latitudes; (j–l) JRA-55 averaged between fixed 30^∘ S–30^∘ N; (m–o) ERA-Interim averaged between fixed 30^∘ S–30^∘ N; (p–r) MERRA averaged between fixed 30^∘ S–30^∘ N. The first column is the correlation of the global diabatic circulation with ${\overline{w}}_{\mathrm{Q}}^{*}$; the second column is the correlation of the global diabatic circulation with ${\overline{w}}_{\mathrm{M}}^{*}$; the third column is the correlation of the global diabatic circulation with ${\overline{w}}^{*}$. The gray dashed line shows the climatological relationship between pressure and potential temperature averaged between 20^∘ S and 20^∘ N. Contours spaced every 0.1.

The highest correlation is found between the global diabatic circulation and ${\overline{w}}_{\mathrm{Q}}^{*}$ (as expected because these are both calculated from the heating rates) for all three reanalyses and the model. In addition, this comparison has the smallest vertical extent, consistent with the narrower extent of vertical autocorrelations seen in Fig. 2. The absolute highest correlations are between the global diabatic circulation and ${\overline{w}}_{\mathrm{Q}}^{*}$ averaged between turnaround latitudes in the WACCM model. Interestingly, when comparing the turnaround latitudes to 30^∘ S–30^∘ N for this cross-correlation in JRA55, the opposite result is seen than for the model. In JRA55 (and for the other two reanalyses; not shown), the correlation between the fixed latitudes is higher. This means that the turnaround-latitude averaging introduces more spurious variability in the reanalysis products, while in the model, using the true turnaround latitudes provides closer agreement with the global diabatic circulation. This seems only natural, since the global diabatic circulation is the average of the total mass flux through the surface instantaneously and therefore itself accounts for the motion of the turnaround latitudes.

For the correlation between the global diabatic circulation and ${\overline{w}}_{\mathrm{Q}}^{*}$ in the other three reanalysis calculations, the 50 hPa ${\overline{w}}_{\mathrm{Q}}^{*}$ variability is captured by the 450–500 K global diabatic circulation. The 10 hPa ${\overline{w}}_{\mathrm{Q}}^{*}$ variability is captured in all three reanalyses by the global diabatic circulation between 800 and 900 K. The climatology of the potential temperature–pressure relationship in the tropics (20^∘ N–20^∘ S) is shown by the dashed line. Note that because the diabatic circulation reflects the global circulation while vertical velocities are calculated only in the tropics, the highest correlations are not necessarily expected to be along this line, but it is a useful visual guide. In all three reanalyses and the model, there is some reflection of the anticorrelation of the upper and lower branches of the circulation that is seen in the global diabatic circulation on isentropic surfaces. The relationship with the other TEM vertical velocities is less clear in the reanalyses, though it is still quite strong in the WACCM model. In the reanalyses, ${\overline{w}}^{*}$ at 70 hPa is not strongly correlated with the global diabatic circulation at any level, with the correlation coefficient only reaching up to r=0.5 (at 550 K for both JRA55 and MERRA and between 550 and 650 K for ERA-Interim). The momentum-derived vertical velocity is the least well correlated, with the lower-level global diabatic circulation having almost no covariability with ${\overline{w}}_{\mathrm{M}}^{*}$ at any level except in WACCM. We conclude from this comparison that the global diabatic circulation is very closely related to the TEM vertical velocity calculated using heating rates with less covariation with ${\overline{w}}^{*}$ and even less with the momentum-derived vertical velocity, ${\overline{w}}_{\mathrm{M}}^{*}$. Similar to Abalos et al. (2015), we generally see as much difference amongst the different estimates of the vertical velocity as between the three reanalyses. The WACCM results demonstrate that the tropical upwelling averaged between 30^∘ S and 30^∘ N and the global diabatic circulation, while closely related, are not equivalent. Although the comparison for the thermodynamic vertical velocity with the global diabatic circulation is in places greater than 0.9, the comparison with the other TEM calculations reveals differences, especially lower in the stratosphere. When turnaround latitudes are used instead, the correlations become worse for the global diabatic circulation with both ${\overline{w}}_{\mathrm{M}}^{*}$ and ${\overline{w}}^{*}$. However, the correlation with ${\overline{w}}_{\mathrm{Q}}^{*}$ suggests that these two are very nearly identical, especially above the middle stratosphere. In a model, they could be used interchangeably, but in reanalysis, they are quite different.

5 The global diabatic circulation and the water vapor tape recorder

As discussed in the Introduction, the water vapor tape recorder can be used to calculate an effective velocity (w_TR) by tracking the seasonal cycle as it is moved along by the BDC and is another way to get at an “observed” circulation. We modify previous approaches (Niwano et al., 2003; Schoeberl et al., 2008) by using four levels (instead of two) for a phase-lagged correlation. This modification appears to better capture interannual variability (e.g., QBO), whereas the two-level method is better at capturing the seasonal cycle (Glanville and Birner, 2017).

We correlate monthly data between three lower levels (z to z+2) and three upper levels (z+1 to z+3) such that the two middle levels overlap. We then calculate the correlation coefficients, shifting the upper-level data from +1 to +9 months while holding the lower level still. The lag with the largest correlation coefficient represents the approximate time needed for the tape recorder signal to ascend from the lower levels to the upper levels. The tape recorder speed, assigned to the midpoints between the levels and the time steps, is simply the distance between the levels divided by the time lag. This modified method was tested on various scenarios within a 1-D model and was found to successfully capture variability but underestimate speeds by 5 %–10 %. This method of calculation improves the representation of interannual variability (like the QBO) compared with a simpler two-level method.

It should be noted that methane oxidation acts as a water vapor source, affecting mean values above 70 hPa (∼ 450 K) but with smaller impacts on the interannual variability up to about 10 hPa (Kawatani et al., 2014). Depending on the seasonal cycle of methane and the speed of the BDC, this can result in an apparent slowdown, speedup, or nothing at all. For example, if oxidation occurs before (after) the wet signal, the effective velocity will appear stronger (weaker). However, if oxidation is concurrent with the wet signal, the velocity calculation will not be affected.

Note that although reanalysis products do output water vapor, the inconsistencies of the water vapor tape recorder with the vertical velocities in reanalysis, likely due to enhanced dispersion from the assimilation process, lead us to omit their analysis (Glanville and Birner, 2017). The results of the water vapor tape recorder comparison to the global diabatic circulation are shown for WACCM at 500 K in the time series in Fig. 7. This figure shows the significant correlation between the water vapor tape recorder vertical velocity and the global diabatic circulation (r=0.57), and certain features stand out. The QBO appears to be related to a significant fraction of the covariation of these two time series, and when examined, the coherence drops off with periods shorter than the annual timescale. The water vapor tape recorder vertical velocities also appear to have greater decadal variability than the global diabatic circulation. The correlation of these improves upon filtering to remove the higher-frequency variability in the global diabatic circulation, which the water vapor tape recorder does not capture. The correlation between these two measures of the circulation is not strong enough for them to be considered equivalent, in the way that the WACCM results above suggest near equivalency between the global diabatic circulation and the tropical residual circulation vertical velocities for considering interannual variability. w_TR results from observations must be understood within this context.

Figure 7 Time series of calculated water vapor tape recorder (blue), the global diabatic circulation (black) on 500 K, and ${\overline{w}}^{*}$ averaged between 30^∘ S and 30^∘ N (red) from the WACCM model. All three time series have been deseasonalized and scaled by their standard deviations. Correlation coefficients between each pair of time series are reported at the bottom.

To examine the correlation more broadly, we show the cross-correlation between the water vapor tape recorder vertical velocity and the global diabatic circulation at every level in Fig. 8a. The correlation is around 0.5–0.6 along the diagonal, with an anticorrelation of up to 0.4–0.5 between the upper branch and lower branch (regardless of the metric). Interestingly, the correlation with the TEM tropical vertical velocity averaged between 30^∘ S and 30^∘ N is considerably weaker, as shown in panel (b). Note that when the correlation between the TEM tropical vertical velocity and w_TR is calculated in pressure coordinates, the magnitude of the correlation is the same as with the isentropic coordinates, except between 5 and 10 hPa where it is much weaker (not shown). When the TEM tropical vertical velocity is averaged between the true turnaround latitudes, however, the correlation becomes stronger than the correlation with the global diabatic circulation (panel c).

This combination of results – i.e., that the coherence drops off at periods less than a year and that the correlation of the water vapor tape recorder with the global diabatic circulation is stronger than with one type of averaging for ${\overline{w}}^{*}$ but weaker than with the other – suggests that the w_TR is mostly recording longer timescale variations, and its correlation with the other vertical velocity metrics is mostly to do with which ones respond the same way to the QBO. The anticorrelation seen in Fig. 8c is as strong as the correlation along the diagonal. Why the ${\overline{w}}^{*}$ averaged between the turnaround latitudes has a response to the QBO that is most similar to that of w_TR is unclear. The turnaround latitudes are the narrowest in latitude near the tropopause, where the w_TR signal is set, and perhaps this geometry matters. A takeaway from this is that, if one were to compare model results to water vapor observations, none of the dynamical vertical velocity metrics from the model would be appropriate comparisons. Instead, the model's water vapor tape recorder velocity would need to be used. This limits the usefulness of w_TR as an observable metric for evaluating reanalyses.

Figure 8 Correlation coefficient of the interannual variability of w_TR with (a) the global diabatic circulation and (b) the residual circulation vertical velocity, ${\overline{w}}^{*}$, at different levels in the stratosphere. The climatological relationship between pressure and potential temperature (averaged between 20^∘ S and 20^∘ N) is shown in (b) and (c) by the dashed gray line.

6 The global diabatic circulation's relationship with ozone

One motivation for studying the BDC is its influence on radiatively important trace gases, such as water vapor and ozone. Water vapor is a quasi-conserved tracer once it enters the stratosphere (in the absence of the aforementioned methane oxidation), so its behavior is comparatively straightforward. In contrast, ozone is both produced and destroyed in the stratosphere in chemical processes that are photochemically and temperature dependent. The ozone maximum is around 7 hPa or 800 K in the tropics (e.g., Paul et al., 1998), where photolysis by wavelengths less than 240 nm dissociates molecular oxygen (Chapman, 1930; Seinfeld and Pandis, 2006). As air moves from the tropics, it advects the ozone to middle and high latitudes. Stratospheric ozone absorbs ultraviolet radiation, creating heat, and thereby influences the thermal structure of the stratosphere (e.g., Andrews et al., 1987) and thus diabatic heating and transport. As the chemistry itself is temperature dependent, ozone, temperature, and circulation are closely connected.

With this interconnectivity in mind, we examine the total column ozone correlation at every latitude with the global overturning circulation at each level within the stratosphere. The correlation of the deseasonalized time series of the monthly mean total column ozone data from the SBUV from 1980 to 2013 and the global diabatic circulation from the three different reanalyses is shown in Fig. 9. Also shown is the correlation of the total column ozone and global diabatic circulation from the WACCM model. Generally, there is a consistent pattern across all three reanalyses and the model. This pattern is consistent with the ozone variability associated with the QBO: an out-of-phase relationship between the lower and upper stratosphere and an out-of-phase equatorial and subtropical pattern (e.g., Zawodny and McCormick, 1991). The ERA-Interim correlation with the SBUV data is much stronger than the correlations of the other two reanalyses with the SBUV data. Note that ERA-Interim assimilates the SBUV data, which MERRA and JRA55 do not, and this is a likely explanation for the higher correlation. Nevertheless, as the same spatial patterns are visible in the correlations with all three reanalyses and the model, we consider them to be robust and seek to understand them, i.e., whether they are due almost entirely to the QBO as with water vapor or whether other dynamical variability is important. We will focus on WACCM, as its dynamics are necessarily consistent with the ozone concentrations.

Figure 9 Correlation coefficient (r) of the interannual variability of total column ozone at every latitude with the total diabatic overturning circulation at every level for (a) ERA-Interim ℳ and SBUV total column ozone, (b) JRA55 ℳ and SBUV total column ozone, (c) MERRA ℳ and SBUV total column ozone, and (d) ozone and ℳ from WACCM.

We see that the high-latitude total column ozone is correlated with the circulation in the lowermost stratosphere, with the correlation explaining up to 25 % of the deseasonalized total column ozone variability in the Northern Hemisphere polar region. The total column ozone in the tropics is strongly anticorrelated with the global diabatic circulation around 500 K. Both of these are qualitatively consistent with transport being the dominant factor driving the relationship between the variability in ozone and in the circulation at these levels. The correlation is strongest in the Southern Hemisphere in the collar region of the polar vortex, around 55^∘ S, and weaker at the pole, while in the Northern Hemisphere, the correlation is stronger at around 70^∘ N. More air is transported by the global diabatic circulation and mixing to the Northern Hemisphere pole than the Southern Hemisphere pole because the Southern Hemisphere polar vortex is a stronger barrier to mixing. The tropical total column ozone is also correlated with the circulation at upper levels above the ozone maximum (800 K). Like water vapor, the correlation is related to the QBO and is strongest at 2-year periods (not shown). Some coherence at higher frequencies can be explained by the anticorrelation of the upper and lower branches of the circulation discussed above. The subtropical total column ozone is anticorrelated with the upper-level circulation strength, with hemispheric asymmetry in which levels relate to subtropical ozone in the different hemispheres. This is consistent with previous results showing that the meridional pattern associated with the QBO at these levels leads to opposite anomalies in the deep tropics and the subtropics (e.g., Randel et al., 1999; Tian et al., 2006).

Figure 10 Correlation coefficient (r) of the interannual variability of local ozone concentration at every latitude and pressure with the global diabatic circulation at (a) 500 K and (b) 1200 K from WACCM. Correlation coefficient (r) of the interannual variability of local temperatures at every latitude and pressure with the global diabatic circulation at (c) 500 K and (d) 1200 K from WACCM. Contours are every 0.1, and correlations are only plotted where they are significant at the 95 % confidence level.

To examine these correlations further, we plot the correlations of the deseasonalized ozone concentrations at each latitude and pressure from WACCM and the deseasonalized global diabatic circulation at two individual levels in Fig. 10. In this way, we try to understand where in the stratosphere the total column ozone correlation patterns are determined. Panel (a) shows the correlation of the local ozone concentration with the global diabatic circulation at 500 K. The strong signal beneath the ozone maximum is consistent with the transport driving the ozone variability – upwelling ozone-poor air from the troposphere and exporting high-ozone tropical air to the midlatitudes and poles in both hemispheres. Panel (b) shows the correlation of the deseasonalized local ozone concentration with the global diabatic circulation at 1200 K. At the Equator at upper levels, the correlation is high, and the strong subtropical signal we see in Fig. 9 is related to the variability of ozone at the uppermost levels and the local ozone concentration on the edge of the tropics in the lower branch. As has been previously reported (Perliski et al., 1989), there is a division between what drives ozone variations in the upper and lower stratosphere. Our results for the global diabatic circulation are consistent with two different processes being responsible for these differing behaviors: near, at, and above the ozone maximum, the ozone distribution is determined by chemistry, while at the lower levels the ozone distribution is determined by transport. Evidence of these two separate processes is discussed below.

The correlation of the upper-level circulation with the lower-level ozone concentrations on the edges of the tropics is consistent with the anticorrelation of the upper and lower branches of the circulation and different characteristics of the transport. In the lower branch, the stratospheric entry latitudes are close to the poleward flanks of the tropics (Birner and Bönisch, 2011), so if the anticorrelation of the upper and lower branches of the circulation is a partitioning between the deep tropical entry latitudes and the more subtropical entry latitudes, the strong upper branch is associated with less upwelling in these flanks and thus less ozone around these turnaround latitudes. This hypothesis for the partitioning of the circulation between upper and lower branches at monthly to seasonal periods and its relationship with trace-gas transport is the subject of further study.

Figure 11 Time series of monthly mean local ozone concentration and the global diabatic circulation from WACCM for the following: (a) 62 hPa, 50^∘ S O₃ in blue and ℳ at 500 K in black; (b) 2 hPa equatorial O₃ in blue, 2 hPa, 50^∘ S O₃ in red (multiplied by −1), and ℳ in black.

Figure 11 shows time series of the local ozone concentrations and total overturning strength based on the maximum correlations shown in Fig. 10. Figure 11a shows the tight coupling between the ozone in the Southern Hemisphere midlatitudes with the global overturning strength at 500 K. Figure 11b shows the very close correlation of the upper-level circulation and the upper-level equatorial ozone and the weaker negative relationship with the upper-level midlatitude ozone. The two time series in (a) and the equatorial ozone and global overturning in (b) are correlated at all timescales, while the anticorrelation between midlatitude ozone and upper-level circulation strength is stronger at short timescales. Ozone variability at upper levels is dominated by photochemical processes (Perliski et al., 1989), resulting in a short chemical lifetime, and this close correlation is due to the relationship of temperature with both ozone and circulation strength. When the circulation is stronger in the tropics at these levels, it is associated with cooling and consequently longer ozone chemical lifetimes. We have therefore plotted the correlation of the temperature with the global diabatic circulation at 500 and 1200 K in Fig. 10c, d. In both (c) and (d), it is evident that at low levels the temperature and ozone respond to the circulation similarly. In (d) in particular, the opposite relationship between circulation and temperature is observed compared to circulation and ozone, which indicates that temperature is driving the chemistry at upper levels. To test this mechanism, we have plotted the natural log of the ozone concentrations against the inverse of temperature at these upper levels and at lower levels, since an exponential dependence on the inverse of temperature is a form that is consistent with the form for many of the reaction rate coefficients for ozone loss processes (Stolarski et al., 2012). These results are shown in Fig. 12. Clearly, the upper level and lower level are behaving differently: because it is dynamically controlled, the lower-level ozone depends as much on latitude as on the inverse of temperature, and the slope is determined by the relative vertical gradients of temperature and ozone. The upper level has little latitudinal dependence and a positive slope, consistent with the chemical control. When the fit is calculated for 45–50^∘ S at 1 hPa, as shown in Fig. 12c, the slope agrees to within error with the slope calculated for the Limb Infrared Monitor of the Stratosphere (LIMS) data used by Stolarski et al. (2012). We have taken the opportunity to show the change in the relationship over time using different colors. Calculating the fit for just the earlier years results in a higher value for the “initial” ozone concentration with a slope that is the same to within error. While we do not investigate the cause for this change here, we hypothesize that it is related to the higher mean ozone concentrations being advected to this region during the initial period of the ozone hole.

Stratospheric transport timescales for even the lower branch are around half a year to a year (Orbe et al., 2014), so instantaneous correlation plots, as in Fig. 9, might seem to be less relevant. The global diabatic circulation necessarily integrates over that transit time, however, as it accounts for variability in both the upwelling region and in the downwelling region simultaneously. Therefore, we do not perform lagged regressions to attempt to understand causality. Rather, we suggest that the use of frequency-dependent correlations, which will have a corresponding phase lag (e.g., Swanson, 2000), would be necessary to look at the causal relationships. However, as we can see from Fig. 11, the correlations are in phase at monthly timescales, so higher-frequency records (minimum daily) will be necessary to diagnose the phase (and thus the implied causality) in these relationships.

Figure 12 The natural log of the ozone volume mixing ratio (in ppm) and the colocated values of 1000 ∕ T for different levels and latitudes in the WACCM model for 1980–2014. (a) 53 hPa and (b) 2 hPa; both are for all latitudes equatorward of 52^∘. (c) 1 hPa for 45–50^∘ S only. In (a) and (b), different colored circles are different latitudes, with the Northern Hemisphere being yellow and the Southern Hemisphere blue. In (c), different colors show different years. In (b) and (c) the best-fit lines are also plotted in black with correlation coefficients of r=0.87 and r=0.96, respectively.

The correlations between ozone and the global diabatic circulation have a resemblance to the pattern of the ozone response to the QBO, but in contrast to the water vapor tape recorder examined above, coherence at higher frequencies suggests that other processes play a role. We have demonstrated the close dependence of ozone variability on the global diabatic circulation variability with ozone data and with reanalysis and model data. The total column ozone at the poles and in the tropics is correlated with transport by the global diabatic circulation in the lower stratosphere. The results of the investigation of the correlation of ozone with the global diabatic circulation have demonstrated consistency with our understanding of the roles of circulation and chemistry, and we suggest that the global diabatic circulation can be adopted in this context with little to no change in interpretation compared to ${\overline{w}}^{*}$.

7 Discussion and conclusions

In this study, we have compared the global diabatic circulation to the more typically used metrics for the strength of the BDC and to tracers. In particular, we have examined the residual circulation vertical velocity, the water vapor tape recorder, and total column ozone concentrations.

We find that the three common methods for quantifying the BDC strength from models and reanalysis data products have somewhat different deseasonalized variability, especially in the lower stratosphere. We also find that the choice of averaging latitudes – whether fixed tropics (30^∘ S–30^∘ N) or turnaround latitudes – has an effect on the deseasonalized variability that depends on the method. These methods also result in different vertical structures; the calculation based on the principle of downward control has a much broader vertical autocorrelation than the calculation from diabatic heating rates. The direct method is somewhere in between, and when the turnaround latitudes are used, it becomes very similar to the calculation from diabatic heating rates. Thus, if the separate evolution of the upper and lower branches of the circulation is of interest, the most appropriate metric is one that uses diabatic heating rates or the direct method with turnaround latitudes. In the model, the relationship of the different TEM calculation methods with fixed tropics are nearly one to one above 70 hPa. For the reanalysis products, however, the differences between calculations of the TEM ${\overline{w}}^{*}$ are quite distinct, especially at lower levels where they are often analyzed. The comparison between the TEM ${\overline{w}}_{\mathrm{Q}}^{*}$ with the diabatic overturning circulation is as favorable as the comparison between the TEM ${\overline{w}}_{\mathrm{Q}}^{*}$ and the ${\overline{w}}^{*}$ from the residual circulation method. In general, consistency between methods is better higher up in the stratosphere, while beneath 70 hPa, the differences between the methods are substantial. These results suggest that the method of calculation could significantly affect comparisons between the residual circulation from reanalysis and any other observed stratospheric variable.

Like the thermodynamically constrained ${\overline{w}}_{\mathrm{Q}}^{*}$, which the global diabatic circulation so closely resembles, the global diabatic circulation requires only monthly mean heating rates, temperatures, and pressures. Its calculation is simpler than that of $\overline{w}{*}_{\mathrm{Q}}$, which requires some assumption about how to enforce mass conservation (Abalos et al., 2012) and which can have complications with convergence when the iterative solving method converges but then occasionally proceeds to diverge after additional iterations. Eddy terms in the thermodynamic equation are neglected in the calculation, which may be a reason for these convergence difficulties. The global diabatic circulation also has an interesting property: the lower and upper branches of the circulation are anticorrelated so that when the lower branch is stronger, less air is flowing through the upper branch. This is even more curious when one takes into account the fact that the vertical velocity in the lower stratosphere is the sum of both branches. This pattern might be expected with the QBO, but as the coherence is not just at QBO frequencies, an additional mechanism is necessary. One explanation is that this could be due to a change in the index of refraction when there is more total wave activity that causes higher-amplitude planetary-scale waves to break lower in the stratosphere; we have yet to test this mechanism. Another possibility is that the meridional location of the wave breaking changes such that when the lower branch is stronger, less wave activity can propagate up into the upper stratosphere. Alternatively, there may be an interaction between planetary and gravity waves. The anticorrelation is consistent with the conclusions of both Ray et al. (2010) and Stiller et al. (2012), who concluded based on observations that the trends in data were best explained by a strengthening in the lower branch of the circulation and a weakening in the middle and upper stratosphere. The ERA-Interim trends in the global diabatic overturning circulation are consistent with this picture, although the upper-level trends are problematic because of the changes in observing systems. The other two reanalyses do not agree.

The global diabatic circulation is correlated with the water vapor tape recorder vertical velocity, especially at intraseasonal and longer timescales. Perhaps unsurprisingly, given its close theoretical relationship with stratospheric tracers, the global diabatic circulation is a predictor for the water vapor tape recorder. However, the overall weakness of the correlation, which explains at most < 40 % of the variability even when both metrics are derived from a model, suggests the inadequacy of using the water vapor tape recorder as a lone observational record of the changing stratospheric circulation. Rather, the water vapor signal should be compared to water vapor in models in order to assess the combined effect of diabatic heating, diabatic diffusion, and adiabatic mixing.

We analyze the impact of the global diabatic circulation on total column ozone using satellite data and the three reanalyses, including examining the dependence of the total column ozone on the circulation through different vertical levels. When we find consistent behavior amongst the three reanalyses, we explore the mechanism using a model that shows the same behavior. We find that tropical ozone is most correlated with the overturning at 500–550 K, Southern Hemisphere ozone is sensitive to the global diabatic circulation at around 480 K, and Northern Hemisphere ozone is most sensitive to the global diabatic circulation between 400 and 450 K. The subtropics are most sensitive to the mid-level circulation at 800–1000 K, related to the dominant role of chemistry at upper levels. Generally, the patterns associated with the ozone correlation with the global diabatic circulation are consistent, with much of this relationship being related to the QBO.

Based on its close relationship with one of the common metrics for the BDC, the ease of calculation, the demonstrated impact on ozone and water vapor, and the constraints provided by tracer observations, we present the global diabatic circulation as a metric for the BDC that should be newly considered. Before the community settled on ${\overline{w}}^{*}$, the global diabatic circulation was used (Pyle and Rogers, 1980; Rosenfield et al., 1987). Some intuition for the behavior of ${\overline{w}}^{*}$ exists, but both Abalos et al. (2015) and this work have demonstrated that the various methods of calculation are not equivalent, especially for reanalyses. Thus, although some variety of TEM ${\overline{w}}^{*}$ is the most common metric at present, its calculation is not held in common amongst different studies. In order to understand models and reanalyses, consistency is critical. For the purposes of reanalysis evaluation, therefore, we advocate using the global diabatic circulation along with a version of the quasi-geostrophic TEM ${\overline{w}}^{*}$ with fixed tropical averaging latitudes (as the turnaround latitudes for the reanalyses are not always well defined, which limits the vertical extent of comparisons). These two metrics rely on different assumptions, and the heating rates from reanalysis might be suspect. For the purposes of model evaluation, the global diabatic circulation should be sufficient. The latitudinal structure of the circulation cannot be examined using the global diabatic circulation, however, so the vertical velocity should be used when meridional structure is of interest.

Apart from the brief analysis with ozone, this paper does not directly address causality. It is an investigation of different metrics for the circulation from an empirical perspective, revealing significant differences in the behavior of these metrics that raise questions about their interchangeability, especially for reanalyses. The inconsistencies reveal the extent to which the reanalyses' momentum and energy budgets are not internally consistent. At upper levels, the different vertical velocities are all nearly equivalent, but at lower levels, especially beneath 70 hPa, the differences are substantial. In particular, using the momentum-based calculation for the residual circulation vertical velocity will mask variability that is not coincident between the upper and lower branches, while the global diabatic circulation emphasizes the difference between the upper and lower branch. This work serves as motivation for additional process-based and theoretical studies that address the causes of these differences between residual circulation metrics and between tracers and the residual circulation.