Inter-comparison of stratospheric mean-meridional circulation and eddy mixing among six reanalysis data sets

. The stratospheric mean-meridional circulation (MMC) and eddy mixing are compared among six meteorological reanalysis data sets: NCEP-NCAR, NCEP-CFSR, ERA-40, ERA-Interim, JRA-25, and JRA-55 for the period 1979–2012. The reanalysis data sets produced using advanced systems (i.e., NCEP-CFSR, ERA-Interim, and JRA-55) generally reveal a weaker MMC in the Northern Hemi-sphere (NH) compared with those produced using older systems (i.e., NCEP/NCAR, ERA-40, and JRA-25). The mean mixing strength differs largely among the data products. In the NH lower stratosphere, the contribution of planetary-scale mixing is larger in the new data sets than in the old data sets, whereas that of small-scale


Introduction
The Brewer-Dobson circulation (BDC), which was discovered by Brewer (1949), Dobson at al. (1929), andDobson (1956), consists of the mean-meridional circulation (MMC) and eddy mixing in the stratosphere.The stratospheric MMC is composed of ascending motions in the tropics, poleward motions toward mid and high latitudes, and descending motions at high latitudes.Planetary waves that propagate from the troposphere break and cause eddy mixing primarily in the stratospheric surf zone surrounding the polar vortex in the winter hemisphere (McIntyre and Palmer, 1983).Because the BDC motion is too slow to measure directly, the detailed structure and long-term variations of the BDC are not well understood.
Stratospheric age-of-air (AoA) (Waugh and Hall, 2002) derived from observation data of long-lived chemical compounds is frequently used as a surrogate for the combined effects of the MMC and eddy mixing in order to investigate the structure and long-term variations in the BDC.Observational studies have found a positive AoA trend in the Northern Hemisphere (NH) mid-latitudes in the middle to Published by Copernicus Publications on behalf of the European Geosciences Union.
upper stratosphere based on balloon measurements of SF 6 and CO 2 for 1975-2005(Engel et al., 2009)), MIPAS SF 6 measurements for 2002-2010 (Stiller et al., 2012), balloonbased measurements of SF 6 and CO 2 for 1975-2012 (Ray et al., 2014), and a merged long-term satellite data record of water vapour for 1986-2010 (Hegglin et al., 2014).In the Southern Hemisphere (SH) mid-latitude lower stratosphere, Stiller et al. (2012) found in contrast a negative AoA trend, a result confirmed by Hegglin et al. (2014).The study by Hegglin et al. (2014) also finds a negative AoA trend in the NH lower stratosphere, in contrast to Stiller et al. (2012), with the difference likely explainable by the different time periods considered.However, these observed trends are not fully consistent with simulated results using general circulation models (GCMs), which show instead an acceleration in the mean BDC throughout the stratosphere (e.g., Butchart et al., 2011).The acceleration of the mean BDC strength in the model simulation is associated with enhanced wave driving (Garcia and Randel, 2008) and shifting critical levels for wave breaking (Shepherd and McLandress, 2011).The reasons for the inconsistency between the measurements and models, especially in the middle and upper stratosphere, have not yet been investigated.
MMC and eddy mixing in the BDC are intrinsically linked (Dunkerton, 1978) and hence both play important roles in determining distributions of long-lived chemical species and AoA in the stratosphere.Interpretation of AoA variations therefore needs to take into account changes in both MMC and eddy mixing.Several recent studies have now quantified the effects of eddy mixing on AoA variations in more detail (Ray et al., 2010;Garny et al., 2014;Ploeger et al., 2015a, b).Garny et al. (2014) stated that eddy mixing causes recirculation of air in the stratosphere, and acts to increase the mean value of AoA.Mixing inside the surf zone modifies the latitudinal distribution of AoA and chemical species, whereas mixing across the subtropical transport barrier is suggested to be important for the mean AoA value above the mixing level (Garny et al., 2014;Ploeger et al., 2015a).
Meteorological reanalyses provide a realistic meteorological field by combining model information with actual observations, and have the potential to provide an alternative tool to study the BDC and AoA variations.Iwasaki et al. (2009) compared four reanalysis data sets (NCEP-NCAR, NCEP-DOE, ERA-40, and JRA-25) and found large differences in their representation of the BDC.Wright and Fueglistaler (2013) showed large differences in the simulated diabatic heat budget in the tropical upper troposphere and lower stratosphere in five reanalysis models (NCEP-NCAR, NCEP-CFSR, JRA-25, ERA-Interim, and MERRA), with substantial implications for representation of transport and mixing.Recently, Abalos et al. (2015) compared the MMC in three reanalyses (ERA-Interim, JRA-55, and MERRA) using three different estimates: from the transformed Eulerian mean (TEM, Andrews and McIntyre, 1976) residual circulation and based on momentum and thermodynamic balances.
They showed a relatively large spread (around 40 %) among the estimates of the magnitude of tropical upwelling.Monge-Sanz et al. (2007, 2012) highlighted the fact that the representation of the BDC (including eddy mixing) has become much more realistic in ERA-Interim than in ERA-40, thanks to large investments made in improving the reanalysis product.They showed that AoA derived using ERA-Interim displays an increasing trend in the NH mid-latitude stratosphere above 25 km in 1989-2010, consistent with the findings by Hegglin et al. (2014) for 1986-2010. Ploeger et al. (2015a) quantified the effects of mixing and residual circulation on mean age trends using ERA-Interim, and showed the importance of eddy mixing in the increasing AoA trend in the NH middle stratosphere in 2002-2012.Diallo et al. (2012) showed that the AoA trend derived using ERA-Interim in the middle stratosphere is positive over the 1989-2010 period.In order to gain confidence in the ERA-Interim results, it is now important to assess their consistency with long-term variations in the BDC obtained from other reanalysis products.
Differences in the forecast model, assimilated measurements, and data assimilation technique used for producing reanalysis data sets can lead to differences in their representation of the BDC.Model simulations without any assimilation produce meteorological fields that follow the dynamical and thermodynamic balance of the forecast model.Data assimilation analysis increments, introduced by using conventional data assimilation techniques such as the threedimensional variational (3D-VAR) one, may upset this balance and degrade the expression of momentum budget and wave structures.This is because they introduce an additional force, without maintaining physical balance, as a result of its isotropic and instantaneous analysis increment.In the 3D-VAR analysis, mean ascending motions in the tropics and mixing in the subtropics in the stratosphere were found to be excessively strong (Schoeberl et al., 2003;Tan et al., 2004;Scheele et al., 2005).Advanced data assimilation techniques such as the four-dimensional variational method (4D-VAR) are capable of assimilating observations at the exact time while maintaining the dynamical balance because of the use of flow-dependent analysis, which are expected to improve the representation of both MMC and eddy mixing.
In this paper, we compare MMC and eddy mixing in the stratosphere for six reanalysis data sets: NCEP-NCAR, NCEP-CFSR, ERA-40, ERA-Interim, JRA-25, and JRA-55.The analysis is conducted for the 34 years from 1979 to 2012 based on mass-weighted isentropic zonal means that allow accurate analysis of Lagrangian-mean motions and eddy mixing.Based on the comparison of the mean and eddy components in the BDC, we discuss whether any of the reanalysis data have the potential to reveal useful information on longterm AoA variations.

Data
The six reanalysis data sets used in our comparison can be described as follows: (1) NCEP-NCAR -the National Centers for Environmental Prediction (NCEP)-National Center for Atmospheric Research (NCAR) reanalysis product (Kalnay et al., 1996), with a model grid resolution of T62L28 produced using a 3D-VAR technique; (2) NCEP-CFSRthe NCEP Climate Forecast System Reanalysis (Saha et al., 2010), with a model grid resolution of T382L64 produced using a 3D-VAR technique; (3) ERA-40 -the 40-year ECMWF Re-Analysis (Simmons and Gibson, 2000), with a model grid resolution of T159L60 produced using a 3D-VAR technique; (4) ERA-Interim -a continuously updated reanalysis since 1979 (Simmons et al., 2007), with a model grid resolution of T225L60 produced using a 4D-VAR technique; (5) JRA-25 -the Japanese 25-year reanalysis product (Onogi et al., 2007), with a model grid resolution of T106L40 provided using a 3D-VAR technique; and (6) JRA-55 -the Japanese 55-year reanalysis product (Kobayashi et al., 2015), with a model grid resolution of T319L60 provided using a 4D-VAR technique.
We here classify NCEP-NCAR, JRA-25, and ERA-40 as old data sets, and NCEP-CFSR, JRA-55, and ERA-Interim as new data sets because of improvements in these latter ones made by using updated forecast models, updated bias correction algorithms, and advanced data assimilation analysis.We use reanalyses results for the 34 years after 1979, when satellite measurements were assimilated into the reanalysis.For ERA-40, the mean state and linear trend are estimated for the 24 years from 1979 to 2002, since the data are not available after 2002.

Analysis framework
The analysis of MMC and eddy mixing is based on massweighted isentropic zonal means (hereafter referred to as MIM analysis; Iwasaki, 1989Iwasaki, , 1992;;Miyazaki and Iwasaki, 2005).The MIM zonal mean is defined as where φ is the latitude, θ is the potential temperature, t is the time, λ is the longitude, and p is the pressure.The asterisks and overbars represent mass-weighting and isentropic zonal means, respectively.Eddies are defined as departures from the mass-weighted zonal means, Their correlations are given by We use isentropic zonal mean pressure for the vertical coordinate, p † ≡ p. (4) Unlike other isentropic coordinate analyses (Andrews, 1983;Tung, 1982Tung, , 1986)), in the MIM analysis, the massweighting is considered not only for meridional circulation, but also for other variables such as zonal wind following Johnson (1980).The MIM analysis thus expresses the conservative nature of momentum, heat, and minor constituents, including the exact lower boundary conditions and non-geostrophic effects (Iwasaki, 1989(Iwasaki, , 1992;;Miyazaki and Iwasaki, 2008).The MIM analysis also exactly specifies the eddy diabatic and adiabatic transport terms (Miyazaki and Iwasaki, 2005).
The TEM (Andrews and McIntyre, 1976) provides a useful framework for understanding mean and eddy transports; however, the estimation of the transport fluxes is limited and complicated.Randel et al. (1994), Strahan et al. (1996), and Abalos et al. (2013) estimated eddy transport terms based on eddy flux vectors for small amplitude eddies following Andrews et al. (1987), while some studies estimated this term as residuals considering the uncertainty and difficulty in the eddy transport term estimations (e.g., Randel et al., 1998).Miyazaki et al. (2008) found a significant (> 30 %) difference in the mean vertical velocity around the Antarctic polar vortex between the TEM residual vertical velocity and the MIM MMC analysis, which can be attributed to the assumptions applied for the TEM Stokes corrections.Most of the disadvantages of conventional analysis methods such as TEM (e.g., complicated and inaccurate representation of transport by both mean and eddy motions, lower boundary conditions, and mass conservations) can be avoided using the MIM analysis (Iwasaki, 1989;Tanaka et al., 2004;Miyazaki and Iwasaki, 2005).
We here focus on two altitudes: 440 K (at approximately 90-80 hPa) and 560 K (40-30 hPa) as representatives of the shallow and deep branches of the BDC, respectively.Birner and Bönisch (2011) reported that the shallow branch extends to about 50 hPa, and the deep branch is located above that altitude.The analysis results are presented on isentropic coordinates for diagnosing adiabatic and diabatic transport components.Note that if the potential temperature at a constant pressure changes with time, there would not be exact agreement between the estimated trends in pressure and isentropic coordinates.For example, in the last 30 years (2008-2012 mean minus 1979-1983 mean), potential temperature at 70 hPa decreased by about 2.5 K at low and mid latitudes in both hemispheres in ERA-Interim.Nevertheless, the general structure of the linear trend was similar between the two coordinates (and hence will not be shown here).The longterm linear trend is estimated based on the least-squares fitting.The statistical significance is determined for the 95 % confidence level using the Mann-Kendall test (Mann, 1945;Kendall, 1975

Mean-meridional circulation (MMC)
The mass stream function x in the MIM analysis is calculated from integrating meridional velocity with respect to p † : where a is the Earth's radius and v is the meridional wind.v * is calculated from meridional wind data with consideration of the mass-weighted isentropic zonal means.The diagnostic form of the zonal mean continuity equation without an eddy term in the MIM analysis confirms that the mean-meridional circulation can be expressed by the nondivergent mass stream function (Iwasaki, 1989).The mean vertical velocity w * † is obtained from the mass stream function x as follows: where ρ 0 is the reference atmospheric density.The mean vertical velocity w * † can be related to the diabatic heating θ * (Eq.2.7 in Tanaka et al., 2004).

Eddy mixing
By assuming a flux-gradient linear relationship, the diffusion coefficient provides a measure of eddy mixing.The isentropic diffusion coefficient K yy can be derived from the eddy meridional flux and meridional potential vorticity (PV) gradient on isentropic surfaces (Tung, 1986;Newman et al., 1988;Bartels et al., 1998;Miyazaki and Iwasaki, 2005;Miyazaki et al., 2010b) by neglecting the influence of slant diffusion: where q is the PV and []l denotes the time average.Under frictionless and adiabatic conditions, the PV acts as an atmospheric passive tracer (Hoskins et al., 1985).Miyazaki et al. (2010a) demonstrated that the diabatic source-sink effect on the PV budget is much smaller than the transport effects in the subtropical and extratropical stratosphere.In the K yy estimation, steady conservative wave motions projected onto a meridional plane can cause apparent diffusion in addition to the true diffusion caused by dissipative wave motions, which may lead to significant differences between the estimated K yy and the true eddy mixing.In order to reduce this effect in our estimates, a time-averaging window t is applied to the eddy PV flux, as in Miyazaki and Iwasaki (2005) and Miyazaki et al. (2010b).As a result, we expect that the estimated K yy will provide information on eddy mixing characteristics similar to estimates of the effective diffusivity (Nakamura, 1996;Haynes and Shuckburgh, 2000).
The absolute value of K yy is influenced by the choice of l (set to one month in this study).We confirmed that, in all data sets, the estimated K yy becomes smaller by increasing l from 6 h to 10 days, and it becomes nearly constant when l is set to be more than 10 days.For instance, K yy (l > 10 days) is approximately 15-20 % smaller than that of K yy (l = 6 h) in the NH middle stratospheric surf zone in December to February (DJF), for the case of ERA-Interim.These constant values are considered to represent the diffusion coefficient due to the true diffusions.
We should note that, even after eliminating the influence of apparent diffusion in K yy estimates, there are limitations in elucidating eddy mixing from these estimates.As discussed by Nakamura (2008), whilst a part of the Eulerian eddy diffusivity (e.g., K yy ) can be attributed to instantaneous, irreversible mixing in a way similar to effective diffusivity, the Eulerian eddy diffusivity and eddy diffusivity are fundamentally different, both qualitatively and quantitatively.This is because of difficulties associated with representing eddy advective transport in the Eulerian formulation.In addition, the results of the K yy analysis and related variables are presented here in the geometric latitude coordinate system (where a zonal average is taken for air parcels with different PVs), whereas effective diffusivity is presented in the equivalent latitude (EL) coordinate system (based on the latitude circle that encloses the same area as the PV contour).Latitudinal variations of zonal-mean eddy mixing and associated fields (e.g., strong eddy mixing outside the polar vortex) are more clearly presented in the EL coordinate system.

The relative importance of mean and eddy transports
The zonal mean equation in the MIM analysis can accurately separate meridional transport into mean transport by Lagrangian-mean circulation and eddy (diffusion) transport.
In the MIM analysis, the mean and eddy PV fluxes are defined as follows: where v * is calculated from meridional wind data (cf.Sect.2.2.1).The mean transport fluxes are parallel to isopleths of the mass stream function, whereas the eddy transport fluxes are parallel to the isentropes under diabatic conditions.Only the meridional components are analyzed in this study.
In order to evaluate the relative importance of eddy transport in each reanalysis product, we estimate the ratio of eddy and total meridional transport fluxes of PV as follows:  If the transport ratio is larger (smaller) than 0.5, the eddy transport (mean transport) dominates the meridional transport on isentropic surfaces.As in the K yy estimation (cf.Eq. 7), a time average (l = 1 month) was applied to the mean and eddy transport fluxes.

Eddy mixing
Figure 1 shows the seasonal mean isentropic diffusion coefficient K yy (l = 1 month) (hereafter referred to as K yy ).Large K yy values reveal strong isentropic mixing in the stratospheric surf zone in both hemispheres, which is stronger in the NH than in the SH.The hemispheric asymmetry can be attributed to the differences in planetary wave activity (e.g., Shepherd, 2000).The small K yy equatorward of about 30 • is indicative of a barrier to horizontal transport between the tropics and mid-latitudes in both hemispheres (Plumb, 1996).At SH high latitudes, K yy becomes large after the break-up of the Antarctic polar vortex (figure not shown).The eddy mixing is strongly suppressed near the core of the subtropical jet stream, but shows clear maxima at its upper flank (equatorwards) and also at its lower flank (poleward, during DJF).These general characteristics are commonly found in all the data sets, and are generally consistent with the analysis of ef-fective diffusivity (e.g., Haynes and Shuckburgh, 2000).It is noted that the counter-gradient transport is found in the negative region of K yy (shaded in white) in the summer hemisphere, which is associated with almost flat PV gradients.
The mean value and linear trend of K yy are estimated in the surf zone (40-60 • ) and subtropics (15-25 • ) around the overturning latitude of the MMC (Table 1 and Fig. 2).Isentropic mixing in the surf zone influences the latitudinal gradient of tracers, whereas the mixing across the subtropics is considered to be important for the stratospheric mean AoA because it recirculates old air from the extratropics to the tropics (Garny et al., 2014;Ploeger et al., 2015a).The anomalously low value in NCEP-NCAR in the middle stratosphere is associated with the low top height of the reanalysis product.
K yy shows an increasing trend in the NH middle stratosphere (i.e., 560 K) surf zone in all the data sets over the time period 1979-2012 (+0.5 to +5.2 % decade −1 ), with relatively weak trends in the newer data sets (+0.5 to +2.2 % decade −1 ).The trend in the NH at 440 K is positive in ERA-Interim and JRA-55 in 1979-2012and also in ERA-40 in 1979-2001.In the SH surf zone, the K yy shows an increasing trend from the lower to middle stratosphere in all the data sets except for NCEP-NCAR and ERA-40.at 560 K (+13.9 % decade −1 ).The intensified surf zone mixing could be associated with changes in the critical level.Climate model simulations demonstrate that long-term changes in zonal wind such as the shift of zero wind line in response to climate change can enhance the upward propagation of westward-propagating waves (Kawatani et al., 2011;Shepherd and McLandress, 2011).Further investigations are required to comprehend the relationship between changes in the critical level, wave forcing, and mixing strength in the reanalysis products.The standard deviation in the middle stratosphere surf zone is smaller in the new data sets than in the old data sets in both hemispheres by about 30-60 % in both hemispheres.
The mean K yy value differs largely among the data products in the NH subtropical lower stratosphere (0.56-0.81).The large difference in the subtropical mixing among the data sets could be associated with different representations of the Quasi-Biennial Ocillation (QBO) in different reanalysis products (Pawson and Fiorino, 1998;Kim and Chun, 2015;Kawatani et al., 2016).In the SH subtropical lower stratosphere, the K yy trend is positive only in ERA-Interim (+1.2 % decade −1 ) and NCEP-NCAR (+8.0 % decade −1 ), and shows a large negative value in NCEP-CFSR (−10.0 % decade −1 ) and ERA-40 (−11.2 % decade −1 ).Because of the large interannual variations, the estimated trends are not statistically significant for most cases in the subtropics and the mid-latitudes.
The K yy trend in the NH middle stratosphere surf zone is nearly zero or positive (−1.8   ERA-40, then changes to large negative values (−5.7 to −31.5 % decade −1 ) in the last 12 years in all the data sets (Table 2).The negative trend in the latter period is larger in the old data sets (−23.5 to −31.5 % decade −1 ) than in the new data sets (−5.7 to −16.8 % decade −1 ).These decadal scale changes in the mixing trend seem to be consistent with those in the tropical upward mass flux and MMC in the NH (cf.Sect.3.2).This suggests that variations in wave forcing lead to decadal scale changes in both the mean and eddy transports as is expected, given that these two features are intrinsically connected with each other (cf.Sect.4.1).In the SH middle stratosphere surf zone, the K yy trend is positive (+4.4 to +14.0 % decade −1 , except for ERA-40) in the first 22 years and becomes negative (−24.4 to −40.4 % decade −1 ) in the last 12 years.

Mean-meridional circulation (MMC)
Figure 3 compares the meridional cross section of the mass stream function averaged over DJF and June to August (JJA) in 1979-2012.Note that ERA-40 has been averaged over the shorter time period 1979-2002 due to limited availability of the data set.Comparison of ERA-40 and ERA-Interim averaged over the same (shorter) time period, however, indicates that the results are comparable.The general structure of the MMC, such as the poleward motion from low latitudes of the summer hemisphere to mid and high latitudes of the winter hemisphere, and the descending motions at high latitudes of winter hemispheres, is commonly found in all the data sets.However, details regarding the structure, intensity, and trend of the MMC differ among the data.In NCEP-NCAR,  the MMC is noisy and unrealistically distorted in the middle stratosphere (altitudes above about 550 K), and also in the SH lower stratosphere during JJA.In ERA-40, it is relatively strong throughout the stratosphere in both hemispheres, as already pointed out by Wohltmann and Rex (2008).Similarly, JRA-25 exhibits an MMC that is somewhat stronger in the middle to upper stratosphere than that found in the other data sets.In NCEP-CFSR, the tropical mean upward motion is distorted in the lower stratosphere, so that mean air trajectories do not ascent in a straight line.Significant diversity of the tropical mean upward motions (not shown) may be related to differences in simulated diabatic heating rates in the forecast models (Fueglistaler et al., 2009).Figure 4 compares the seasonal variation of the mass stream functions at 560 K.The strong mean poleward motions are present from low latitudes to mid and high latitudes from autumn to spring.They are relatively suppressed at the polar vortex edge.In NCEP-NCAR, the cross-equatorial mean-meridional flow is relatively weak throughout the year.Also, the strong poleward flow at mid-latitudes does not exhibit the seasonality apparent in the other reanalyses.In ERA-40 and JRA-25, the mean poleward flow from the subtropics to the mid-latitudes is relatively strong in both hemispheres throughout the year.
The MMC long-term trend differs greatly among the reanalysis data sets (Fig. 5).In NCEP-NCAR, the spatial structure of the linear trend slope is very noisy, likely reflecting inhomogeneities in the reanalysis product.In NCEP-CFSR, the MMC shows a strengthening trend from the tropics to the mid-latitudes in the lower and middle stratosphere in both winter hemispheres.JRA-25 reveals a strengthening trend in the lower and middle stratosphere equatorward of 60 • N in the NH during winter and throughout the lower stratosphere in the SH in both seasons, whereas it shows a weakening trend above 500 K in the SH in DJF.JRA-55 shows a weak positive trend in the lower stratosphere in winter.ERA-40 reveals a strengthening trend throughout the lower and middle stratosphere over 1979-2002 in the SH in JJA, whereas the trend pattern is noisy in the NH in both seasons.Only in ERA-Interim, the MMC in the NH middle stratosphere (or above about 480 K) tends to weaken during winter.A negative trend is also found from the tropics to the mid-latitudes in both hemispheres in JJA in the middle stratosphere in ERA-Interim.The BDC shallow branch mostly shows a strengthening trend in ERA-Interim as in other data sets.Abalos et al. (2015) revealed that the acceleration in the MMC is a qualitatively robust result across different estimates (from the TEM residual circulations and based on momentum and thermodynamic balances).Obvious structural changes (i.e., that the shallow and deep branches changed differently) in the wintertime MMC can be found only in ERA-Interim and to some extent also in JRA-25.There are only a few statistically significant regions in the BDC trends.The low statistical significance could be partly attributed to the fact that large variances associated with ENSO and QBO, as pointed out by Abalos et al. (2015), are not removed in our trend estimates.

Wave decomposition
To gain a better understanding of the relative contributions of the wave phenomena on various scales to isentropic mixing, the zonal eddy PV flux is decomposed into individual zonal wavenumber components s using fast Fourier transform (FFT) and then normalized by the background meridional PV gradient as follows: The eddy PV flux is separated into three groups: planetary waves (associated with the zonal wavenumbers 1-3), synoptic-scale waves (zonal wavenumbers 4-7), and smallscale waves (zonal wavenumbers 8 and higher).For this analysis, l is set to 1 day to include short-term variations in eddy mixing, although this may cause apparent diffusion in the K yy estimates.
The relative contribution of these groups is summarized in Table 6.The relative contribution of the planetary-scale mixing to the eddy flux is larger in the new data sets than in the old data sets for all the data sets at 440 K in the NH midlatitudes, whereas that of zonal disturbances at both synoptic and small scales are somewhat smaller.With conventional data assimilation techniques such as 3D-VAR, physical consistency cannot be maintained during the data assimilation analysis because of lack of flow-dependent background error information, and this probably together with the lower forecast model resolution may cause spurious disturbances and excessive mixing, especially at small scales in the old data sets.For planetary-scale and synoptic-scale waves, forecast model configurations may also be important for the obtained differences (e.g., Biagio et al., 2014).
Each wave group reveals different long-term variations.The planetary-scale disturbance shows increasing trends in the NH surf zone at both 440 K (+0.6 to +4.2 % decade −1 in 1979-2012 and +4.8 % decade −1 in 1979-2001 for ERA-40) and 560 K (+0.6 to +11.5 % decade −1 in 1979-2012 and +11.1 % decade −1 in 1979-2001 for ERA-40).At 560 K, all the groups reveal positive trends in all the data sets, with smaller linear trends in the new data sets in most cases, consistent with the result of the K yy trend analysis except for ERA (cf.Sect.3.1), although the estimated trends are not statistically significant.In the subtropical lower stratosphere (not shown in table), as commonly found in the NH surf zone, the planetary-scale disturbance is stronger, and the smallscale disturbance is weaker in the new data sets than in the old data sets.

Relative importance of eddy meridional transport
Atmospheric wave breaking drives both the mean and eddy transports, but in a different manner; it drives the MMC below the breaking level through downward control (Haynes et al., 1991), whereas it causes eddy mixing at the level of wave breaking (e.g., Waugh et al., 1994).Therefore, changes in MMC and eddy mixing do not necessarily occur in the same way in response to changes in atmospheric wave forcing.These wave-induced processes may be represented differently among the reanalysis products.
We compare the relative ratio of the eddy and total meridional PV fluxes (Eq.10) to investigate the relative importance of these two transport processes.As shown in Fig. 8, the mean-meridional transport is dominant equatorward of about 40 • throughout the year, and in the extratropics during the summer.The eddy transport is dominant in the extratropics of the winter hemisphere.Air masses entering the stratosphere are thus dominantly transported by the meanmeridional transport from low to mid-latitudes and then transported mainly by the eddy mixing from the mid-to high latitudes in the winter hemisphere.The strong westerly coincides with the large contribution of the eddy transport.These characteristics are found in a consistent way in all the data sets.The mean vertical transport may be pronounced in some regions (e.g., inside the polar vortex), but the relative importance of vertical transport is not evaluated in our analysis.
The 34-year mean value of the transport ratio is larger (i.e., the contribution of eddy transport is larger) in the new data sets than in the old data sets both in the extratropics (by 5-16 % at 440 K and by 11-27 % at 560 K) and subtropics (by 5-14 %) in the NH (Table 7 and Fig. 9).In the NH surf zone at 440 hPa, the trend in the transport ratio is different among the data sets.At 560 hPa, the transport ratio shows a negative trend (i.e., the contribution of the mean transport becomes more important over time) in all the data sets.In the SH surf zone, the trend in the transport ratio varies among the data set in 1979-2012 (−10.8 to +0.4 % decade −1 at 440 K and −3.2 to +2.6 % decade −1 at 560 K, except ERA-40).The standard deviation is generally smaller in the new data sets than in the old data sets in both hemispheres.
In the NH subtropics at 440 K, the large mean value of the transport ratio in the new data sets when compared to the old data sets suggests that eddy-induced recirculation of old air masses from the extratropics to the tropics is more effective than the mean poleward transport of fresh air masses from the tropics into the extratropics in the new data sets, which could have led to older AoA as derived for example in the study by Monge-Sanz et al. (2012), as discussed further in Sect.4.2.The transport ratio in the NH subtropics shows a positive trend in most data sets except for NCEP-NCAR and NCEP-CFSR.In the SH subtropics, the trend is negative Table 7. Same as in Table 1, but for the relative importance of the eddy transport to the total meridional transport averaged over 40-60 • N at 440 and 560 K and 15-25 • N at 440 K in DJF, and over 40-60 • S at 440 and 560 K and 15S-25 • S at 440 K in JJA, and its linear trend slope (in % decade −1 ).The linear trend with a confidence level greater than 95 % is shown in bold.The transport ratio in the NH subtropics reveals a larger positive trend in the last 12 years (+3.5 to +6.8 % decade −1 ) than in the first 22 years (−2.8 to +1.8 % decade −1 ) in ERA-Interim and NCEP-NCAR (Table 8).This positive trend implies that the eddy transport has become more important over time with this tendency becoming even more important during the later time period considered.On the other hand, the trend changes from positive to negative between the two periods for NCEP-CFSR, JRA-25, and JRA-55.In the SH subtropics, the mean transport becomes more important in all data sets except ERA-Interim in the earlier time period.The trend in the later time period is strongly positive in ERA-Interim (+8.3 % decade −1 ) and JRA-25 (+21.9 % decade −1 ).

Dynamical consistency
The relationship between the Eliassen-Palm (E-P) flux divergence and the MMC is expressed through the downward control principle for steady-state conditions.The eddy PV flux in K yy can also be related to the mass flux and E-P flux under the quasi-geostrophic assumption (Andrews et al., 1987;Schneider, 2005).Conventional data assimilation may degrade the wave mean flow and wave mixing relationships.To evaluate the short-term dynamical balance in the reanalysis products, we computed the temporal correlation between the E-P flux divergence and the mass stream function, and between the E-P flux divergence and K yy (l = 1 day).The E-P flux was estimated based on the MIM zonal mean momentum equation, in which forcings by sub-grid processes were not considered.Momentum changes through sub-grid processes such as gravity wave drag (GWD) parameterization need to be taken into consideration for strict momentum budget analysis.However, these data were not provided in all the data sets, and hence could not be considered in the analysis.The wave-driven adjustment time for the wave mean flow interactions is typically several weeks in the lower stratosphere (Haynes et al., 1991); this effect was also not considered.Instantaneous analysis fields with no time lag were used for the correlation calculation.We here discuss the relationship between the E-P flux divergence at each vertical level and the MMC at 440 K and between the E-P flux divergence at each vertical level and K yy (l = 1 day) at 560 K.The MMC at 440 K and K yy at 560 K are driven by wave forcings at mostly the same vertical levels (approximately between 500 and 650 K), as discussed below.
As shown in Fig. 10 for ERA-Interim, the mass stream function at 440 K at NH mid-latitudes shows a large correlation with the E-P flux divergence at higher levels, because of downward control.The correlation coefficient averaged be-tween 500 and 650 K at the NH mid-latitudes is larger in the new data sets (approximately 0.5-0.6)than in the old data sets (approximately 0.4) except for JRA (Table 9).This result suggests that the physical balance associated with the wave mean flow interaction is represented more strictly in the new data sets for the NCEP and ERA data sets.JRA-25 and JRA-55 both reveal large correlation coefficients (approximately 0.6).
K yy (l = 1 day) at 560 K shows strong correlations with the E-P flux divergence around 500-650 K in ERA-Interim (Fig. 10).This confirms that the reanalysis products represent local wave forcings that drive eddy mixing.The correlation coefficient averaged between 550 and 650 K in the NH mid latitudes is larger in the new data sets than in the old data sets by 0.04-0.07(i.e., by 24-37 %, Table 9).This result suggests that the wave-mixing relationship is also more accurately represented in the new data sets.Note that not only the dynamic balance in analysis increments, but also characteristics of atmospheric diffusivity in the forecast model (e.g., associated with choice of transport scheme and numerical diffusion) could influence the wave-mixing relationship in the reanalysis products.Also note, that an alternative diagnosis of eddy mixing and its trends using effective diffusivity based on ERA-Interim and JRA-55 has recently become available by Abalos et al. (2016).

Implications for AoA long-term variations
Based on analyses of observational data, several observational studies (e.g., Engel et al., 2009;Stiller et al., 2012;Ray et al., 2014;Hegglin et al., 2014) revealed a positive AoA trend at the NH mid-latitudes from the middle to upper stratosphere, but for different periods (1975-2005in Engel et al., 2009, 2002-2010in Stiller et al., 2012, 1975-2012in Ray et al., 2014, and 1986-2010in Hegglin et al., 2014).Hegglin et al. (2014) also found a reversed trend in the NH lower stratosphere in 1986-2010.Note that the analysis by Hegglin et al. (2014) was based on observed changes in stratospheric water vapour and not the typical age tracers used in the other studies (SF 6 and CO 2 ), and therefore may only be seen as an approximation of age-of-air.
Using ERA-Interim, Monge-Sanz (2012) and Diallo et al. (2013) have demonstrated an AoA decreasing trend in the NH lower stratosphere and an increasing trend in the middle stratosphere in 1990-2009 and 1989-2010, respectively, consistent with the available observational estimates.As discussed by Garny et al. (2014) and Ploeger et al. (2015a, b), it is essential to take the effects of mixing along the transport pathway into consideration, on top of the effects of the MMC for understanding AoA variations.Interpreting AoA variations simply in terms of changes in the local MMC and eddy mixing effects can be misleading, because the AoA of a given air parcel is the result of the integrated effect of the various tendencies along the parcel trajectory.Although the local analysis given in this study does not allow for the air par-   cel history and transport pathway variations to be analyzed, it can provide useful information on the general structure of the transport intensity and its potential impacts on AoA variations.
To investigate whether any reanalysis data have the potential to reveal useful implications of long-term AoA variations, we summarize long-term variations in the three important transport processes: the tropical upward mass flux, the relative contribution of eddy mixing in the subtropics, and the mid-latitudes mean poleward motions (v * ).We consider that, rather than the mixing strength in the subtropics itself, the relative ratio of the mean and eddy transport is essential to understanding AoA variations, since it determines the ratio of fresh air entering and old air recirculating between the tropics and extratropics.Eddy mixing across the MMC overturning may recirculate old air masses from the extratropics to the tropics, whereas mean poleward flows carry fresh air masses from the tropics to the extratropics, and may return old air masses that were recirculated from the extratropics by eddy transport.Here, we focus on transport pro-cesses during winter to explain their possible influences on annual mean AoA variations.In winter, atmospheric waves are active, and strongly induce both MMC and eddy mixing.Abalos et al. (2015) showed that the DJF trends make a major contribution to the overall structure of the annual mean trends of BDC in the NH.However, the influence of transport processes in other seasons cannot be neglected (Konopka et al., 2015).Therefore, the implication obtained from the local analysis of the wintertime transport processes in this study is limited.Explicit calculations of AoA using the reanalysis data sets are required to provide further insights into the role of each transport process on AoA variations.
As shown by Fig. 11, during the 34 years, from 1979 to 2012, the tropical upward mass flux shows a positive trend from the lower to middle stratosphere, except for ERA-Interim in DJF (consistent with the results of Seviour et al., 2011 andAbalos et al., 2015), and NCEP-NCAR and ERA-Interim in JJA (Fig. 11).The relative importance of the eddy transport compared to the mean transport shows an increasing trend in ERA-Interim and JRA-55 from the lower to   middle stratosphere, which may act to increase AoA, especially in the BDC deep branch.The value of v * at the NH mid-latitudes shows a strengthening trend in the lower stratosphere in most data sets, except for NCEP-CFSR below 440 K.The v * trend varies largely with height in ERA-Interim, showing a sharp positive trend peak around 400-430 K, and a negative trend above 520 K. To summarize, AoA derived using the reanalysis products in the NH middle stratosphere for the 34 years can be expected to increase in ERA-Interim, and probably also in JRA-55, because of the large increase in the contribution of the eddy mixing in the subtropics and the weakened mean poleward motion in the deep branch.AoA in the NH lower stratosphere can be expected to decrease in all the data sets, except in NCEP-CFSR.
In the SH, the relative importance of the eddy transport to the mean transport in the subtropics shows an increasing trend only in ERA-Interim in the lower stratosphere, and in ERA-Interim and NCEP-NCAR in the middle stratosphere (above 460 K).The value of v * at the SH mid-latitudes shows a strengthening trend in all the data sets, except for ERA-Interim above 500 K.These changes suggest the possibility of an AoA decreasing trend for both the lower and middle stratosphere in the SH in most data sets except for ERA-Interim during the 34 years.The analysis of ERA-Interim may suggest the possibility of a weaker decreasing or a weak positive AoA trend in the middle stratosphere.
The decadal scale variation in the transport processes will also cause changes in the AoA trend on similar timescales.In the last 12 years (Fig. 12), the upward mass flux trend becomes negative in the lower stratosphere in most data sets, and shows larger increasing trends in the middle stratosphere in all the data sets in DJF.The contribution of the eddy transport in the NH subtropics becomes even more important in the last 12 years in ERA-Interim.The value of v * at the NH mid-latitudes tends to weaken in most data sets, except for ERA-Interim below 430 K, and all the data sets above that level in the last 12 years.Because of the larger negative trend in v * and the increased contribution of the eddy transport in the subtropics, a larger increase in AoA in the NH is expected to be derived using ERA-Interim during the last 12 years than during the 34 years, as suggested by Ploeger et al. (2015b).
In the SH, ERA-Interim reveals a large decreasing contribution of the eddy transport in the last 12 years, which may indicate a decrease in AoA in the SH stratosphere in the last 12 years, as consistently revealed by Stiller et al. (2012) for the 2002-2010 period.JRA-25 reveals substantial decadal scale changes in the transport processes in both hemispheres

Implications for future developments of reanalyses
Although the reanalysis systems have been individually updated at each operation center, similar aspects were encountered regarding the effects of the system updates on the MMC and eddy mixing, as follows.
1.Both the BDC shallow and deep branches reveal weaker MMC mean intensity (by 15-60 % in the NH) and smaller interannual variability (by 10-60 % in the SH) in the new data sets than in the old data sets.This tendency generally will tend to increase AoA throughout the NH stratosphere.
2. The contribution of planetary-scale mixing to eddy mixing was stronger in the new data sets than in the old data sets in the NH lower stratosphere, whereas that of smallscale mixing was smaller.The weaker small-scale mixing in the new data sets may be a result of reduced spurious eddies associated with analysis increments using the flow-dependent analysis (but not in NCEP-CFSR) and the use of the higher forecast model resolution.
3. The relative importance of the eddy transport to the mean transport in the NH is larger in the new data sets by 5-27 % in the extratropics, and by 5-14 % in the subtropics.The larger eddy contribution in the subtropical lower stratosphere may cause an older AoA in the entire stratosphere in the new data sets.
4. The wave mean flow relationship and the wave-mixing relationship in the NH surf zone are more accurately represented in the new data sets for the NCEP and ERA data sets.
Updated systems can be expected to provide better representations of both the MMC and eddy mixing, because of reduced systematic errors of the forecast model, the application of improved bias correction algorithms, and advanced data assimilation techniques.Nevertheless, it is not straightforward to identify which changes are mainly responsible for the differences in the MMC and eddy mixing.Monge-Sanz et al. (2007, 2012) investigated that ERA-Interim benefits using an omega-equation balance operator in the background constraint, which is likely to have reduced spurious propagation of eddy motion associated with analysis increments.They also suggested other important factors, such as reduced stratospheric temperature biases of the forecast model and improved bias corrections for satellite radiance data.In ERA-40 and JRA-25, systematic analysis increments were introduced to compensate for model temperature biases, which are thought to cause an overly strong BDC (Uppala et al., 2005).Kobayashi and Iwasaki et al. (2015) found that in JRA-25, the temperature analysis increment evolves from the 1990s to the 2000s, associated with changes in forecast model biases and lack of effective bias corrections, either for the forecast model and assimilated measurements.As a result, the MMC structure in JRA-25 was unrealistically distorted in a particular period corresponding to the changes in external forcing.
The decadal scale variations of the transport processes in the reanalysis could be introduced artificially because of measurement discontinuity.Stratospheric temperatures in ERA-Interim are known to be affected by the introduction of AMSU-A data in 1998 and radio occultation data at the end of 2006 (Dee et al., 2011).Improved bias correction schemes have been implemented in recent reanalysis systems (e.g., variational bias correction in JRA-55, Kobayashi et al., 2015, andERA-Interim, Dee andUppala, 2009).However, it is still difficult to distinguish between forecast models and measurements biases, and to correct them properly.
In addition, the reanalysis quality could be strongly affected by the performance of the forecast model.Any changes in wave propagation and breaking lead to changes in both the MMC and mixing, but current models have large uncertainties in representing various waves, including forcings from GWD parameterization (e.g., Sigmond and Shepherd, 2014).Further efforts on model development are still important.

Conclusions
We compared the characteristics of the stratospheric MMC and eddy mixing in the six reanalysis products for the period 1979-2012 based on mass-weighted isentropic zonal means.Both of the mean and eddy transport processes play important roles in determining distributions and variations of chemical tracers and AoA, thus it is important to clarify whether there are any reanalysis data that can be used to investigate long-term BDC variations.
The mean K yy value differs largely among the data products both in the extratropics and subtropics.The contribution of planetary-scale mixing in the NH lower stratosphere is larger in the new data sets, whereas that of small-scale mixing is smaller.The weaker small-scale mixing in the new data sets is considered to be a result of reduced spurious eddies associated with analysis increments in the 4D-VAR analysis.Isentropic mixing shows a strengthening trend in the NH middle stratosphere surf zone in all data sets, associated with strengthened large-scale mixing.

Atmos
There are large differences in the MMC strength among the data sets, but are were some similar characteristics.In the NH, the MMC is stronger in the new data sets (NCEP-CFSR, ERA-Interim, and JRA-55) than in the old data sets (NCEP-NCAR, .All the data sets show a strengthening trend in the BDC shallow branch in the winter hemispheres, consistent with model simulations (e.g., Butchart et al., 2010).The MMC in the BDC deep branch showed a weakening trend only in ERA-Interim at the NH mid-latitudes and in ERA-Interim and JRA-25 in the SH mid-latitudes in winter.The vertical structure in the changes as seen in ERA-Interim that would lead to a decrease in AoA in the lower stratosphere and an increase in AoA in the middle (and upper) stratosphere are broadly consistent with inferences made from changes observed in stratospheric water vapour (Hegglin et al., 2014).
We also evaluated the relative importance of the mean and eddy transports.The contribution of eddy mixing was generally stronger in the new data sets than in the old data sets in the NH.In the subtropical lower stratosphere, the relative contribution of eddy transport showed an increasing trend in data sets, except for NCEP-NCAR and NCEP-CFSR in the NH, and only in ERA-Interim in the SH.These trends reveal changes in the relative effectiveness of eddy-induced recirculation from the extratropics to the tropics, and the mean poleward transport of fresh air into the extratropics; they are considered to be important in understanding AoA variations.
The transport analysis suggests that ERA-Interim provides a consistent transport picture, with AoA trends derived based on observations in both the deep and shallow branches, whereas other data sets would provide different implications.The increasing trend in AoA in the NH middle stratosphere can be understood as a result of the weakening MMC in the deep branch (found only in ERA-Interim), together with the increased contribution of eddy transport compared with the mean transport in the subtropical lower and middle stratosphere (found in ERA-Interim and JRA-55).The decreasing trend in AoA in the SH lower stratosphere could be a result of the strengthening MMC trends in the shallow branch (large in ERA-Interim and JRA-25).All the reanalysis data sets revealed decadal scale variations in the strength of both MMC and eddy mixing, which may result in significant AoA trend variations.For instance, the increased contribution of the eddy transport in the subtropics and the weakened mean poleward motion in the middle stratosphere during the last 12 years  compared to the first 22 years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) in the NH suggests larger increasing trends in AoA in ERA-Interim and NCEP-NCAR in the last 12 years.
Differences in data assimilation schemes and forecast models are thought to cause differences in the expression of the MMC and eddy mixing in the BDC among reanalysis products.Our analysis suggests that advanced reanalysis products are potentially useful for studying long-term BDC variations, because of the flow-dependent background error covariance, assimilation of the observations at the exact time, balance operators, and improved bias correction algo- rithms used in these reanalyses.However, there still seem to be problems that cause unrealistic atmospheric variations (e.g., large differences in long-term variations between ERA-Interim and JRA-55) associated with discontinuities in the assimilated measurements and large uncertainties in the forecast models.Further efforts are essential to improve the reanalysis systems so that the reanalysis products become more consistent among each other and that they can be used to study the details in BDC variations.

Figure 7 .
Figure 7. Temporal variations of the annual mean total mean upward mass flux (in 10 10 kg s −1 ) at (top) 440 K and (bottom) 560 K.

Figure 8 .
Figure 8. Same as Fig. 1 but for the relative importance of the eddy transport to the total meridional transport.The value larger (smaller) than 0.5 indicates that the eddy transport (mean transport) is dominant in the meridional transport.

Figure 9 .
Figure 9. Same as Fig. 2 but for the ratio of the relative importance of the eddy transport to the total meridional transport.

Figure 10 .
Figure 10.Latitude and potential temperature cross section of the correlation (left) between the 6-hourly mass stream function at 50 • N, 440 K and the 6-hourly E-P flux divergence at each grid point, and (right) between Kyy (l = 1 day) at 50 • N, 560 K, and the 6-hourly E-P flux divergence at each grid point in DJF in 2000-2012.Correlation coefficients at the 95 % confidence level (following t test) are indicated by hatching.

Figure 11 .
Figure11.Vertical profiles of the linear trend of the tropical upward mass flux (left, in 10 10 kg s −1 decade −1 ), relative contribution of the eddy transport to the total meridional transport averaged over 15-25 • (second left, in decade −1 ), and the mean-meridional velocity v * averaged over 30-50 • (m s −1 decade −1 ) during DJF in the NH (upper panels) and during JJA in the SH (lower panels) for the34 years,  1979-2012.

Table 1 .
Mean value (Mean) over 34 yearsof K yy (l = 1 month) (in 10 6 m 2 s −1 ) averaged over 40-60 • N at 440 and 560 K and over 15-25 • N at 440 K in DJF, and averaged over 40-60 • S at 440 K and 560 K and over 15-25 • S at 440 K in JJA.The standard deviation (SD) is shown in brackets.The linear trend slope (Trend, in % decade −1 ) is also shown.The linear trend with a confidence level greater than 95 % is shown in bold.

Table 2 .
Linear trend slope of K yy

Table 3 .
Same as in Table1, but for the mass stream function (in 10 10 kg s −1 ) averaged over 40-60 • N in DJF, and over 40-60 • S in JJA at 440 and 560 K.The linear trend with a confidence level greater than 95 % is shown in bold.

Table 4 .
Same as in Table1but for the annual mean total upward flux (in 10 10 kg s −1 ) and its linear trend slope (in % decade −1 ) at 440 and 560 K.The linear trend with a confidence level greater than 95 % is shown in bold.

Table 5 .
Linear trend of the annual mean tropical total mean upward flux (in % decade −1 ) during the first 22 years, and during the last 12 years (2001-2012) at 440 and 560 K.The linear trend with a confidence level greater than 95 % is shown in bold.

Table 6 .
Relative contribution of different wave components to the zonal eddy PV component (in %) and the linear trend slope (in % decade −1 in brackets) averaged over 40-60 • N at 440 and 560 K in DJF in 1979-2012.None of the linear trend has a confidence level of 95 %.
Interim corresponds to the increasing trend in K yy that only appeared in ERA-Interim and NCEP-NCAR.

Table 9 .
Correlation coefficient between the 6-hourly mass stream function (Mean) at 50 • N, 440 K, and the 6-hourly E-P flux divergence at each grid point averaged between 500 and 650 hPa and 50-60 • N, and between K yy (l = 1 day) at 50 • N, 560 K, and the 6-hourly E-P flux divergence at each grid point averaged between 550 and 650 hPa and 45-55• N in DJF in 2001-2012 (1990-2001 for ERA-40).