Future Changes in the Quasi-Biennial Oscillation Under a Greenhouse Gas Increase and Ozone Recovery in Transient Simulations by a Chemistry-Climate Model

In the equatorial stratosphere, the zonal wind direction changes from westward to eastward and vice versa at intervals of from less than two to about three years and centered at about 28 months (Reed et al., 1961; Veryard & Ebdon, 1961). This phenomenon is called the quasibiennial oscillation (QBO) from the length of its period (Angell & Korshover, 1964). The QBO initiates in the upper stratosphere, then descends gradually at a rate of about 1 km per month with an amplitude of about 20 ms-1, and diminishes near the tropopause (~16 km). The QBO structure for zonal wind and temperature is highly symmetric in longitude and latitude with a half-width of about 12° latitude, and the QBO dominates the annual and semiannual cycles in the equatorial stratosphere (e.g., Baldwin et al., 2001; Pascoe et al., 2005). The QBO can also be seen in changes in the abundances of ozone (e.g., Angell & Korshover, 1964; Randel & Wu, 1996) and other trace gases (e.g., Luo et al., 1997; Dunkerton, 2001) because their transport and chemistry are strongly dependent on the wind and temperature fields, respectively.


Introduction
In the equatorial stratosphere, the zonal wind direction changes from westward to eastward and vice versa at intervals of from less than two to about three years and centered at about 28 months (Reed et al., 1961;Veryard & Ebdon, 1961). This phenomenon is called the quasibiennial oscillation (QBO) from the length of its period (Angell & Korshover, 1964). The QBO initiates in the upper stratosphere, then descends gradually at a rate of about 1 km per month with an amplitude of about 20 ms -1 , and diminishes near the tropopause (~16 km). The QBO structure for zonal wind and temperature is highly symmetric in longitude and latitude with a half-width of about 12° latitude, and the QBO dominates the annual and semiannual cycles in the equatorial stratosphere (e.g., Baldwin et al., 2001;Pascoe et al., 2005). The QBO can also be seen in changes in the abundances of ozone (e.g., Angell & Korshover, 1964;Randel & Wu, 1996) and other trace gases (e.g., Luo et al., 1997;Dunkerton, 2001) because their transport and chemistry are strongly dependent on the wind and temperature fields, respectively.
The QBO plays a very important role in the stratosphere not only because the QBO is the largest variation in the equatorial stratosphere but also because its effects significantly extend to the northern extratropical stratosphere (e.g. Holton & Tan, 1980, with the result that the stratospheric polar vortex tends to be weaker (stronger) during the QBO easterly (westerly) phase in the lower stratosphere. In addition, the QBO effect can also be seen in the extratropical troposphere at or near the surface not only in the Northern Hemisphere (NH) (e.g., Coughlin & Tung, 2001;Marshall & Scaife, 2009) but also in the Southern Hemisphere (SH) (Kuroda & Yamazaki, 2010), indicating that it has dynamical impacts on the tropospheric circulation even though the physical mechanism is still unknown.
The QBO effect on the polar vortex, which has been confirmed by many statistical analyses, is referred to as the Holton-Tan relation (e.g., Labitzke, 1982;Naito & Hirota, 1997;Lu et al., 2008). The Holton-Tan relation is ascribed to modulation by the QBO of the strength of the effective waveguide for the mid-latitude planetary waves propagating through the winter 356 stratosphere (e.g., Holton & Tan, 1980. Although the Holton-Tan relation has been reproduced by general circulation models (GCM) (e.g., Niwano & Takahashi, 1998;Hamilton 1998, Anstey et al., 2010 and chemistry-climate models (CCM) (Naoe & Shibata, 2010), the details of its mechanism are not yet well understood because the Eliassen-Palm (EP) flux diagnostics do not show a more poleward propagation in the mid-latitude stratosphere during the easterly QBO phase. Rather, planetary waves propagate more equatorward as well as more upward during the easterly QBO phase in both 25-year (1980-2004) data of the reanalysis and a CCM ensemble simulation (Naoe & Shibata, 2010) under the REF-1 scenario of CCM validation phase 1 (CCMVal-1) activity for stratospheric processes and their role in climate (SPARC) (Eyring et al., 2005). Furthermore, the polar vortex intensity depends not only on the QBO but also on the 11-year solar cycle (e.g., Labitzke & van Loon, 1999); as a result, the Holton-Tan relation varies irregularly or nonstationarily with the solar activity (e.g., Gray et al., 2001;Labitzke, et al., 2006;Lu et al., 2009), making the circumstances of the Holton-Tan relation more complicated than what Holton & Tan (1980 first suggested. The mechanism of the QBO is generally explained with a wave-mean flow interaction, wherein gravity waves and/or equatorial waves excited in the troposphere deposit their westward or eastward momentum, depending on the mean flow direction during their vertical propagation in the stratosphere, as has been demonstrated with numerical models (e.g., Lindzen & Holton 1968;Plumb 1977). This mechanism, which has also been reproduced in a laboratory experiment of Plumb & McEwan (1978), is one reason why the oscillation does not have a narrow annual or semiannual period but a peculiar, broad 28month period. The QBO is hence strongly dependent on the wave source and the propagation conditions, that is, the wave generation in the troposphere and the background wind field in the troposphere and stratosphere.
Under the changing climate resulting from increases in greenhouse gases (GHG), the stratosphere is projected to continue to cool because of strong CO 2 longwave emission, in a sharp contrast to the tropospheric warming caused by the sea surface temperature (SST) rise driven by the enhancement of downward longwave radiation by CO 2 and H 2 O. This raises the question as to how the QBO might behave in the future under global warming caused by increased GHGs. To date, two equilibrium GCM simulations under doubled CO 2 levels have projected the future QBO. Giorgetta & Doege (2005) investigated the effect of variations in the source strength of parameterized gravity wave forcing (GWF) with a middle atmosphere GCM (Roeckner et al., 2003) at the T42L90 resolution and showed that the QBO period is shortened when the source strength of GWF increases. The future QBO periods are 26, 22, and 17 months for a 0%, 10%, and 20% increase in source strength, respectively, against a 29month period under the current CO 2 conditions. Kawatani et al. (2011) carried out a similar experiment with a GCM (Hasumi & Emori, 2004) at the T106L72 resolution. They did not employ a parameterized GWF, so that the wave forcing of the QBO was by resolved waves alone. In their simulation, the future QBO has a period that is 1-5 months longer than that of the control QBO, which has a period of 24 months. Moreover, the future QBO amplitude is smaller than the control; in particular, the amplitude is much smaller below 50 hPa. Thus, the two GCM results show an opposite change in the QBO period, although both Giorgetta & Doege (2005) and Kawatani et al. (2011) attributed the future QBO change mainly to changes in the residual circulation (mean tropical upwelling) and to the enhanced wave forcings. Aside from the use of a parameterized or resolved GWF, the difference in the QBO In the polar regions, ODS increases led to severe ozone depletion in spring through the heterogeneous reactions on the polar stratospheric clouds (PSCs), which manifested as the Antarctic ozone hole (Chubachi, 1984;Farman et al., 1985), defined as an area with less than 220 Dobson Units (DU) of total ozone. The ozone hole developed rapidly from about 1980 to the mid-1990s; since then, the area with severely depleted ozone has remained nearly constant, except during a few years, when wave activity was unusually active (e.g., WMO, 2011). Severe ozone depletion also occurs in the northern polar region if the winter stratospheric temperature is cold enough for PSCs to form (e.g. Rex et al., 2006), though not as frequently as in the southern polar region.
Here, to address future QBO behavior under climate changes due to increasing GHGs and decreasing ODSs, Meteorological Research Institute (MRI) CCM transient simulation data from 1960 to 2100 are analyzed focusing on the QBO trend and its forcings. The future QBO effect on the polar vortex, that is, the future Holton-Tan relation, in MRI-CCM simulations is described in another paper (Naoe & Shibata, 2012). Three simulations, using the three scenarios of CCMVal phase 2 (CCMVal-2, Eyring et al., 2008;SPARC, 2010) for the future projection of the stratospheric ozone and related species, are used. The scenarios used in this study specify only anthropogenic forcings of GHGs/SST and/or ODSs, and include no natural forcings, that is, neither the 11-year solar cycle nor volcanic aerosols. The simulation data are thereby free from the complicated overlapping effects among the QBO, solar cycle, and volcanic aerosols (e.g., Lee & Smith, 2003;Yamashita et al., 2011), which would otherwise reduce the confidence level of each detected signal.
The rest of this paper is organized as follows. Section 2 describes the model and simulation conditions, and section 3 presents past and future climatologies of zonal wind and temperature. Section 4 describes the QBO and its forcings in the past, and section 5 treats the evolutions of the QBO and of its forcings. The results are discussed in section 6, and section 7 is a conclusion.

MRI-CCM
MRI-CCM is an upgraded version of the MRI chemistry transport model (Shibata et al., 2005), and its details are described by Shibata & Deushi (2008a, b) and references therein. The dynamics module of MRI-CCM is a spectral global model with triangular truncation, a maximum total wavenumber 42 (about 2.8° by 2.8° in longitude and latitude), and 68 layers in the eta-coordinate extending from the surface to 0.01 hPa (about 80 km). The vertical spacing is about 500 m in the stratosphere between 100 and 10 hPa, with tapering below and above. Non-orographic gravity wave (GW) forcing by Hines (1997) is incorporated, and GW source strength is enhanced symmetrically with respect to latitudes between 30°S and 30°N by superimposing a Gaussian function (15 e-folding latitude) source (0.7 ms -1 ) on an isotropic source (2.3 ms -1 ). The gravity wave spectrum is isotropically launched toward eight equally spaced azimuths (north, northwest, west etc) at the lowest level with horizontal wavenumber k * = 5.010 -6 (m -1 ).
Biharmonic (∆ 2 ) horizontal diffusion is weakened only in the middle atmosphere from that in the troposphere in order to spontaneously reproduce the QBO in zonal wind while minimizing the changes in the troposphere (Shibata & Deushi, 2005a, b). The e-folding time at the maximum total wavenumber (n = 42) is 18 h below 100 hPa, 100 h above 150 hPa, and intermediate values in between. This e-folding time of 100 h in the middle atmosphere results in an e-folding time of about 3.6 years at n = 10; almost all of the power (>99%) is included within this total wavenumber for a representative equatorial wave of Gaussian meridional structure centered at the equator with a 15° e-folding latitude (Shibata & Deushi, 2005a). In addition, vertical diffusion is not applied in the middle atmosphere in order to keep the sharp vertical shear in the QBO.
The chemistry-transport module employs a hybrid semi-Lagrangian transport scheme compatible with the continuity equation to satisfy the mass conservation. The horizontal form is an ordinary semi-Lagrangian scheme, and the vertical form is equivalent to a massconserving flux form in the transformed pressure coordinate specified by the vertical velocity. The vertical procedure employs the piecewise rational method (Xiao & Peng, 2004), and the horizontal procedure uses quintic interpolation. The chemistry scheme treats 36 long-lived species including 7 families, and 15 short-lived species with 80 gas-phase reactions, 35 photochemical reactions and 9 heterogeneous reactions on PSCs and sulfate aerosols. Since MRI-CCM does not include all of the chlorine species, the chlorine species (CCl 4 , CFCl 3 , CF 2 Cl 2 , and CH 3 Cl) are lumped, except for BrCl, such that the modeled total chlorine corresponds to the input total chlorine at the surface.

Simulations with the CCMVal-2 scenarios
Two scenarios in CCMVal-2 are used for reference simulations (Eyring et al., 2008): REF-B1 is used for simulations of the past, from 1960 to 2006, with all observed forcings of GHGs, ODSs, SST, solar spectral irradiance variability, and (volcanic) aerosols in the stratosphere, and REF-B2 is used for both past and future simulations, from 1960 to 2100, with time-evolving forcings of GHGs, ODSs, and SST and with fixed solar minimum and background aerosol conditions. Fixed conditions are used because plausible specification of future solar variability and volcanic aerosols is very difficult. SST/sea-ice data for the entire period from the spin-up (covering the 1950s) to 2100 are taken from a transient simulation (1%/year CO 2 increase) by the MRI coupled ocean-atmosphere model MRI-CGCM2.3.2 (Yukimoto et al., 2006), in which the atmospheric GCM is nearly the same as in MRI-CCM. Note that even the past simulation from the 1950s uses not observed but simulated SST/sea-ice data in order to avoid gap in SST/sea-ice data, stemming from model SST biases, in the vicinity of the present condition.
GHG and ODS abundances are assumed to be geographically uniform and specified at the surface (Eyring et al., 2008). The future scenarios for GHGs and ODSs are taken from the Special Report on Emissions Scenarios (SRES) A1B GHG scenario (IPCC, 2000) and the adjusted A1 halogen scenario. The A1B GHG scenario represents a future world of very rapid economic growth, a global population that peaks in the mid-twenty-first century and declines thereafter, and the rapid introduction of new and more efficient technologies, balanced across all GHG sources. On the other hand, the adjusted A1 halogen scenario includes the earlier phase-out of hydrochlorofluorocarbons (HCFCs) that was agreed to by the Parties to the Montreal Protocol in 2007, with other species such as CFCs, halons, and other non-HCFCs remaining identical to the original A1 scenario (WMO, 2007). In the original A1 scenario, future production of ODSs is determined from the Montreal Protocol limitations or the most recent annual estimates, whichever is less; future emissions are determined in order to yield banks consistent with those of IPCC/TEAP (2005); and the future annual fractional bank release is adjusted so as to attain the IPCC/TEAP (2005) bank for 2015.
Beside the two reference simulations, two other simulations were performed by fixing either ODSs or GHGs at those values at 1960 levels; these simulations correspond to two of the eight sensitivity simulations of CCMVal-2 (Eyring et al., 2008). The first simulation (SCN-B2b, i.e., fixed halogens) uses fixed ODSs and the same time-evolving GHGs and SST/seaice data as REF-B2, whereas the second (SCN-B2c, i.e., no-climate-change) uses fixed GHGs and annually repeating SST/sea-ice data and the same time-evolving ODSs as REF-B2. Note that the SST/sea ice at 1960 levels is the 10-year average centered on 1960.   CH 4 , and N 2 O) in the A1B scenario and ODSs (as total chlorine, CCl y , and total bromine, CBr y ) in the A1 adjusted scenario for 1960 to 2100. CO 2 and N 2 O abundances increase monotonically throughout the entire period, reaching maxima of about 700 ppmv and 370 ppbv, respectively, in 2100, whereas CH 4 abundances attains a maximum of about 2.4 ppmv around the 2050s and then decreases to about 2.0 ppmv in 2100. We henceforth refer to this GHGs evolution simply as the GHG increase, neglecting the non-monotonic behavior of CH 4 . CCl y increases rapidly from 0.8 ppbv around 1960 to a maximum of about 3.7 ppbv in the first half of the 1990s and then gradually decreases to about 1.3 ppmv in 2100. CBr y also behaves qualitatively like CCl y but it starts to increase much later (from about 1980) and attains a maximum of about 17 pptv around 2000.

Data processing
As observed data, for comparison with past simulations, we used reanalysis data sets compiled by the European Centre for Medium-Range Weather Forecasts (ECMWF): namely, ERA40 (Uppala et al., 2005) from 1958 to 1988, and ERA-Interim (Dee et al., 2011) from 1989 to 2009. The two reanalysis data sets are smoothly merged such that the data for 1989 overlap, without any gap in any month, thus producing a long, continuous time series of data from 1958 to 2009, which is henceforth referred to simply as ERA-Interim. Monthlymean data are used in this study for both the observation and the simulations.
A Lanczos filter (e.g., Duchon, 1979) was applied to the monthly-mean time series to get band-passed data within the QBO period (20-40 months) or low-passed data in lowfrequency ranges, with periods much longer than the QBO period. The cut-off for the low frequency data period was set at 120 months. The running-mean was used as an alternative way of deriving low-frequency data. Student's t-test was used to evaluate the linear trend in the QBO and related quantities over the entire integration period, and results were considered statistically significant at the 95% or 99% confidence level.

Zonal mean climatology
Simulated annual-and zonal-mean temperature and zonal wind in the past 30 years  are displayed for B21, Sb, and Sc in Figs. 2 and 3, together with observed values. Although the overall features of temperature are very similar between the observation and the simulations, the simulated temperatures show common biases in the mid-and highlatitude stratosphere in both the NH and SH, which are substantially due to biases during the winter. In the NH, the simulated temperatures show local minima at around 200 hPa in the polar lower stratosphere, above which temperatures increase monotonically upward to the stratopause. On the other hand, observed temperatures exhibit a local warm maximum at around 200 hPa and a local cold minimum at around 30 hPa, in the polar lower stratosphere, indicating that MRI-CCM simulations show a cold bias in the lower stratosphere and a warm bias in the upper stratosphere at high latitudes. Fig. 2. Annual-and zonal-mean temperature. Latitude-pressure cross sections for annualand zonal-mean temperature (K) in the past 30-year (1970-1999) are for (a) ERA-Interim, (b) B21, (c) Sb, and (d) Sc. The contour interval is 10 K, and blue colors represent lower temperatures and red colors higher temperatures. Fig. 3. Annual-and zonal-mean zonal wind. Same as Fig.2 except for zonal wind (ms -1 ). The contour interval is 5 ms -1 , and blue colors represent positive (westerly) winds and red colors negative (easterly) winds.
Consistent with these temperature biases, the simulated zonal winds also show common biases in the subtropical jet and the polar night jet as a result of the thermal wind balance. The core velocity of the subtropical jet is too strong, and its separation from the polar night jet is too weak; that is, the zonal wind at around 50 hPa and 60°N is too strong, and the polar night jet in the upper stratosphere is too weak. In the SH, there are also cold biases in the mid-and high-latitude lower stratosphere at around 200 hPa. Similar to the NH, in the SH the simulated subtropical jet is stronger than the observed one, but the simulated polar night jet is weaker.

Future climate change in temperature and zonal wind
Figures 4 and 5 display the changes in the future 30-year period (2070-2999, 100 years later than the past 30-year period ) in annual-and zonal-mean temperature and zonal wind. During this future period, the ODS abundances are almost restored to pre-1970s levels ( Fig. 1), so that the forcing is almost entirely due to GHGs alone. The temperature changes in the three B2 members (B21, B22, and B23) and in the fixed-halogen simulation Sb exhibit very similar CO 2 -(GHG-) induced patterns (e.g., Yukimoto et al., 2006;IPCC, 2007): warming in the troposphere through SST warming (indirect CO 2 effect) and cooling in the stratosphere (direct CO 2 effect) (Kodama et al., 2007). In addition, weak but substantial cooling is seen in the tropical lower stratosphere. The tropospheric warming reaches a maximum of about 4 K in the tropical upper troposphere, reflecting the enhanced latent heat release resulting from intensified convective activity caused by the SST rise. In the B2 simulations, the maximum stratospheric cooling, to about -6 K, in the upper stratosphere is caused both by enhanced CO 2 terrestrial radiation and by recovered ozone solar radiation as described below.
(c) (a) (b) Fig. 4. Future changes in temperature. Latitude-pressure cross sections of future changes (difference between 2070-2099 and 1970-1999) in annual-and zonal-mean temperature (K) for (a) B21, (b) Sb, and (b) Sc. The contour interval is 1 K for B21 and Sb and 0.5 K for Sc, and blue colors represent cooling and red colors warming.
In the no-climate-change simulation Sc, on the other hand, strong warming occurs in the upper stratosphere above about 5 hPa in both hemispheres, except in the southern high latitudes, where prominent warming occurs in the lower stratosphere at 300-30 hPa, centered at 100 hPa, and the cooling above 30 hPa. This temperature change pattern is completely different from the CO 2 -induced pattern and is due solely to the ODS decrease, or equivalently, the ozone recovery, because the average ODS concentration was higher in the past 30 years  than in the future 30 years (2070-2099) (Fig. 1). In fact, stratospheric ozone abundances were projected to recover to or exceed pre-1980 levels globally in the future, except for the tropical lower stratosphere due to the increased inflow of ozone-poor tropospheric air mass (e.g., SPARC, 210). In other words, the upper stratospheric cooling in the past 1970-1999 was due to an ozone decrease caused by gas-phase reactions with reactive chlorine and bromine atoms photodissociated from ODSs under strong ultraviolet solar radiation, whereas the pattern in the southern high latitudes was associated with severe ozone depletion (i.e., the ozone hole) caused by heterogeneous reactions on PSCs (e.g., SPARC, 2010). This ODS-decrease effect is also included in the B2 simulations; thus, temperature changes similar to those in Sc can be obtained by subtracting the changes in Sb from those in the B2 members (not shown). Zonal wind changes in the three B2 members and Sb display a very similar pattern in the polar night jet regions in both hemispheres and in the equatorial lower stratosphere. Westerly wind intensification extends from the upper region of the subtropical jet upward along the equatorward flank of the polar night jet axis, reaching a maximum at around 70 hPa and 40° in both hemispheres. In the SH, the westerly wind intensification also extends down to the surface at around 55°S, resulting in poleward spreading of the subtropical jet.
In the tropics, westerly wind strengthening (eastward acceleration) occurs above 30 hPa and weakening (westward acceleration) at 50 hPa, below which very weak westerly wind strengthening occurs down to the surface. In the no-climate-change simulation Sc, on the other hand, the zonal wind change pattern is completely different. In the SH, weakening of the polar night jet and both weakening of the poleward and strengthening of the equatorward flank of the subtropical jet down to the surface occur, changes which are almost opposite to those in B2 and Sb in the troposphere above southern high latitudes. This change in the subtropical jet in Sc means that the equatorward shift of the jet is associated with ozone recovery (e.g., Son et al., 2008). In the NH, the poleward flank of the polar night jet is slightly intensified above the middle stratosphere. In the tropics, there is almost no change in the troposphere and stratosphere, consistent with the temperature results.

Evolution of the Brewer-Dobson circulation and the ozone hole area
The GHG increase affects the mass exchange between the troposphere and the stratosphere, which can be evaluated as the upward mass flux across the tropopause. Most GCMs and CCMs (e.g., Butchart et al., 2006;Austin & Li, 2006;Garcia & Randel, 2008) show an acceleration of BDC resulting from an increase in GHGs, but balloon measurements do not necessarily support these model results, even given the large uncertainties and sparse observations in both time and space (Engel et al., 2009). Figure 6 depicts the evolution of the global upward mass flux anomaly (49-month running mean) at 100 hPa, calculated from the residual mass stream function. It is evident that the simulations with the GHG increase scenarios (B21, B22, B23, and Sb) result in a nearly linear increase of the upward mass flux across 100 hPa, and the trend slopes are very similar at about 8×10 9 kgs -1 per century. On the other hand, the simulation results for the no-climate-change scenario Sc exhibit no trend at all. This feature can be also analyzed in the equatorial stratosphere, as displayed in Fig. 7, which depicts the low-frequency (period longer than 10 years) components of the residual vertical velocity * w (in units of 10 -4 ms -1 , averaged from 10°S to 10°N) from 70 to 1 hPa. Note that some high-frequency components remain, due to leakage of the Lanczos filter. The magnitude of * w is very small in the lower stratosphere, with a minimum at 50 hPa; above 50 hPa it increases with altitude, gradually up to about 5 hPa and steeply above that. This vertical structure does not change in the GHG-increase simulations or in the no-climatechange simulation during the whole integration period, consistent with equilibrium CO 2doubling simulations (e.g., Kawatani, et al., 2011).  ODSs will continue to induce an ozone hole over the Antarctic from late winter to spring in the future, as long as their abundances exceed certain values (e.g., SPARC, 2010). Figure 8 exhibits the evolutions of the annual maximum ozone hole area in the ODS time-evolving simulations (B21, B22, B23, and Sc). Different from these simulations, no ozone hole appears in the fixed-ODS simulation Sb, even if the GHGs-induced climate change progresses. Even though the inter-simulation variation is not especially small, the features of all of these simulations are very similar with regard to the long-term evolution of the ozone hole area: (1) initial appearance around 1980, (2) followed immediately by a rapid increase until the later 1990s, (3) then nearly saturated values until about 2010, and (4) followed by a gradual decrease, with disappearance around 2070. Features 1-3 qualitatively correspond to the observed features (e.g., WMO, 2011), aside from the short-term variations. The fact that the evolution of the ozone hole area in Sc is similar to that in the three B2 members indicates that the effect of the GHG increase on the ozone hole in the SH is limited.  Figure 9 shows the observed and simulated zonal-mean zonal wind (averaged from 10°S to 10°N) between 100 and 5 hPa for the 20 years from 1990 to 2009. Monthly averages are subtracted from the values of each month to extract anomalies from the seasonal cycle, so that the semi-annual oscillation (SAO) in the upper stratosphere becomes smaller than or similar to the QBO. It is evident that all of the simulations reproduced the overall features of the observed QBO, such as the stronger westerly shear than easterly shear during the phase transition as well as the stalling of the descent of the easterly winds. Empirical orthogonal function (EOF) analysis revealed that EOF 1 and 2 explain major parts of the variance as in observations (e.g., Wallace et al.,1993) (not shown). The simulated QBOs realistically descend to as low as 100 hPa, whereas their amplitudes are smaller than the observed amplitude, particularly below about 50 hPa in the lower stratosphere, a common bias among three-dimensional models (Takahashi, 1999;Hamilton et al., 1999;Scaife et al., 2000;Giorgetta et al., 2002;Butchart et al., 2003;Shibata & Deushi, 2005a;Watanabe et al., 2008;Kawatani et al., 2011). The differences among the simulations are not necessarily discernible in Fig. 9.  QBO has less power than the simulated QBO in the upper stratosphere. On the basis of the power spectrum distribution, we defined the QBO as the components with a period between 20 and 40 months. In the higher frequency domain (periods shorter than 16 months), the power spectrum is a narrow monomodal distribution concentrated at an annual frequency below 10 hPa, and a bimodal distribution sharply concentrated at annual and semiannual frequencies above 10 hPa, in both the observation and the simulations (not shown). The QBO power is the largest of these two or three peak powers below the middle stratosphere.  Figure 11 displays the composite relations within the QBO cycle at 30 and 70 hPa for B21 between the QBO zonal wind and the three major forcings, parameterized gravity-wave forcing (GWF), resolved-wave forcing (Eliassen-Palm flux divergence, EPD), and vertical advection of zonal momentum ( * (/) w uz   , WUZ), calculated for the entire period of 1960-2099. Note that all of the forcing terms and the zonal wind are band-passed (20-40 months) values with the QBO components extracted by a Lanczos filter. In making the composite figures, we assumed that each cycle expands in 2 phase space with a peak power period of 28 months, even though each cycle had a different period. Although the fine structures such as steep gradient during the phase transitions disappear due to the band-pass filtering, basic configuration among the QBO and its forcings are represented in Fig. 11.

Forcings of the QBO
AT 30 hPa, GWF is retarded by about one-seventh cycle relative to zonal wind, so that positive (eastward acceleration) GWF occurs approximately from after the initial stage of the westerly shear descent (equivalently, u/t > 0) to the initial stage of the easterly shear descent (u/t < 0); conversely, negative (westward acceleration) GWF continues from after Fig. 11. Composite relations among the QBO and its forcings. Composite relations among zonal wind (U, red solid curve with closed circles), parameterized gravity wave forcing (GWF, black dashed curve), resolved wave forcing (EPD, green solid curve), and vertical advection (WUZ, cyan solid curve) at (upper) 30 hPa and (lower) 70 hPa over two cycles for an average QBO cycle (28 months). All quantities are averaged from 10°S to 10°N and bandpass filtered (20-40 months). Units are ms -1 for U, 10 -2 ms -1 day -1 for forcing/advection. the initial stage of the easterly shear descent to the initial stage of the westerly shear descent. The amplitude of WUZ is three-quarters that of GWF amplitude and is almost of opposite phase, thus partially cancelling GWF. Therefore, GWF promotes the descent of the vertical shear, and WUZ hinders it. EPD lags GWF by about one-seventh cycle, and thus lags the zonal wind by approximately one-quarter cycle. Of these three forcings, GWF and WUZ are the largest and second largest terms, but as a result of the partial cancellation due to their being of almost opposite phase, the sum of GWF and WUZ has a similar amplitude to that of EPD. At other altitudes, a similar relation quantitatively holds true, particularly at 20 hPa and above, although the deviation from the relation at 30 hPa increases in the lower stratosphere (50-70 hPa). AT 50 hPa, GWF is nearly in-phase with zonal wind, whereas WUZ and EPD lag the zonal wind by about four-and three-sevenths cycle, respectively, and are thus out of phase with zonal wind. In addition, the amplitudes of GWF and WUZ decrease from around 0.13 and 0.10 ms -1 day -1 , respectively, at 30 hPa to 0.04 and 0.03 ms -1 day -1 , respectively, at 70 hPa, resulting in similar amplitudes among GWF, WUZ, and EPD at 70 hPa. The Differences in the phase relation among all of the simulations were very small.
The vertical advection WUZ consists of two terms. One is associated with low-frequency (<< the QBO frequency) vertical wind, and the other with the vertical wind from the secondary circulation of the QBO (e.g., Plumb & Bell, 1982). This is because WUZ is the product of the residual vertical wind and the vertical shear of zonal wind, both of which have two major peaks in power spectrum at the QBO and annual frequencies below the middle stratosphere, (1) where the superscripts Q and L represent the QBO and the low-frequency components, respectively. The first term on the right-hand side of the approximation sign represents the contribution of the lower-frequency (background) vertical wind, corresponding to the interaction with BDC, and the second term is the interaction between the background zonal wind shear and the QBO secondary circulation. Henceforth the first term is referred to as BDT and the second term as SCT. Just below the maximum power altitude at 30 hPa, the magnitude of each term on the right hand side of Eq. (1)  where |(X) Q | represents the QBO amplitude of X.
It is evident that the residual vertical velocity is one order smaller than the vertical shear of zonal wind for both the QBO and the low-frequency components. Because the QBO vertical shear of zonal wind is almost of opposite phase to the QBO vertical wind (e.g., Plumb & Bell, 1982;Andrews et al., 1987), and the low-frequency term of vertical shear of zonal wind and that of vertical wind are negative-and positive-definite, respectively, at 30 hPa throughout the whole period from 1960 to 2100, BDT and SCT exert acceleration in approximately the same direction with a ratio of about 8 to 1, indicating that BDT is the major acceleration term of the vertical momentum advection. The magnitude of BDT is also much larger than that of SCT at other altitudes, in particular in the lower stratosphere, because the vertical shear is larger in the QBO than in the low-frequency and the residual vertical wind is larger in the low-frequency than in the QBO. However, other features such as acceleration direction and trend differ depending on altitude and latitude. In the lower stratosphere at 50 hPa and below, for example, BDT and SCT do not always exert acceleration in the same direction because the low-frequency vertical shear did not hold but, in the middle of the integration period , changes its sign from negative to positive at 50 hPa and from positive to negative at 70 hPa over the equator as a result of the climate change caused by the GHG increase.

Future change of the zonal wind forcing in the tropical stratosphere
Changes in the zonal wind forcings affect the background zonal wind, which in turn modifies the QBO through the interactions between the background zonal wind and the QBO; thus, such changes are one of the crucial factors affecting the future QBO. Figures 12-14 depict annual mean GWF, WUZ, and EPD in the past (30 years, 1970-1999) and their changes after 100 years (30 years, 2070-2099) in the low-latitude (30°S-30°N) stratosphere from 70 to 5 hPa. GWF and EPD generally exert negative (westward) acceleration in this domain, except in some regions: namely, in the upper stratosphere around 15° in both hemispheres, and at 50 hPa in the deep-tropics (10°S-10°N) for GWF, and at 30 hPa south of the equator for EPD. On the other hand, WUZ exerts positive (eastward) acceleration almost everywhere in this domain.
It is evident that the distributions of the forcing magnitude are commonly very different between the mid-tropics (15°S-15°N) and higher latitudes. In the mid-tropics, the magnitude of GWF is minimal below 50 hPa, and above 50 hPa it is maximal and more concentrated south of the equator. The magnitude of WUZ shows a similar pattern with a broader maximum above 50 hPa. The magnitude of EPD is also minimal in the mid-tropics, and the width of the minimum area becomes narrower with altitude; the magnitude is largest at 50 hPa, and then the minimum area becomes more concentrated downward toward the equator. The EPD minimum over the tropics reflects the suppression of the propagation of Rossby waves from the mid-latitudes and subtropics against the easterly background zonal wind (Fig. 3).    Fig.12 except for EPD. Units are 10 -2 ms -1 per day. The contour interval is 2 for the values greater than -12, except that the outer contours in the left panels have the values -16, -32, and -64. The contour interval is 1 for values greater -2, but for values less than -2, the contour interval is 2 in the right panels.
Future changes in GWF and EPD are larger outside of the mid-tropics below the middle stratosphere in GHG-increase simulations, owing to the shift of the zero line of zonal wind (e.g., Deushi & Shibata, 2011). In contrast, future GWF and EPD changes are small in the mid-tropical middle stratosphere, in particular at around 30 hPa, in GHG-increase simulations. On the other hand, in the mid-tropical middle stratosphere, the future WUZ changes are larger negative values with a large cell structure. In the no-climate-change simulation, the future changes in the QBO forcings are very small everywhere. Figure 15 exhibits the QBO zonal wind (averaged from 10°S to 10°N) at 30 hPa for the whole period from 1960 to 2099 in simulations for B21, Sb, and Sc. It is evident that the QBO zonal winds show long-term decreasing trends in all of the GHG-increase simulations, whereas in the no-climate-change simulation Sc, the QBO zonal wind shows no long-term trend at all. To quantify the trends, the QBO amplitude was calculated by two independent methods, that is, direct and wavelet methods. The direct method assumes that the QBO amplitude is √2σ, analogous to a monochromatic wave, where σ is the root mean square of the QBO time series over three cycles, even though the QBO time series is not a single sinusoidal wave (Pascoe et al., 2005;. The cycles are defined as periods from one zero point to the third consecutive zero point, and the QBO amplitude is assigned to the center time of the three cycles. In this way, the QBO amplitudes are first calculated discretely in time, about a half (14 months) cycle apart, and then monthly amplitudes are obtained through cubic interpolation with Lagrange polynomials. This procedure results in a QBO amplitude that is a low-passed (> 3 cycles~7 years) data with a temporal resolution of about 14 months. The wavelet method used a Morlet mother wavelet (plane wave modified by a Gaussian envelope) with non-dimensional frequency ω 0 =6 (e.g., Torrence & Compo, 1998), so the result of the wavelet analysis incorporates a mean information within approximately three cycles centered at time and frequency concerned.   Different from the distinct decreasing trend of the QBO amplitude, the QBO period does not exhibit any significant trend as depicted by the local wavelet power spectrum of the QBO zonal wind at 20 hPa, where the power maximizes (Fig. 17). The QBO powers are mainly confined to between 22 and 32 months, and the peak power frequencies (periods) stabilize at 28 months, in spite of some excursions toward 20-24 months in the GHG-increase simulations. Similarly to the direct evaluation, the QBO powers show a common decreasing trend in the GHG-increase simulations. The no-climate-change simulation Sc, on the other hand, shows a very stable QBO without a substantial trend in either frequency or power magnitude.   Figure 18 exhibits the evolution of the QBO local wavelet power from 70 to 10 hPa for B21, Sb, and Sc simulations. Below 30 hPa, the QBO powers show an apparent decreasing trend in the GHG-increase simulations, whereas they do not at 10 hPa. The no-climate-change simulation Sc, on the other hand, does not seem to show a discernible trend for the QBO zonal wind at any altitude. To quantify linear trends in the QBO amplitude, the linear regression is calculated not for the monthly interpolated amplitude but for the original about half-cycle (~14 months) interval amplitude using the least squares method at each altitude. Figure 19 depicts the vertical profile of the linear trend in the QBO zonal wind amplitude from 70 to 10 hPa for all simulations. The GHG-increase simulations all exhibit a statistically significant (at the 99% confidence level) negative trend below 20 hPa, except for that of B23 at 20 hPa. The maximum negative trends at 30 hPa are about 0.3 ms -1 per decade, and negative trend below and above 30 hPa are 0.2-0.1 ms -1 per decade. In the upper stratosphere, the trends are near zero or small positive, but the differences among the simulations are not small. The no-climate-change simulation exhibits a near-zero trend in the lower stratosphere below 50 hPa, and a small but statistically significant positive trend above 30 hPa. These results demonstrate that the GHG increase is responsible for the decrease in the QBO amplitude below 20 hPa and that ODS changes alone cause the increase above 30 hPa. Fig. 19. Trends of the QBO amplitude. Vertical profiles of the linear trends (ms -1 per decade) in the QBO amplitude for the zonal-mean zonal wind (averaged from 10°S to 10°N) from 70 to 5 hPa in all simulations. Open circles and diamonds represent the trend being different from zero at the 99% and 95% confidence levels (t-test), respectively. The QBO amplitude was calculated by the direct method.

Trends in the forcings of the QBO zonal wind
The QBO trends of the zonal wind amplitude are statistically significant below 20 hPa in the GHG-increase simulations and above 30 hPa in the no-climate-change simulation; therefore, the forcings of the QBO zonal wind themselves show similar trends. Figure 20 exhibits the evolutions of the QBO forcing terms, GWF, WUZ, and EPD, at 30 hPa in all simulations. As shown by the composite relations among the QBO and its forcings (Fig. 11), GWF has the largest amplitude and EPD the smallest amplitude during the whole period .
Although the short-term variations are considerable, it is evident that the forcings show decreasing trends in the GHG-increase simulations and that they do not in the no-climatechange simulation. Similarly to the quantification of the trend of the QBO zonal wind amplitude, linear regression was used to calculate the trends of the forcings. In addition, to include different responses at different frequencies, that is, a larger response at a smaller frequency, even in the QBO narrow-frequency domain (period of 20-40 months) and to get more direct contribution to the trend of the QBO zonal wind, each raw forcing was integrated over time and band-passed to yield the corresponding zonal wind response (see Eq. (1)). For example, u GWF is the integrated QBO zonal wind (response wind) due to GWF. Following this procedure, the linear regressions were calculated. Note that the trends in the forcing amplitudes are not additive, even though forcings are additive in the temporal space, owing to their different phases (Fig. 11).  The GHG-increase simulations exhibit the same statistically significant negative trend at 30 ant 70 hPa and, in most cases, at 50 hPa as well, with a maximum of about 0.2x10 -2 ms -1 day -1 per decade for GWF at 70 hPa. At 20 hPa the trends are dispersed around zero, and at 10 and 5 hPa, the GWF shows statistically significant positive values of nearly the same magnitude. The trends of integrated zonal wind u GWF are qualitatively similar to those of GWF but the statistically significances are not necessarily so. For example, while most GHG-increase simulations show a statistical significant negative trend at 70 hPa and positive trend above 10 hPa both in u GWF and in GWF, the number of statistically significant trends decreases at 50 and 30 hPa and increases at 20 hPa in u GWF compared with GWF. The no-climate-change simulation exhibits statistically significant positive trends of about 0.1x10 -2 ms -1 day -1 per decade at 20 and 30 hPa, and a negative trend of similar magnitude at 5 hPa for GWF.
The vertical transport of zonal momentum WUZ shows a statistically significant trend at almost all altitudes in the GHG-increase simulations (Fig. 22). A negative trend of about 0.1x10 -2 ms -1 day -1 per decade is seen at 30 and 50 hPa, and a positive trend of nearly the same magnitude at 70 hPa and above 20 hPa. The no-climate-change simulation exhibits almost no trends at any altitude in both WUZ and u WUZ , and most apparent trends are thus not statistically significant. The resolved-wave forcing EPD and u EPD exhibit a vertically flat structure (Fig. 23), and the trend magnitudes are smaller than those of the other forcings in all simulations; nevertheless, most of the trends are statistically significant.

Discussion
The GHG-increase simulations have thus far demonstrated that the QBO amplitude decreases linearly below 20 hPa with rates of 0.1-0.3 m s -1 per decade and increases above 10 hPa with rates of 0-0.1 m s -1 per decade (Fig. 19). As shown by the difference between the B2 and Sb simulations, the trends of the B2 members tend to be smaller below 20 hPa than the Sb trend, but the difference between them is not substantial. That is, the inter-member variations are too large to ascribe the difference between the B2 members and Sb to the forcing difference, that is, ODSs. On the other hand, the small but statistical significant increasing trends of less than 0.1 m s -1 per decade at 30 hPa and above exhibited by the noclimate-change simulation is likely an effect of the ODS decrease, although to confirm this result , an Sc ensemble simulation should be performed.
The cause of the QBO trend should reflect the trends of the forcings, GWF, EPD, and WUZ. However, because the forcings include not only external quantities but also internal ones related to the QBO, separating the forcings from the responses is not straightforward. For example, the propagation of the resolved waves and that of the parameterized gravity waves strongly depend on the background QBO wind, as does the deposition of their momentum on the QBO. The vertical advection WUZ is also affected by the QBO, because WUZ is explicitly expressed as the product of the QBO and the background wind field (Eq. (1)). In addition, the trends in the forcing amplitude are not additive because of the phase differences as stated before; thus, evaluation of their total effects is not straightforward either, if the signs of the forcing trends are not the same. Almost all of the trends in the forcings and their response winds in the GHG-increase simulations exhibit negative trends at 30 hPa, where the trend of the QBO wind is the largest negative; therefore, the forcings and responses are qualitatively consistent. At 50 hPa, the magnitude of the negative trend in WUZ is significantly larger than that of the GWF (small negative) or the EPD trend (nearly zero or small positive), which is also consistent with the negative trend in the QBO wind. At 20 hPa, the trends are on average nearly zero in GWF, have small positive values in WUZ, and small negative values in EPD, as are those of their responses, which is not consistent with the small but significant negative trend in the QBO wind. At 70 hPa, GWF exhibits a statistically significant negative trend, and its magnitude is largest compared with the magnitudes of the trends at other altitudes. WUZ, on the other hand, exhibits a statistically significant positive trend, the magnitude of which is about three-fourths that of GWF. EPD tends to have a small negative trend, but its response is nearly zero in the ensemble mean, indicating that EPD plays a minor role at 70 hPa. WUZ trends in the GHG-increase simulations are statistically significant below the middle stratosphere but they show sharp contrasts in sign, being negative at 70 hPa and positive at 50 and 30 hPa, although their magnitudes are almost the same at 0.1x10 -2 ms -1 day -1 per decade (Fig. 22). To explore the cause of this, we evaluated the trend of each component, that is, BDT and SCT, of WUZ (Eq. 1). Figure 24 depicts the evolution of the residual vertical velocities of the low-frequency and the QBO at 30 hPa in simulations for B21, Sb, and Sc, and Fig. 25 displays the evolution of the vertical shears. The tendencies of the residual vertical velocity and the vertical shear at 50 hPa are very similar to those at 30 hPa, aside from the background levels. Figures 26 and 27  decreasing trend comparable to that of SCT. At 50 hPa the tendencies of the four terms are qualitatively similar to those at 30 hPa, but the SCT tendency is much smaller than the BDT tendency, in which the decrease in (/) Q uz   is much larger than the increase in *L w . As a result, BDT shows a slightly larger decreasing trend at 50 hPa than at 30 hPa.
At 70 hPa, the BDT components are a rapidly increasing *L w and a very slowly decreasing (/) Q uz  , leading to an increasing trend in BDT. The SCT components are both increasing but very small, so SCT is much smaller than BDT. At 20 hPa, both SCT components are increasing, so SCT has a positive trend. Although (/) Q uz   is decreasing, BDT has a positive trend comparable to the SCT trend. These results indicate that the BDT trend is of similar magnitude to the SCT trend at 30 and 20 hPa, while the SCT trend is much smaller than the BDT trend at 50 and 70 hPa. This altitude dependence of SCT derives mainly from the low-frequency vertical shear, which is very small as well as its tendency in the lower stratosphere.
In this study the source strength of GWF was fixed in all simulations, irrespective of the ODS evolution and/or the GHG-induced global warming, wherein convective precipitation increased (e.g., Yukimoto et al., 2006;IPCC, 2007) as manifested in the strong warming in the tropical upper troposphere (Fig. 4). This fixation of source strength is different from the precedent simulations of equilibrium CO 2 doubling. Giorgetta & Doege (2005) varied the GWF source by 0% to 20% and showed that a statistically significant reduction in the QBO period occurred in runs with a 10% or 20% increase. Kawatani et al. (2011) demonstrated that wave momentum fluxes increase by 10%-15% in the equatorial lower stratosphere, though those relevant to the QBO forcing do not increase as much. Since tropical convective precipitation is closely linked to gravity and equatorial wave generation (e.g., Horinouchi et al., 2003), we assume that the fixation of the GWF source differs more or less from reality. Thus, to specify appropriate changes in the GWF source, observations of relations between the QBO and other independent quantities such as SST/precipitation or outgoing longwave radiation are needed. propagation under El Niño conditions, when the west-east contrast of SST and precipitation drastically changes in the tropical Pacific (e.g., Philander, 1990). However, since the tropical upwelling is also enhanced (Randel et al., 2009) during El Niño conditions, isolation of the SST/precipitation effect is required for evaluation of the GWF source in the troposphere. On the other hand, according to Fischer and Tung (2008), the QBO period does not show any trend but has been stable within 24-32 months for five decades (from 1953 to 2007) despite period variations due mainly to stalling of the westerly phase. These findings indicate that further investigation is required to identify the tropospheric factors affecting the QBO.

Conclusion
To investigate future changes in the QBO in the tropical stratosphere, transient simulations were carried out for the period from 1960 to 2100 with MRI-CCM by inputting observed and projected GHGs and/or ODSs at the surface. SST/sea ice, generated by an MRI coupled ocean-atmosphere model, was also specified so as to be compatible with the GHGs abundances. Three types of simulations were performed to evaluate the effects of a GHGs increase and an ODSs decrease, both in combination and separately: The first type used the REF-B2 scenario of CCMVal-2, in which both GHGs and ODSs evolve following observed values for past simulations, and specified SRES A1B scenario for GHGs and the adjusted A1 scenario for ODSs for future simulations. The second type of simulation used an SCN-B2b scenario, in which GHGs were the same as in REF-B2 but ODSs were fixed at 1960 levels. The third type used an SCN-B2c scenario, in which ODSs were the same as in REF-B2 but GHGs/SST were fixed at 1960 levels. The REF-B2 simulation was carried out with three members, and the other two were performed with a single member.
The future climate change due to increasing GHGs was characterized as tropospheric warming as a result of SST warming (indirect CO 2 effect), with maximum warming in the tropical upper troposphere, and cooling in the stratosphere (direct CO 2 effect), with maximum cooling in the upper stratosphere. Zonal wind changes were characterized by westerly wind intensification from the upper region of the subtropical jet upward along the equatorward flank of the polar night jet axis, reaching a maximum at around 70 hPa and 40° in both hemispheres. In the SH, the westerly wind intensification also extends down to the surface at around 55°S, resulting in poleward spreading of the subtropical jet. In the tropics, westerly wind strengthening (eastward acceleration) occurred above 30 hPa and weakening (westward acceleration) at 50 hPa, below which very weak strengthening occurred down to the surface. The Brewer-Dobson circulation was similarly intensified in all of the GHGincrease simulations.
The future climate change in the no-climate-change simulation was characterized by strong warming in the upper stratosphere above about 5 hPa in both hemispheres, except in the southern high latitudes, where prominent warming occurred in the lower stratosphere at 300-30 hPa, centered at 100 hPa, with cooling above 30 hPa. This pattern of temperature change was due solely to the ODSs decrease, or, equivalently, ozone recovery. The zonal wind change characterized by slight intensification of the polar night jet in the poleward flank above the middle stratosphere in the NH. In the SH, weakening of the polar night jet and both weakening of the poleward and strengthening of the equatorward flank of the subtropical jet down to the surface occurred. In the tropics, there was almost no change, including in the Brewer-Dobson circulation.
The QBO was spontaneously generated with a period of about 28-month period in the equatorial stratosphere by MRI-CCM under both past and current conditions in all simulations. The QBO frequency was very stable without any substantial trend in all simulations, and the QBO amplitude showed a small but statistically significant decreasing trend below 20 hPa with a maximum (~0.3 ms -1 per decade) at 30 hPa in GHG-increase simulations. ODS decreases in the future (i.e., ozone recovery) led to very small increasing trend (<0.1 ms -1 per decade) in the QBO amplitude above 20 hPa. The causes of the QBO amplitude trend were largely changes in the parameterized gravity wave forcing and the vertical advection of zonal wind momentum. In the QBO trend, the vertical advection was brought about not only by the Brewer-Dobson circulation but also by the QBO secondary circulation above the middle stratosphere, while the Brewer-Dobson circulation was dominant in the lower stratosphere, mainly because of the very weak vertical shear and the very small vertical shear trend of the background wind.
The QBO effects globally extend from the stratosphere to the troposphere, and thus the projected decreasing trend of the QBO amplitude in the CCM simulations would likely cause, more or less, changes in various phenomena in the whole atmosphere. However, since the current CCM simulations included only forcings of SST, GHGs, and ODSs, other forcings such as solar irradiance variability and volcanic aerosols would induce different responses in the QBO, depending on the timescale of these forcings. In addition, different scenarios in GHGs result in different SSTs, and thereby likely in different QBO trends. Therefore, there are larger uncertainties in the effects of the QBO in the future than in the QBO itself.