Multitimescale variations in modeled stratospheric water vapor derived from three modern reanalysis products

Abstract

Stratospheric water vapor (SWV) plays important roles in the radiation budget and ozone chemistry and is a valuable tracer for understanding stratospheric transport. Meteorological reanalyses provide variables necessary for simulating this transport; however, even recent reanalyses are subject to substantial uncertainties, especially in the stratosphere. It is therefore necessary to evaluate the consistency among SWV distributions simulated using different input reanalysis products. In this study, we evaluate the representation of SWV and its variations on multiple timescales using simulations over the period 1980–2013. Our simulations are based on the Chemical Lagrangian Model of the Stratosphere (CLaMS) driven by horizontal winds and diabatic heating rates from three recent reanalyses: ERA-Interim, JRA-55 and MERRA-2. We present an intercomparison among these model results and observationally based estimates using a multiple linear regression method to study the annual cycle (AC), the quasi-biennial oscillation (QBO), and longer-term variability in monthly zonal-mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ mixing ratios forced by variations in the El Niño–Southern Oscillation (ENSO) and the volcanic aerosol burden. We find reasonable consistency among simulations of the distribution and variability in SWV with respect to the AC and QBO. However, the amplitudes of both signals are systematically weaker in the lower and middle stratosphere when CLaMS is driven by MERRA-2 than when it is driven by ERA-Interim or JRA-55. This difference is primarily attributable to relatively slow tropical upwelling in the lower stratosphere in simulations based on MERRA-2. Two possible contributors to the slow tropical upwelling in the lower stratosphere are suggested to be the large long-wave cloud radiative effect and the unique assimilation process in MERRA-2. The impacts of ENSO and volcanic aerosol on ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry variability are qualitatively consistent among the three simulations despite differences of 50 %–100 % in the magnitudes. Trends show larger discrepancies among the three simulations. CLaMS driven by ERA-Interim produces a neutral to slightly positive trend in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values over 1980–2013 ($+\mathrm{0.01}$ ppmv decade${}^{-\mathrm{1}}$), while both CLaMS driven by JRA-55 and CLaMS driven by MERRA-2 produce negative trends but with significantly different magnitudes ($-\mathrm{0.22}$ and $-\mathrm{0.08}$ ppmv decade${}^{-\mathrm{1}}$, respectively).

1 Introduction

Water vapor is one of the most influential greenhouse gases , modulating not only the surface radiative forcing  but also stratospheric ozone loss e.g.,. The extreme dryness of the stratosphere results from “freeze-drying” of air entering the stratosphere, as initially explained in relation to the Brewer–Dobson circulation BDC;. Mixing ratios of water vapor in the lower stratosphere are extremely low as a result, but nonetheless vary substantially in time and space.

Substantial inconsistencies among observational estimates of stratospheric water vapor (SWV) remains, including both balloon-based and satellite-based instruments e.g.,. Artificial jumps introduced by changes in the observing system influence the reanalysis temperatures through the assimilation, which affect the reanalysis water vapor. The inconsistencies are even more pronounced among atmospheric reanalyses and independent observations . Uncertainties in balloon- or aircraft-based observations arise primarily due to insufficient spatiotemporal coverage  and the unreliability of operational sondes at stratospheric altitudes e.g.,. Inconsistencies in satellite observations of SWV reflect limited observational periods and short overlap times , which make it difficult to control for platform-specific biases or differences in temporal or spatial sampling patterns . Observations of SWV are rarely assimilated in reanalysis systems, for which estimates of SWV are effectively model products, in some cases nudged to climatologies . Differences in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ observations in the upper troposphere and stratosphere and the unreliability of reanalysis estimates of SWV have motivated several efforts to merge observational records from different satellites . Such homogenized records facilitate analysis of long-term changes in SWV across the most recent 2–3 decades.

Chemical transport models (CTMs) provide another approach to understanding the multitimescale variability and global distribution of SWV e.g.,. Water vapor values entering the stratosphere are determined primarily by the lowest temperature along their advective transport pathway in the tropical tropopause layer (TTL) . Lagrangian approaches provide accurate records of the temperature histories of air parcels that Eulerian models cannot provide, and therefore provide more reliable representations of entry mixing ratios in SWV e.g.,. This distinction is particularly relevant around the tropopause, where temperature gradients are both large and highly variable . In addition to transport across the tropical tropopause, photochemical oxidation of methane (${\mathrm{CH}}_{\mathrm{4}}$) is an important source of SWV, especially in the middle and upper stratosphere. Previous studies have concluded that recent increases in ${\mathrm{CH}}_{\mathrm{4}}$ have substantially contributed to long-term variability in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the stratosphere e.g.,. In this study, we use a forward Lagrangian transport model with implanted methane chemistry to study the climatological features of SWV, and compare the results against observational estimates from the Stratospheric Water and OzOne Satellite Homogenized (SWOOSH) dataset version 2.5  and the Aura Microwave Limb Sounder (MLS) version 4 .

Key aspects affecting the performance of CTMs with respect to SWV are the meteorological fields selected to drive the transport and freeze-drying in TTL. Modern meteorological reanalyses, as “best estimates” of the historical evolution of the atmospheric state , are commonly used to provide the necessary meteorological variables for CTMs. Current widely used reanalysis products include the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis ERA-I;, the Japanese 55 year Reanalysis JRA-55;, the Modern-Era Retrospective Analysis for Research and Applications version 2 MERRA-2;, and the Climate Forecast System Reanalysis CFSR;. Substantial uncertainties among these reanalyses have been identified by previous studies, including significant differences in temperature and wind structures e.g.,, diabatic heating rates in the tropical upper troposphere and lower stratosphere (UTLS) e.g.,, and representations of the BDC e.g.,. Such differences among reanalysis products are critical for simulations of atmospheric composition but have rarely been discussed in this context.

evaluated SWV simulated using the same trajectory model with meteorological fields taken from MERRA the predecessor of MERRA-2;, ERA-I and CFSR. The trajectory model based on ERA-I produced a low bias of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the lower stratosphere and a “tape-recorder” signal  that was 30 % too fast in the 17–22 km layer. By contrast, use of MERRA resulted in reasonable ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values but a tape-recorder signal that was 15 % too slow within a similar altitude range, while use of CFSR resulted in a wet bias in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values but a reasonable propagation of the tape-recorder signal. An earlier study by reported similar biases in simulated water vapor entering the stratosphere via the Asian monsoon, with trajectories based on MERRA indicating systematically larger entry mixing ratios relative to trajectories based on ERA-I. Further evaluation of how uncertainties in current reanalysis products impact simulations of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the stratosphere is lacking, particularly with respect to the more recently released JRA-55 and MERRA-2 products. The former is currently the most recent full-input reanalysis (which assimilate the full observing system) to provide coverage of the pre-satellite era, while the model used for the latter has been reported to have a more realistic spontaneous quasi-biennial oscillation (QBO) than its predecessor MERRA .

In this study, we provide an intercomparison of SWV produced using a Lagrangian transport model driven by three recent reanalysis products: ERA-I, MERRA-2 and JRA-55. We also present comparisons against SWV estimates from satellite observations. Our evaluations of simulated SWV focus on the climatological annual mean, the annual cycle (AC) and the QBO. We also discuss some key sources of variability, including the El Niño–Southern Oscillation (ENSO) and variations in volcanic aerosol, as well as long-term linear trends in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios. The objective of this work is to evaluate the sensitivity of simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ variability to uncertainties in the driving meteorological fields among recent reanalysis products. The study also sheds light on the quality of reanalysis products in the stratosphere, especially their representations of dynamical fields (e.g., winds, heating rates and temperatures). Few independent observational datasets are available to evaluate these dynamical variables in the upper troposphere and stratosphere, especially with respect to their influences on SWV.

2 Data and methods

2.1 CLaMS simulations

Our simulations of SWV cover the period 1980–2013 using the Chemical Lagrangian Model of the Stratosphere CLaMS;. Modeled water vapor mixing ratios entering the stratosphere are based on a simplified dehydration scheme as described in detail by and . The saturation volume mixing ratio is calculated for each trajectory, where the saturation pressure ($p_{\mathrm{s}}$) is given by $p_{\mathrm{s}}={\mathrm{10}}^{-\mathrm{2663.5}/T+\mathrm{12.537}}$ . If saturation occurs along a CLaMS air parcel trajectory, the amount of water vapor in excess of the critical saturation mixing ratio (100 % with respect to ice) is instantaneously transformed to the ice phase. The sedimentation fall speed is calculated by assuming a mean ice particle radius. Ice sedimentation is then determined by comparison of the sedimentation length over the model time step against a characteristic length ($\sim \mathrm{300}$ m) . Furthermore, if ice exists in a subsaturated parcel, this ice is instantaneously evaporated until saturation is reached.

Temperatures and horizontal winds are prescribed from the ERA-I , JRA-55 and MERRA-2 reanalyses. The cross-isentropic vertical velocity ($\dot{\mathit{\theta }}$) is derived from the total diabatic heating rates, which includes all-sky radiation, latent heat release and the diffusive turbulent heat transport from each reanalysis product e.g.,. Water vapor mixing ratios at pressures greater than $\sim$500 hPa are set equal to water vapor mixing ratio products from the corresponding reanalysis. The CLaMS runs driven by ERA-I, JRA-55 and MERRA-2 are respectively referred to as CLaMS-ERA, CLaMS-JRA and CLaMS-MRA. Additional information regarding the structures of these reanalysis systems has been provided by . We exclude CFSR because the model and data assimilation system used to produce this reanalysis changed abruptly at the beginning of 2011 . This and other stream transitions produced substantial discontinuities in stratospheric variables based on CFSR during the 1980–2013 analysis period .

We use the critical saturation rate of 100 % with respect to ice, although supersaturation with respect to ice is frequently observed at these altitudes e.g.. Within CLaMS, we have the freedom to “optimize” supersaturation thresholds to produce better agreement with observed ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ values. The effects of tuning the critical supersaturation threshold in CLaMS have an effect similar to that of applying a frost point offset to the Lagrangian dry point temperature in that an increase in the supersaturation threshold enhances both the mean value and the amplitude of the annual cycle in simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$. Due to uncertainty in the appropriate value of the supersaturation threshold, differences among the modeled values of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratio or between modeled values against observations cannot be unequivocally interpreted as errors in Lagrangian dry point temperatures in these reanalyses. Our intercomparison among modeled and observed ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ values is thus limited to uncertainties related to using different reanalyses under the standard configuration of the model.

Besides the freedom in optimization of supersaturation threshold, additional degrees of freedom are in the parameterized small-scale mixing in the model. The modeled water vapor mixing ratios are influenced by parameterized small-scale mixing driven by large-scale deformation . tested the sensitivity of modeled ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratio to mixing parameters in CLaMS and found that the range of mixing strength they considered was representative of actual uncertainty in small-scale mixing. Our study uses the reference configuration of mixing parameters, which has been shown to produce the best agreement with MLS observations by .

The model calculations include methane oxidation, which is a source of water vapor mainly in the middle and upper stratosphere, with concentrations of hydroxyl, atomic oxygen and chlorine radicals taken from model-generated climatologies . We explicitly diagnose the fraction of water vapor supplied by methane oxidation. Methane-supplied water vapor (${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{CH}}_{\mathrm{4}}}$) at any given location is calculated as

1$${\displaystyle }{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{2}\cdot ({\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}-{\mathrm{CH}}_{\mathrm{4}}),$$ where ${\mathrm{CH}}_{\mathrm{4}}$ is modeled methane accounting for loss by oxidation and ${\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}$ is the passively transported methane assuming the same source and circulation without photochemical loss. Reconstructed methane values (${\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}$) are calculated as the mean tropospheric ${\mathrm{CH}}_{\mathrm{4}}$ prescribed at the lower boundary of the model (see Fig.  in Appendix ), lagged by the mean age of air ($\mathrm{\Gamma }$) in each run : 2$${\displaystyle }{\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}(x,t)={\mathrm{CH}}_{\mathrm{4}}^{\mathrm{LB}}(t-\mathrm{\Gamma }(x,t\left)\right).$$ The fractional release factor of ${\mathrm{CH}}_{\mathrm{4}}$ ($\mathit{\alpha }$) is expressed as 3$${\displaystyle }\mathit{\alpha }=({\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}-{\mathrm{CH}}_{\mathrm{4}})/{\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}.$$ Note that showed that chemical loss can increase the time dependence of $\mathit{\alpha }$. Our calculation of $\mathit{\alpha }$ neglects this effect, as also suggested that this influence should be very small for tracers with weak tropospheric trends, like ${\mathrm{CH}}_{\mathrm{4}}$ ($\sim \mathrm{0.2}$ % yr${}^{-\mathrm{1}}$–0.3 % yr${}^{-\mathrm{1}}$).

Stratospheric water vapor without the contribution from methane oxidation (${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{nCH}}_{\mathrm{4}}}$) is then diagnosed as 4$${\displaystyle }{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{nCH}}_{\mathrm{4}}}={\mathrm{H}}_{\mathrm{2}}\mathrm{O}-{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{CH}}_{\mathrm{4}}}.$$

Our intercomparison of the three reanalysis-driven CLaMS runs covers the period from January 1980 through December 2013 (hereafter referred to as the “CLaMS period”) within the vertical range from $\mathit{\theta }=\mathrm{350}$ K to $\mathit{\theta }=\mathrm{2000}$ K. A reference monthly-mean SWV is calculated by averaging the three reanalysis-driven runs (hereafter referred to as the multireanalysis mean or MRM). The MRM is used as a benchmark to more effectively highlight differences among the three model runs.

2.2 Observational estimates

The Stratospheric Water and Ozone Satellite Homogenized (SWOOSH) database provided by the NOAA Chemical Sciences Division (CSD) contains vertically resolved ozone and water vapor data from the SAGE-II/III, UARS HALOE, UARS MLS, and Aura MLS satellite instruments starting from 1984 . We use the SWOOSH version 2.5 zonal-mean monthly-mean time series of merged water vapor mixing ratios with 2.5${}^{\circ }$ resolution on 21 isentropic levels from 300 to 650 K. The homogenization process, which has been described by , is designed to minimize artificial jumps in time and account for intersatellite biases. The merged SWOOSH data thus provide a long-term time series with reliable representations of interannual to decadal variability.

We compare monthly mean water vapor from CLaMS runs with SWOOSH water vapor on the same vertical grid (21 $\mathit{\theta }$ levels from 350 to 650 K). For each latitude–$\mathit{\theta }$ grid location, comparison between CLaMS and SWOOSH starts from the first month when SWOOSH has more than 12 months of available ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ data within the following 2 years and ends in December 2013. In the following, we refer to this period as the “SWOOSH period”. Note that while CLaMS provides continuous temporal coverage, the SWOOSH data may include some gaps. Additional details regarding SWOOSH data coverage are provided in Appendix .

In addition to SWOOSH, we use Aura MLS version 4 retrievals of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$  for comparison with CLaMS simulations during the period 2005–2013 (the “MLS period”). MLS provides over 3000 profiles per day, with water vapor estimates at 30 pressure levels from 316 to 1 hPa. The vertical resolution in the stratosphere is approximately 3 km (2.5–3.5 km). Uncertainties in the water vapor retrievals are on the order of $\sim \mathrm{10}$ % in the lower stratosphere and $\sim \mathrm{5}$ % in the upper stratosphere. The relatively high frequency of horizontal sampling and high quality of vertical profiles allows MLS ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ to reliably represent the zonally and monthly-averaged distribution of SWV. We interpolate MLS ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ profiles to 26 isentropic levels, chosen to span the range 350–2000 K at a vertical resolution close to that of the original retrievals.

The SWOOSH dataset is also based in part on Aura MLS version 4 retrievals, particularly during the MLS period. Thus, the comparisons of CLaMS against SWOOSH and comparisons of CLaMS against MLS during the MLS period are not independent. Differences between the comparisons are due to the homogenization procedure applied in SWOOSH. We did not apply the MLS averaging kernels to CLaMS ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ simulations. It is because the application of MLS averaging kernels to CLaMS ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ simulations can potentially produce artifacts, especially at high latitudes more details provided bytheir Fig. 2. Our main focus is on the comparison among the reanalyses and to apply averaging kernels that could smear out the differences among the CLaMS runs. However, it is worth noting that artifacts may appear in MLS observations, especially in high latitudes, when MLS retrievals are compared to model results without MLS weighting functions.

2.3 Extraction of variability at multiple timescales

The objective of this study is to examine and compare the climatological features of variability in model-based and observationally based estimates of SWV. Pronounced periodic and quasi-periodic signals in SWV include the annual cycle (AC), the quasi-biennial oscillation (QBO) and semiannual harmonic (SAH). The AC, QBO and SAH are therefore considered explicitly in our regression model:

5$${\displaystyle }\mathit{\chi }\left(t\right)=\overline{\mathit{\chi }}+{\mathit{\chi }}_{\mathrm{AC}}\left(t\right)+{\mathit{\chi }}_{\mathrm{SAH}}\left(t\right)+{\mathit{\chi }}_{\mathrm{QBO}}\left(t\right)+{\mathit{\chi }}_{\mathrm{res}}\left(t\right).$$ After determining the mean value calculated over the whole considered time series ($\overline{\mathit{\chi }}$), we extract AC (${\mathit{\chi }}_{\mathrm{AC}}\left(t\right)$), SAH (${\mathit{\chi }}_{\mathrm{SAH}}\left(t\right)$) and QBO (${\mathit{\chi }}_{\mathrm{QBO}}\left(t\right)$) signals following the harmonic regression method used by . In this method, one sine and cosine pair are used for each regression (e.g., AC and SAH). The zonal wind records at 50 and 30 hPa over Singapore provided by Freie Universität Berlin (FUB) are quasi-orthogonal, which are used for the QBO regression. The AC, QBO and SAH variables are assessed from two aspects of periodic variations: amplitudes ($A_{\mathrm{AC}}$ and $A_{\mathrm{QBO}}$ shown in the first half of Figs.  and , $A_{\mathrm{SAH}}$ in top of Fig. ) and phases ($P_{\mathrm{AC}}$ and $P_{\mathrm{QBO}}$ shown in the second-half of Figs.  and , $P_{\mathrm{SAH}}$ not shown). Note that (1) the amplitude represents half the corresponding variation from maximum to minimum, (2) $P_{\mathrm{AC}}$ and $P_{\mathrm{SAH}}$ are defined as the month of the annual maximum in the regression fit and (3) $P_{\mathrm{QBO}}$ is defined as the month (within 0–28) with the largest lag-correlations between the QBO fit and the 50 hPa Singapore wind.

The linear trend ($C_{\mathrm{trd}}$) is estimated by applying a least-squares linear fit to the residuum of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ variability (${\mathit{\chi }}_{\mathrm{res}}\left(t\right)$ in Eq. ) after removing the AC, QBO and SAH signals. Note that the linear trend for ${\mathit{\chi }}_{\mathrm{res}}\left(t\right)$ is identical to the trend of the original ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ time series since periodic signals such as AC have zero long-term trend when the full cycle is taken. The quasi-periodic signal like QBO does not show long-term trend in the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry over the considered period. Additional variability after subtracting the AC, QBO and SAH signals may result from the influences of ENSO e.g.,, and variations in stratospheric aerosol e.g.,. The contributions of such variations to ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios are discussed in Sect. .

Figure 1

Multiple linear regression results for monthly-averaged tropical mean (20${}^{\circ }$ S–20${}^{\circ }$ N) water vapor on the $\mathit{\theta }=\mathrm{400}$ K isentropic surface from SWOOSH. The top panel (a) shows anomalies of SWOOSH data relative to the mean value over the 1993–2013 period (red dots) and the regression fit (black line) based on Eq. () (correlation coefficient $r=\mathrm{0.91}$). The AC, QBO and SAH components of the regression are shown in panels b, c and d, respectively. Additional variability beyond the three periodic and quasi-periodic components (residuum of the regression) is shown in the bottom panel (e). The linear trend in this time series is also shown with $\pm \mathrm{2}\mathit{\sigma }$ uncertainty shown as gray shading. Note the different $y$-axis ranges across the different panels. The amplitudes and phases of the AC, QBO and SAH signals are listed in the upper-right corners of the corresponding panels. The standard deviations corresponding to each contribution are listed near the bottom of each panel, together with correlation coefficients for each term against the original time series. The standard deviation of the original SWOOSH time series (${\mathit{\sigma }}_{\mathrm{org}}$) is 0.73 ppmv.

[Figure omitted. See PDF]

Figure  provides example regression results for SWV on the 400 K isentropic surface in the tropics (20${}^{\circ }$ S–20${}^{\circ }$ N) based on the SWOOSH monthly merged satellite dataset. This time series approximates that of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ values entering the stratosphere (${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$) at the base of the “tropical pipe” , and is characterized by a climatological annual mean of 3.83 ppmv and a negative trend of 0.22 ppmv decade${}^{-\mathrm{1}}$. The regression model explains over 80 % of the variance in ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$. According to the correlation coefficient of each variable, the periodic terms listed in descending order of influence are AC, QBO and then SAH. The AC represents the most pronounced variation in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ at this level, with a correlation coefficient of 0.88, an amplitude of $A_{\mathrm{AC}}=\mathrm{0.89}$ ppmv and an annual maximum in $P_{\mathrm{AC}}=$ October. The QBO is also a significant factor, with an amplitude of 0.21 ppmv and a maximum lag-correlation (correlation coefficient $=\mathrm{0.16}$) against the 50 hPa Singapore winds of $P_{\mathrm{QBO}}=\mathrm{1}$ month. The bottom panel of Fig.  shows other variability in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios (the residuum of the regression in Eq. ). This residual variability represents a substantial component of the total variability, and contains intraseasonal variability, interannual variability and the long-term trend. Features of the AC, QBO and the other variability are compared and discussed in more detail in Sects.  to , respectively. Since the contribution of SAH to the total variance is generally smaller than that of other variables, we defer discussion of the SAH influence to Appendix .

3 Climatological annual mean of stratospheric water vapor

Figure  shows climatological annual-mean simulated SWV from the three CLaMS runs, as well as the components of SWV with and without ${\mathrm{CH}}_{\mathrm{4}}$ oxidation (${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{nCH}}_{\mathrm{4}}}$). The climatological mean is based on the period 1985–2013 because ${\mathrm{CH}}_{\mathrm{4}}^{\mathrm{rec}}$ cannot be diagnosed during the first few years of the simulation. Zonal-mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ matches zonal-mean ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{nCH}}_{\mathrm{4}}}$ well in the lower stratosphere but consists mainly of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{CH}}_{\mathrm{4}}}$ in the upper stratosphere.

Figure 2

Comparison of climatological annual mean SWV during the CLaMS period (1985–2013). Top row (a1–c1): annual-mean zonal-mean simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ based on (left to right) CLaMS-ERA, CLaMS-JRA and CLaMS-MRA. Dashed lines mark the zonal-mean locations of isobaric surfaces. Middle row (a2–c2): as in (a1–c1) but for simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ without methane oxidation. These distributions highlight the effects of variations in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values (tropopause temperature). Black contours show zonal-mean temperatures of 192.5, 195 and 200 K (from thin to thick lines) near the tropical tropopause. Bottom row (a3–c3): as in (a1–c1) but for simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ from ${\mathrm{CH}}_{\mathrm{4}}$ oxidation. The distributions reflect the effect of the BDC on water vapor produced by methane oxidation. Black contours show the fractional ${\mathrm{CH}}_{\mathrm{4}}$ oxidation ratio ($\mathit{\alpha }$).

[Figure omitted. See PDF]

The driest stratosphere is simulated by CLaMS-ERA, while the wettest is simulated by CLaMS-JRA. This difference is mainly attributable to ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{{\mathrm{nCH}}_{\mathrm{4}}}$, as shown in Fig. a2–c2, which is in turn controlled primarily by entry mixing ratios at the tropical tropopause (with the notable exception of dehydration in the Southern Hemisphere polar vortex). These differences reflect differences in tropical tropopause temperatures among the reanalyses (black contours in Fig. 2a2, b2 and c2; see also Fig. ). Tropical tropopause temperatures differ by about 1 K between ERA-I and JRA-55, consistent with a difference of $\sim \mathrm{0.6}$ ppmv in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ mixing ratios assuming the same 100 % saturation threshold applied in the CLaMS model.

Figure 3

Comparison of climatological annual-mean zonal-mean age-of-air (AoA) from the CLaMS-ERA, CLaMS-JRA and CLaMS-MRA simulations (a1–c1), along with corresponding cross-isentropic velocities ($\dot{\mathit{\theta }}$) based on total diabatic heating rates for each reanalysis product (a2–c2). Color shading in the lower panels indicates zonal-mean diabatic upwelling rates that exceed 0.5 K d${}^{-\mathrm{1}}$. The thick gray line shows the zero-contour of diabatic heating rates. Lighter thin gray lines show diabatic downwelling.

[Figure omitted. See PDF]

The amount of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ from ${\mathrm{CH}}_{\mathrm{4}}$ oxidation is related to the time an air parcel has spent in the stratosphere. This time is measured by the stratospheric “age-of-air” (AoA) as shown in Fig. a1, b1 and c1. Values of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ from the CLaMS-MRA simulation are systematically larger than those from the CLaMS-ERA or CLaMS-JRA simulations (see Fig. a3–c3). The most pronounced high biases in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and $\mathit{\alpha }$ in CLaMS-MRA relative to the other two simulations are colocated with the sharpest gradients in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and $\mathit{\alpha }$. Figure  provides useful context for interpreting these differences, as high values of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ are accompanied by systematically older stratospheric air in CLaMS-MRA. Note that consistent climatological mean age differences within the same model framework are presented by . The lower panels of Fig.  show differences in reanalysis heating rates, which are used to drive vertical transport in the CLaMS model. Tropical upwelling in the shallow branch of the BDC is weaker in MERRA-2 than in ERA-I or JRA-55, thus producing larger gradients of AoA and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ between the lower and middle stratosphere. Quantitatively, relative differences in time-mean zonal-mean diabatic heating rates at 400–500 K are as large as 50 % between MERRA-2 and ERA-I, resulting in large differences in the vertical advection of SWV and other tracers. It is worth noting that differences among the reanalyses in the lower stratosphere are larger in the NH subtropics than in the deep tropics or the SH subtropics (not shown).

Although diabatic upwelling in the tropical lower stratosphere is weakest in MERRA-2, upwelling in the tropical middle stratosphere ($\mathit{\theta }$ equal to 600–1000 K) is stronger and covers a wider range of latitudes in MERRA-2 than in ERA-I or JRA-55. Above $\mathit{\theta }=\mathrm{1000}$ K, the deep BDC as represented in JRA-55 is noticeably weaker than that in MERRA-2 or ERA-I. Previous studies have pointed out that the mean magnitudes of tropical upwelling in the lower stratosphere are substantially different among different reanalysis products .

4 Annual cycle

The annual cycle is the largest contributor to total variability in SWV. Moreover, seasonal variations in SWV give essential insight into the overall behavior of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ as these variations reflect the effects of seasonal changes in both the stratospheric circulation and temperatures at the tropical tropopause.

4.1

The ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ tape recorder

Before discussing the global features of the annual cycle in SWV, we first report the simulated climatological seasonalities of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios (${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$) from the three simulations. Values of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ are mainly controlled by “freeze-drying” around the tropical tropopause and the upward propagation of signals imprinted by variations in the conditions under which this freeze-drying takes place. The upward propagation of these signals produces the well-known stratospheric tape-recorder structure described by . In this study, we define water vapor entry mixing ratios (${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$) as averages over the 20${}^{\circ }$ S–20${}^{\circ }$ N tropical band on the $\mathit{\theta }=\mathrm{400}$ K isentropic surface.

Figure 4

Relationships between tropical tropopause temperatures and water vapor entry mixing ratios (mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ mixing ratios on $\mathit{\theta }=\mathrm{400}$ K averaged over 20${}^{\circ }$ S–20${}^{\circ }$ N). Here tropical tropopause temperatures refer to the minimum in the mean temperature profile within the tropics (20${}^{\circ }$ S–20${}^{\circ }$ N) between $\mathit{\theta }=\mathrm{360}$ and 420 K. The upper set of lines shows the climatological annual cycles for tropical tropopause temperatures, while the lower set of lines shows the climatological annual cycles for water vapor entry mixing ratios. Note that the model results of annual cycles is based on “CLaMS period” (1980–2013), which is longer than “SWOOSH period” starting around 1990.

[Figure omitted. See PDF]

Figure  summarizes the relationship between the climatological annual cycles of tropical tropopause temperatures (upper set of lines) and ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ (lower set of lines). The annual cycles in tropical mean tropopause temperatures based on the three reanalyses (upper set of lines) are similar in both amplitude and phase, but with $\pm \mathrm{1}$ K differences in mean value. Mean ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ differences among the three CLaMS runs are roughly consistent with the corresponding differences in mean tropopause temperatures, with a 1 K difference in temperature corresponding to a $\sim \mathrm{0.6}$ ppmv difference in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios. Although the differences in the mean tropical tropopause temperature are consistent with the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ differences, it is actually the Lagrangian cold point that controls the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values and makes the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ differences, which will be indicated by the next figure. Compared with SWOOSH, simulated annual cycles of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ from CLaMS-ERA and CLaMS-MRA are relatively reasonable in terms of both absolute value and amplitude (i.e., generally within uncertainties). By contrast, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ from CLaMS-JRA is systematically larger than that indicated by SWOOSH or the other two runs. Simulated values based on CLaMS-JRA are $\sim \mathrm{0.5}$ ppmv larger than SWOOSH estimates in February and $\sim \mathrm{1}$ ppmv higher than SWOOSH estimates in September. Thus, CLaMS-JRA produces $\sim \mathrm{0.5}$ ppmv larger peak-to-peak annual amplitude than other estimates. This is qualitatively consistent with nonlinearity in the Clausius–Clapeyron equation, which mean that the effects of temperature offsets on ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ values are larger when the average temperature is higher . Thus, a larger annual amplitude in entry mixing ratios is expected when the average tropopause temperature is larger.

The annual cycle in tropical tropopause temperature leads that of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ by 1–2 months because ascent within the tropical tropopause layer is relatively slow. Given the strong consistency among tropopause temperature annual cycle phases as represented in the reanalyses, we find the modeled phases of entry ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ annual cycle are in a good agreement with only small discrepancies. CLaMS-ERA and CLaMS-JRA indicate very consistent phases of the AC, with the annual minimum occurring in February and the annual maximum occurring in September. CLaMS-MRA shows the annual minimum occurring from February to March and the annual maximum in October. SWOOSH data indicates that the annual minimum occurs in February and the maximum in October. The maximum value of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ based on both of these estimates occurs in October, 1 month later than the maximum values based on CLaMS-ERA and CLaMS-JRA. This phase shift is related to the fact that the ascent in the tropical tropopause layer is slower in MERRA-2 compared to JRA-55 or ERA-I see also comparison of MERRA and ERA-I diabatic heating rates by. Another relevant fact is that the strength of vertical mixing in tropical lowermost stratosphere is important for reproducing the observed seasonality of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$, especially during boreal summer . We did not specify the strength of vertical mixing in this study but intend to examine its impacts in future work.

Figure 5

Climatological mean cold point tropopause temperatures during boreal winter (DJF, a, c, e) and boreal summer (JJA, b, d, f) from MERRA-2 (a, b), along with differences in ERA-I (c, d) and JRA-55 (e, f) relative to MERRA-2. Isolines in the middle and lower rows show absolute values of tropopause temperature. Cold point tropopause temperatures are defined here as minimum temperatures between $\mathit{\theta }=\mathrm{360}$ K and $\mathit{\theta }=\mathrm{420}$ K.

[Figure omitted. See PDF]

In addition to mean tropical tropopause temperature, the spatial distribution of climatological tropopause temperatures is also an important factor in determining the mean entry mixing ratio. Figure  shows this spatial pattern based on MERRA-2 during boreal winter (DJF) and boreal summer (JJA), along with differences in tropopause temperatures from ERA-I and JRA-55 relative to those in MERRA-2. Differences in tropopause temperatures among the three reanalysis products vary considerably with both location and season. ERA-I shows relatively large negative differences ($\le -\mathrm{0.5}$ K) relative to MERRA-2 over the tropical Indian and Pacific oceans during DJF. Differences are particularly notable over the tropical western Pacific, where tropopause temperatures based on ERA-I are substantially colder than those based on MERRA-2 during boreal winter. Tropopause temperatures in this region are widely believed to play a key role in determining the value of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ during boreal winter . Warmer tropopause temperatures in JRA-55 during JJA are located mainly in the Southern Hemisphere subtropics, as well as over the Bay of Bengal and tropical Pacific, with differences relative to MERRA-2 reaching magnitudes of $\sim \mathrm{1}$ K in these regions. The latter two regions are also widely recognized as playing influential roles in troposphere-to-stratosphere transport and the final dehydration of air parcels entering the stratosphere during boreal summer .

Figure 6

Climatological structure of the stratospheric tape-recorder signal based on SWOOSH and the three CLaMS runs. The tape recorder is defined as anomalies in tropical (20${}^{\circ }$ S–20${}^{\circ }$ N) mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ relative to the climatological mean at each level (color shading). The phase of upward propagation (solid circles connected by a line) is defined by the largest correlation with the layer below. For ease of comparison, propagation based on SWOOSH is marked in each panel using gray circles connected by a gray line. Red and blue contours indicate positive and negative contributions of ${\mathrm{CH}}_{\mathrm{4}}$ to ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ anomalies (in units of ppmv, at intervals of 0.02 ppmv).

[Figure omitted. See PDF]

A complementary picture of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values is provided in Fig. , which depicts the tape-recorder signal in SWV averaged over 20${}^{\circ }$ S–20${}^{\circ }$ N. Upward propagation of the tape-recorder signal between 450 and 600 K is 0.5–1.5 months faster in CLaMS-ERA and CLaMS-JRA relative to SWOOSH, and 1–1.5 months slower in CLAMS-MRA than in SWOOSH. Similarly, the amplitude of the tape-recorder signal is systematically stronger than SWOOSH in CLaMS-ERA and CLaMS-JRA, but weaker above 450 K in CLaMS-MRA. These differences are attributable in part to the relatively slow upwelling in MERRA-2 (as shown in Fig. ). Slower upwelling not only delays the propagation of the signal but also allows more time for horizontal advection and mixing of the middle latitude air into the tropics, which tend to damp the signal. We can also see the remarkably strong contribution of ${\mathrm{CH}}_{\mathrm{4}}$ oxidation in CLaMS-MRA, which is shown by the blue and red contour lines in Fig. . The contribution of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ to the tape-recorder signal is substantially larger in CLaMS-MRA than in the other two runs. This feature is a joint effect of slower tropical upwelling and stronger in-mixing from the extratropics, resulting in a relatively pronounced seasonal cycle in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ in CLaMS-MRA with a maximum amplitude of $\sim \mathrm{0.05}$ ppmv near the 450 K isentrope. The amplitude of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ in CLaMS-MRA is twice as large as that in CLaMS-JRA. Meanwhile, CLaMS-ERA shows virtually no anomalies in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ at these levels due to relatively rapid rates of ascent in the lower branch of the BDC. Seasonal variations in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ are opposite in phase relative to seasonal variations in ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$, and account for $\sim \mathrm{20}$ % of the reduced amplitude of the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ tape recorder above 450 K in CLaMS-MRA.

4.2

Global features of the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ annual cycle

We now consider differences in the representation of the annual cycle in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ throughout the global stratosphere. The amplitude of the simulated annual cycle, $A_{\mathrm{AC}}$, and its phase, $P_{\mathrm{AC}}$, are shown in Fig. . Here, we examine the CLaMS simulations in the context of Aura MLS observations during the MLS period (2005–2013). In general, these comparisons show good agreement between the CLaMS simulations and MLS. Results for comparisons against SWOOSH ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ between $\mathit{\theta }=\mathrm{350}$ K and $\mathit{\theta }=\mathrm{650}$ K are qualitatively similar to those against MLS and are not shown here.

Figure 7

AC amplitudes ($A_{\mathrm{AC}}$) and corresponding phases ($P_{\mathrm{AC}}$) during the MLS period ($\sim \mathrm{2005}$–2013, the first and third rows) and the entire CLaMS period (1980–2013, the second and fourth rows). The leftmost panels show $A_{\mathrm{AC}}$ and $P_{\mathrm{AC}}$ based on Aura MLS observations (a1, a3) and the mean of the three reanalysis-driven models (MRM, a2, a4). The corresponding $A_{\mathrm{AC}}$ values from CLaMS-ERA, CLaMS-JRA and CLaMS-MRA are shown in the rightmost three panels of each row. Overlaid “$+$” and “$-$” symbols in panels (b1)–(d1) and (b2)–(d2) demonstrate major differences by more than $\pm \mathrm{30}$ % in simulated $A_{\mathrm{AC}}$ relative to the benchmark estimate shown in the corresponding leftmost panel. Arrows in panels (b3)–(d3) and (b4)–(d4) indicate the major region with phase differences by more than 1 month against the benchmark estimate of $P_{\mathrm{AC}}$ shown in the corresponding leftmost panel. The upward arrows show phase lagging behind of the benchmark phase while the downward arrows show phase ahead of the benchmark phase.

[Figure omitted. See PDF]

The enhanced amplitude of the AC in “region 1” (see Fig. a1 and a2) reflects the tape-recorder signal in the tropical lower stratosphere, which propagates from tropopause to the middle stratosphere and extends into subtropics and midlatitudes. Significant differences in both the amplitude and phase of the AC are evident in the transition layer between the tropical lower stratosphere and the tropical middle stratosphere ($\mathit{\theta }$ equal to 450–600 K). These differences are again linked to discrepancies in the strength of tropical upwelling described above (Fig. a2, b2 and c2 and corresponding discussion). Consistent with our intercomparison of tape-recorder amplitudes shown in Sect. , the amplitude of the AC between $\mathit{\theta }=\mathrm{450}$ K and $\mathit{\theta }=\mathrm{600}$ K is larger in CLaMS-JRA and smaller in CLaMS-MRA relative to that inferred from MLS observations (see the signs of “$+$” and “$-$” overlaid on Fig. b1–d1). Examination of the phase of the AC likewise confirms that AC-related signals propagate faster in CLaMS-ERA and CLaMS-JRA and slower in CLaMS-MRA than indicated by MLS observations (see arrows overlaid on Fig. b3–d3). Although we use different methods to estimate the AC amplitude in this section relative to Sect. , both approaches produce similar results with respect to differences among the simulated AC signals in “region 1”.

MLS observations further indicate a hemispheric asymmetry in the AC amplitude in “region 1” (NH $>$ SH), possibly related to the transport of relatively moist air from the Asian summer monsoon (ASM) anticyclone to the stratosphere during boreal summer e.g.,. This feature is also evident in all three model-based estimates of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$. However, all three CLaMS runs overestimate the amplitude of the AC in the SH subtropics (around 380–450 K). One potential reason is that water vapor in the Southern Hemisphere subtropical lower stratosphere is highly sensitive to small-scale mixing processes that must be parameterized in the model .

Additional local maxima in the AC amplitude are found in the subtropical upper stratosphere of both hemispheres, marked as “region 2” in Fig. . The AC phase in region 2 indicates that this feature propagates quasi-meridionally. Both the enhanced amplitude of the AC and the hemispheric asymmetry (SH $>$ NH) in this part of the stratosphere are related to seasonal variations in the deep branch of the BDC . The amplitude of the feature in region 2 is largest in CLaMS-JRA, followed by CLaMS-ERA and CLaMS-MRA. The amplitude implied by the MLS observations falls between those produced by CLaMS-JRA and CLaMS-ERA and is approximately double the amplitude produced by CLaMS-MRA. Thus, the relative strength of this AC signal in the subtropical upper stratosphere is consistent with the respective strength of diabatic upwelling in the tropical lower and middle stratosphere (see Fig. ), which is strongest in CLaMS-JRA and weakest in CLaMS-MRA. Weaker upwelling in the lower part of the BDC ascending branch in CLaMS-MRA produces strong gradients in tracer concentrations (including both ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and ${\mathrm{CH}}_{\mathrm{4}}$) in the lower stratosphere and weak gradients in the upper stratosphere (see also Fig. ).

5 The quasi-biennial oscillation

We now shift our focus to another quasi-periodic oscillation in SWV: the quasi-biennial oscillation (QBO). Zonal-mean distributions of the climatological amplitudes and phases of the QBO are shown in Fig. . As above, metrics based on MLS observations are provided in the first and third rows and metrics based on the MRM are provided in the second and fourth rows. Two regions show clear peaks in $A_{\mathrm{QBO}}$: the tropical lower stratosphere (region 1) and the tropical upper stratosphere (region 2). A peak in the tropical middle stratosphere (region 3) is evident in CLaMS-ERA and CLaMS-JRA but not in MLS or CLaMS-MRA.

Figure 8

Same as Fig. but for QBO amplitudes and phases. Note that the scale of color bars are different from that in Fig. . Overlaid “$+$” and “$-$” symbols in panels (b1)–(d1) and (b2)–(d2) demonstrate major differences by more than $\pm \mathrm{30}$ % in simulated QBO amplitudes relative to the benchmark estimate shown in the corresponding leftmost panel. The large and small arrows show the major region with phase differences by more than 1 and 3 months, respectively, against the benchmark $P_{\mathrm{QBO}}$ shown in the corresponding leftmost panel. The upward arrows show phase lagging behind of the benchmark phase while the downward arrows show phase ahead of the benchmark phase.

[Figure omitted. See PDF]

The large values of $A_{\mathrm{QBO}}$ in region 1 are due to QBO effects on tropical tropopause temperature. Since the tropopause is colder (warmer) during the easterly (westerly) phase of the QBO (referring here to the 50 hPa wind) , the phase of the QBO effect at the tropical tropopause shown in Fig.  is related to the phase of the 50 hPa QBO wind. This water vapor signal then propagates upward at a speed similar to that of the canonical tape-recorder signal. Further details about the QBO modulation of the tape recorder were suggested by . Similar to what we found for AC, CLaMS-ERA and CLaMS-JRA both overestimate $A_{\mathrm{QBO}}$ and show upward phase propagation $\sim \mathrm{1}$ month faster than MLS at isentropic levels between 450 and 550 K. The QBO amplitude in region 1 based on CLaMS-MRA compares better with $A_{\mathrm{QBO}}$ based on MLS observations. These amplitude differences among CLaMS runs evokes similar differences in $A_{\mathrm{AC}}$ in the tropical lower stratosphere. Both biases can be traced back to strong upward transport in the lower to middle stratosphere in ERA-I and JRA-55 (Fig. ), which results in weaker damping in the amplitudes of tropopause-originating signals (such as $A_{\mathrm{AC}}$ in Fig.  region 1 and $A_{\mathrm{QBO}}$ in Fig.  region 1) in the vertical direction by mixing.

The large values of $A_{\mathrm{QBO}}$ in region 2 are mainly linked to QBO-related modulation of the stratospheric circulation and references therein. The corresponding phase is effectively simultaneous with the 50 hPa QBO wind, with very little spatial shift. Satellite-derived estimates of this QBO signal indicate that amplitudes are in the range of 0.2 to 0.6 ppmv . MLS ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ is on the lower end of this range with an amplitude of $\sim \mathrm{0.2}$ ppmv. However, CLaMS runs during the MLS period produce even smaller amplitudes of 0.1 to 0.2 ppmv, below the range implied by satellite measurements. The smallest amplitude is produced by CLaMS-ERA. This discrepancy may be related to lower values of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and $\mathit{\alpha }$ in the upper stratosphere as shown in Fig. a3, which result from a systematically faster stratospheric circulation in ERA-I.

The amplitude and phase of the QBO signal in SWV show pronounced discrepancies in the middle stratosphere among the CLaMS simulations driven by different reanalyses (region 3 in Fig. ). This region coincides with the strongest QBO signal in zonal wind. During the MLS period (2005–2013), both CLaMS-ERA and CLaMS-JRA overestimate $A_{\mathrm{QBO}}$ in this region, particularly within the 0–20${}^{\circ }$ N latitude band. These overestimates can be attributed to strong upward transport of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{e}}$ from the lower stratosphere. Although CLaMS-MRA underestimates $A_{\mathrm{QBO}}$ in this region, the spatial pattern produced by CLaMS-MRA in the middle stratosphere is more consistent with that indicated by MLS. The estimate of phase in this region from CLaMS-MRA shows $\sim \mathrm{6}$ months difference when years before MLS period are included or excluded (see Fig. d3–d4, where phase shifts are indicated by arrows). Uncertainties in QBO-related ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ anomalies arise not only from different circulation responses but also from disagreements in zonal wind at the Equator among the reanalyses. Equatorial winds at 10 hPa ($\sim \mathrm{800}$–900 K) in MERRA-2 are clearly different from the FUB records, ERA-I and JRA-55 during the 1980s and early 1990s, when the assimilated observations are less able to pull the forecast model away from its own internally generated QBO signal . Considering the consistency in QBO wind between MERRA-2 during MLS period and the FUB records, the QBO phase for CLaMS-MRA (d3) is more reliable than that during the whole CLaMS period (d4).

6 Trends and residual sources of natural variability

As mentioned in Sect.  and observed in the SWOOSH record (see Fig. e), residual variability beyond the AC, QBO and SAH (${\mathit{\chi }}_{\mathrm{res}}\left(t\right)$ in Eq. ) also contributes substantially to variations in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios. Figure a shows that the three CLaMS simulations produce similar interannual variability in the ${\mathit{\chi }}_{\mathrm{res}}\left(t\right)$ term. Within the SWOOSH period (1993–2013), representations of interannual variability among the simulations and the SWOOSH analysis are more consistent during the Aura MLS period (2004–2013) than the pre-MLS period (1993–2004), likely due to significant improvement in SWOOSH data once MLS observations were incorporated as well as the contemporaneous improvements in the quantity and quality of satellite observations assimilated by the reanalysis systems . However, significant differences are evident in the simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratio trends, as discussed later in this section.

Figure 9

Comparison of residual variability after removing the AC, QBO and SAH from ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values (mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ mixing ratios on $\mathit{\theta }=\mathrm{400}$ K averaged over 20${}^{\circ }$ S–20${}^{\circ }$ N) based on CLaMS simulations and SWOOSH. Panel (a) shows the residual term (${\mathit{\chi }}_{\mathrm{res}}$) after applying a 3-month running mean. The long-term linear fit is plotted as a solid line with $\pm \mathrm{2}\mathit{\sigma }$ uncertainties. Panel (b) shows additional components of the multiple linear regression, including the influences of variations in volcanic aerosol (upper set of lines) and ENSO (lower set of lines). Panel (c) shows the residuum of (a) after subtracting all contributors shown in (b), along with the corresponding trend.

[Figure omitted. See PDF]

have discussed the causes of interannual variability in SWV using a regression model, which considered variability in the QBO, BDC and tropical tropospheric temperature. Although this regression model is successfully captures most of the variability, it is difficult to disentangle the effects from each regressor. We intend to explore the processes that ultimately control the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ variability, i.e., the processes that might influence SWV by modulating BDC or tropical tropospheric temperatures. Among the potential controlling processes, ENSO and stratospheric volcanic aerosols have been shown to modulate both the tropical ascending branch of the BDC e.g., and tropical tropopause temperatures e.g.,. Consequently, variations in ENSO and stratospheric volcanic aerosols have significant influences on ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios e.g.,. To clarify the ENSO and volcanic aerosol impacts on the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios in each simulation, we analyze the residuum of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values (${\mathit{\chi }}_{\mathrm{res}}\left(t\right)$ in Eq. ) using the following multiple regression model: $$\begin{array}{cc}{\displaystyle }{\mathit{\chi }}_{\mathrm{res}}\left(t\right)=& {\displaystyle }a\cdot P_{\mathrm{volcano}}(t-{\mathit{\tau }}_{\mathrm{volcano}})+b\cdot P_{\mathrm{ENSO}}(t-{\mathit{\tau }}_{\mathrm{enso}})\\ \text{6}& {\displaystyle }& {\displaystyle }+{\mathit{\chi }}_{\mathrm{res}}^{\prime }\left(t\right),\end{array}$$ where $P_{\mathrm{ENSO}}$ is the normalized Multivariate ENSO Index (MEI; ) and $P_{\mathrm{volcano}}$ is the aerosol optical depth (AOD) as recorded in satellite data . The coefficients are the amplitude $a$ and lag ${\mathit{\tau }}_{\mathrm{volcano}}$ associated with volcanic aerosols and the amplitude $b$ and lag ${\mathit{\tau }}_{\mathrm{enso}}$ associated with ENSO. We determine the parameters $a$, $b$, ${\mathit{\tau }}_{\mathrm{volcano}}$ and ${\mathit{\tau }}_{\mathrm{enso}}$ as the parameter set that minimizes the residual (${\mathit{\chi }}_{\mathrm{res}}^{\prime }$) in the least-squares sense. Additional details of the method and its application have been summarized by .

Figure b shows the response of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ to stratospheric volcanic aerosol (upper set in panel b) and ENSO (lower set in panel b) based on the simulations and the observations. Beyond some differences in magnitude, the volcanic aerosol-induced changes in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values agree well in terms of response signs and lags. The increase in modeled ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ associated with volcanic aerosols arises mainly from aerosol-induced warming of the TTL . Positive ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry anomalies induced by the two major volcanic eruptions during the CLaMS period are twice as large in CLaMS-JRA and CLaMS-MRA (0.6 ppmv for El Chichón and 0.8 ppmv for Pinatubo) as in CLaMS-ERA (0.3 ppmv for El Chichón and 0.4 ppmv for Pinatubo). Note that among these three reanalyses, only MERRA-2 explicitly includes perturbations to the stratospheric aerosol burden following volcanic eruptions , although the effects of these perturbations may also enter all of the reanalyses indirectly through the assimilation of temperature, ozone and other affected quantities.

As with the volcanic aerosol-induced effects, ENSO impacts on ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values from CLaMS and SWOOSH agree well with respect to the main characteristics of the response. All four estimates show positive ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry anomalies during El Niño and negative anomalies during La Niña, consistent with previous studies e.g.,. Lags (${\mathit{\tau }}_{\mathrm{enso}}$) in ENSO-induced tropical mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ anomalies at $\mathit{\theta }=\mathrm{400}$ K relative to the ENSO signal (here represented by the MEI index) are likewise consistently between 10–11 months among all three simulations and the observations. However, there are notable quantitative differences in the magnitude of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratio changes induced by ENSO variability. CLaMS-JRA and CLaMS-MRA show increases in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios about $\sim \mathrm{0.1}$ ppmv larger than CLaMS-ERA during strong El-Niño events (1982/1983, 1986/1987 and 1997/1998). However, the response in the latter exhibits the best agreement with SWOOSH during the 1997/1998 El Niño.

Figure c shows the residual terms (${\mathit{\chi }}_{\mathrm{res}}^{\prime }$) and the trends related to the residual terms after removing ENSO and volcanic aerosol-induced variability. We notice that the significant differences among the trend estimates shown in Fig. a largely remain in the residual variability. Further quantification of contributions to the trends is listed in Table . We also notice a substantial quasi-decadal variability remaining in the residuals. After removing the high-frequency variability and the linear trends, this quasi-decadal variability agrees well among the three simulations as well as with the SWOOSH data. Further interpretation of this consistent quasi-decadal signal is out of the scope of this study, which aims mainly to intercompare different reanalyses.

Table 1

Trend estimates in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios and contributions from different sources of variability based on the three CLaMS simulations and SWOOSH over the SWOOSH period (1993–2013) and the longer CLaMS period (1980–2013). Trends are reported in units of parts per million by volume per decade (left of the slash) and percent per decade (right of the slash). Statistically significant trends are marked in bold.

	SWOOSH	CLaMS-ERA	CLaMS-JRA	CLaMS-MRA
$C_{\mathrm{trd}}$
1993–2013	–0.22/–5.8	–0.06/–1.6	–0.39/–8.6	–0.06/–1.5
1980–2013		0.01/0.2	–0.22/–4.9	–0.08/–2.0
$\mathrm{\Delta }{C}_{\mathrm{trd}}$ (ENSO)
1993–2013	–0.05/–1.3	–0.06/–1.5	–0.07/–1.6	–0.05/–1.3
1980–2013		–0.02/–0.6	–0.03/–0.6	–0.02/–0.5
$\mathrm{\Delta }{C}_{\mathrm{trd}}$ (volcano)
1993–2013	0.00/0.0	0.02/0.4	0.03/0.7	0.03/0.6
1980–2013		0.00/0.0	0.00/0.0	0.00/0.0
$\mathrm{\Delta }{C}_{\mathrm{trd}}$ (res.)
1993–2013	–0.18/–4.6	–0.02/–0.6	–0.35/–7.8	–0.03/–0.8
1980–2013		0.03/0.7	–0.19/–4.3	–0.06/–1.5

Table  lists the trends in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios along with the estimated contributions of ENSO and volcanic aerosols to trends in the three reanalysis-driven simulations and SWOOSH over the shorter period (1993–2013) and the longer-term period (1980–2013). Trends are calculated by analyzing the residual of the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry anomalies with and without the ENSO and volcanic aerosol signals. The trend produced by CLaMS-ERA is very consistent with that produced by CLaMS-MRA after 1985 (not shown), as is the residual variability shown in Fig. c. By contrast, CLaMS-JRA produces a strong negative trend. As nonperiodic variability, such as ENSO variability and the volcanic aerosol burden, may contribute to these trends, their estimated contributions are listed in Table . ENSO variability contributes to negative trends in both modeled and observed ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios over the 1993–2013 SWOOSH period. For the entire simulation period (1980–2013), the ENSO contribution to trends is also negative but with a reduced magnitude. The variability due to volcanic eruptions has a zero-trend contribution to the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values over the longer period (1980–2013), which is consistent with the volcanic-aerosol-induced zero effect on the trends of mean age of air at the tropical lower stratosphere shown in . However, it has a positive trend contribution during the 1993–2013 SWOOSH period, likely due to the increased frequency of minor volcanic eruptions after 2008 . Conversely, the volcanic aerosol contribution to the trend from 1980 to 2008 (not shown) is negative due to the forcing of the two major eruptions. This analysis suggests that the contribution of volcanic aerosols to interannual variability in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry values is highly period dependent. Although we find near compensation over the CLaMS period used in this work (1980–2013), trends in SWV can be substantially influenced by periods of enhanced or reduced supplies of volcanic aerosol to the stratosphere.

Differences among the trend estimates persist in the residual variability even after accounting for the effects of ENSO and volcanic aerosol (Fig. c and last row in Table ), indicating that uncertainties in the model-based trends do not emerge from different responses to major volcanic eruptions or strong ENSO events. The model-based ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ trends are sensitive to the representation of the tropical tropopause temperature from reanalyses, especially for the early years of the satellite era (e.g., 1980–1985) . Differences in the observations assimilated by the three reanalyses or in the representation of the BDC may also contribute to differences in the trend estimates. Reliable attribution of differences in simulated SWV trends will require disentangling complex interactions among changes in tropopause temperature, the stratospheric circulation and anthropogenic factors such as the amount of methane entering the stratosphere at global scale.

7 Discussion

In the sections above, we relate differences in representations of the climatological mean, AC and QBO signals of SWV to differences in upwelling rates in the shallow branch of the BDC among the reanalysis products used to drive the model. One remaining question is whether these differences in upwelling strength are also responsible for differences in AC or QBO variance. Analysis of the total standard deviation and its contributors in each simulation can shed light on this question. If different upwelling strengths are the main factor, then these different rates of vertical advection should influence each contributor similarly to the total variance. In other words, a larger (or smaller) magnitude of the AC or QBO variance should be approximately in proportion to a larger (or smaller) magnitude of the total variance. We use the standard deviation from the full ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ time series (${\mathit{\sigma }}_{\mathit{\chi }}$, where $\mathit{\chi }$ is the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ volume mixing ratio) to quantify the magnitude of total variability. This metric is shown for each of the three simulations in Fig. a–c. We then check the fraction of the total variance attributable to the AC and QBO, respectively, ${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{AC}}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathit{\chi }}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{QBO}}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathit{\chi }}^{\mathrm{2}}$. The variance fractions attributable to the AC and QBO components of the total signal are shown in Fig. d–f as gray shading and purple contours, respectively.

Figure 10

The standard deviation in monthly-mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ (${\mathit{\sigma }}_{\mathit{\chi }}$; upper row) and the fraction of AC and QBO variance relative to total variance (in %; d, e, f) for CLaMS-ERA (a, d), CLaMS-JRA (b, e) and CLaMS-MRA (c, f). Shading in the lower panels shows the variance fraction attributed to the annual cycle (${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{AC}}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathit{\chi }}^{\mathrm{2}}$), while purple contours show the variance fraction attributed to the QBO (${\mathit{\sigma }}_{{\mathit{\chi }}_{\mathrm{QBO}}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathit{\chi }}^{\mathrm{2}}$).

[Figure omitted. See PDF]

Figure  indicates that zonal-mean distributions of AC and QBO variance fractions are qualitatively consistent across all three simulations (Fig. d–f), although the magnitudes of total variance differ to some extent among the three CLaMS runs (Fig. a–c). For example, CLaMS-MRA produces the smallest amplitudes of both the AC and QBO signals in the lower and middle stratosphere among the three simulations. Figure. c (corresponding to CLaMS-MRA) shows that these smaller amplitudes are in turn associated with weaker total variance in the middle stratosphere. By contrast, the fractional contributions attributed to the AC and QBO (Fig. f) are comparable to if not larger than the corresponding fractional contributions based on CLaMS-ERA or CLaMS-JRA. This similarity among variance fractions implies that differences in the strength of the tropical upwelling introduced via the prescribed reanalysis diabatic heating rates can adequately explain differences in the amplitudes of the AC and QBO signals produced by the CLaMS simulations.

Figure 11

Joint distributions of daily-mean gridded diabatic heating rates ($\dot{\mathit{\theta }}$) in the tropical lower stratosphere ($\mathit{\theta }=\mathrm{380}$ K to $\mathit{\theta }=\mathrm{460}$ K) against the corresponding long-wave cloud radiative effect (LWCRE) based on ERA-I (a), JRA-55 (b) and MERRA-2 (c) during 1980–2013. The LWCRE is calculated as clear-sky minus all-sky outgoing long-wave radiation. Also shown are means and standard deviations of ${\dot{\mathit{\theta }}}_{\mathrm{LS}}$ composited into LWCRE bins in intervals of 25 W m${}^{-\mathrm{2}}$ starting from zero. Estimates based on ERA-I (blue), JRA-55 (purple) and MERRA-2 (red) are shown on all three panels to facilitate comparison. Distributions are based on variables archived on $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ (ERA-I) and $\mathrm{1.25}{}^{\circ }\times \mathrm{1.25}{}^{\circ }$ (JRA-55 and MERRA-2) latitude–longitude grids between 30${}^{\circ }$ S and 30${}^{\circ }$ N.

[Figure omitted. See PDF]

The reason that diabatic heating rates in the tropical lower stratosphere are smaller in MERRA-2 than in ERA-I or JRA-55 is as yet unclear. Although a full attribution is beyond the scope of this study, we consider here two possible contributors that may be informative for other users of these data products: (1) the larger long-wave cloud radiative effect (LWCRE) in MERRA-2 than ERA-I or JRA-55 and (2) the unique assimilation process used in MERRA-2.

Diabatic heating rates immediately above the tropical tropopause are dominated by radiative heating $Q$ e.g.,, which, under the plane-parallel assumption applied in the reanalysis models, indicates the net vertical convergence of energy in the form of radiation. Adopting the Newtonian cooling approximation $Q\sim -\mathit{\alpha }(T-T_{\mathrm{eq}})$ (with $\mathit{\alpha }$ the inverse of the radiative relaxation time), radiative heating rates $Q$ can be treated as inversely proportional to the difference between the local temperature $T$ and a local radiative equilibrium temperature $T_{\mathrm{eq}}$ e.g.,. The latter depends on the vertical profile of temperature, the chemical composition, and the radiative effects of aerosol and clouds within the atmospheric column. Figure  illustrates the relationship between daily-mean gridded diabatic heating rates in the tropical lower stratosphere ${\dot{\mathit{\theta }}}_{\mathrm{LS}}$ (30${}^{\circ }$ S–30${}^{\circ }$ N; $\mathit{\theta }=\mathrm{380}$ K to $\mathit{\theta }=\mathrm{460}$ K) and the corresponding LWCRE (defined as clear-sky minus all-sky outgoing long-wave radiation at the nominal top of atmosphere) based on ERA-I, JRA-55 and MERRA-2. Daily-mean diabatic heating rates and associated upwelling in the tropical lower stratosphere are negatively correlated with LWCRE in all three reanalyses, with $r$ ranging from $-\mathrm{0.18}$ in JRA-55 to $-\mathrm{0.38}$ in ERA-I and $-\mathrm{0.54}$ in MERRA-2. This negative relationship may be explained by noting that clouds occur almost exclusively within the troposphere, so that a larger LWCRE corresponds to a smaller upward flux of long-wave radiation across the tropopause. All else remaining equal, a reduced upward flux across the tropopause acts to reduce $T_{\mathrm{eq}}$ and the net convergence of radiant energy, and therefore implies a smaller diabatic heating rate. Moreover, values of the LWCRE in MERRA-2 (median value: 15.9 W m${}^{-\mathrm{2}}$) are systematically larger than those in ERA-I (10.0 W m${}^{-\mathrm{2}}$) or JRA-55 (6.8 W m${}^{-\mathrm{2}}$), especially at the upper end of the range (cf. 90th percentile values of 82.9, 50.1, and 34.8 W m${}^{-\mathrm{2}}$, respectively). It is thus unsurprising that MERRA-2 produces weaker diabatic upwelling near the base of the tropical pipe. However, this effect does not appear to account for the entire difference as lower stratospheric heating rates composited for the same ranges of LWCRE are still systematically smaller in MERRA-2 than in ERA-Interim or JRA-55 (Fig. ).

Another potential contributor is the unique assimilation process used in MERRA-2 (and its predecessor MERRA), which includes an additional “corrector” step after the initial “3D-Var” analysis i.e., an “iterative predictor-corrector approach”; see also. During the corrector step, all analysis increments are applied gradually over time as additional tendency terms. This procedure improves the internal consistency among analyzed variables and other model variables. However, it also means that the diabatic tendency and other physical tendency terms in MERRA-2 are archived during the corrector step, whereas other reanalyses produce these terms during an initial forecast analogous to the predictor step in MERRA-2. This may have the unintended effect of systematically damping (or amplifying) the tendency terms produced by the physical parameterizations. Again adopting the Newtonian cooling approximation, systematic biases in $T$ or radiatively active constituents could result in the analysis tendency constantly acting to reduce $T-T_{\mathrm{eq}}$ in the lower stratosphere, thereby damping radiative heating rates and associated diabatic upwelling. Using “replay” simulations, which mimic the corrector forecast using a modified version of the atmospheric model and data assimilation system used for MERRA-2, showed that a simulation constrained by time-averaged assimilated fields produced approximately 30 % slower ascent in the tropical lower stratosphere than one based on instantaneous analysis fields. Although concluded that the slower ascent in the simulation based on time-averaged assimilated fields produced output in better agreement with available observations, their result nonetheless highlights the potential impact of the predictor–corrector approach on upwelling in the tropical lower stratosphere. Additional evidence that data assimilation may suppress lower stratospheric upwelling in MERRA-2 is provided by comparison against the MERRA-2 AMIP dataset , a 10-member ensemble of free-running simulations generated using the same model and boundary conditions as MERRA-2 but without data assimilation. Time-mean zonal-mean diabatic heating rates in the tropical lower stratosphere are about 0.1–0.2 K d${}^{-\mathrm{1}}$ smaller in MERRA-2 than in MERRA-2 AMIP over the period 1980–2017 (see Fig.  in Appendix ).

8 Conclusions

We evaluate representations of SWV and its variations at multiple timescales using the Chemical Lagrangian Model of the Stratosphere (CLaMS) driven by horizontal winds and diabatic heating rates from three recent atmospheric reanalyses: ERA-I, JRA-55 and MERRA-2. The analysis is based on CLaMS simulations of monthly-mean zonal-mean ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ from 1980–2013. We present an intercomparison of simulated variations in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the context of observational estimates from Aura MLS and the Stratospheric Water and Ozone Satellite Homogenized (SWOOSH) database, focusing on the annual cycle (AC), the quasi-biennial oscillation (QBO), and long-term variability and trends. Based on the results of this intercomparison, we reach the following conclusions.

The climatological means of SWV, which represents a combination of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entering the stratosphere through the tropical tropopause (dominant in the lower stratosphere) and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ supplied by ${\mathrm{CH}}_{\mathrm{4}}$ oxidation (dominant in the upper stratosphere), are in a good agreement (within $\pm \mathrm{10}$ % differences) among the three simulations. The “age-of-air” (AoA) and the ratio of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ from ${\mathrm{CH}}_{\mathrm{4}}$ oxidation in the lower stratosphere are larger in CLaMS-MRA than in the other two runs, consistent with relatively weaker diabatic upwelling in the lower and middle tropical stratosphere in MERRA-2 relative to ERA-I or JRA-55. Mean tropopause temperatures in ERA-I are approximately 1 K colder than those in JRA-55, resulting in a tropical lower stratosphere that is $\sim \mathrm{0.6}$ ppmv drier in CLaMS-ERA than in CLaMS-JRA. CLaMS-MRA produces a moderate ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratio close to observed as a result of moderate mean tropopause temperature among the three reanalyses.

Differences in amplitudes of modeled ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratio ACs can largely be explained by differences in mean tropical tropopause temperatures. Tuning of the dehydration scheme could possibly eliminate differences in both mean entry values and the AC amplitude, but not the phase see also. The spatial distribution of tropopause temperatures is also important to understanding differences in seasonal ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios. The relative dryness of CLaMS-ERA is most pronounced in boreal winter, consistent with especially large negative biases in tropopause temperatures over the western Pacific in ERA-I relative to MERRA-2. By contrast, CLaMS-JRA produces particularly large values of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios during boreal summer, consistent with systematically warmer tropopause temperatures over the Bay of Bengal and tropical Pacific.

The AC and QBO signals in simulated SWV are robust across the three simulations and reasonably consistent with observations. The main discrepancies in the AC and QBO components of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ variability are located in the tropical lower and middle stratosphere. The observed AC and QBO signals are generally located between those produced by CLaMS-ERA and CLaMS-JRA (too strong) and those produced by CLaMS-MRA (too weak). That is mainly linked to the fact that the upwelling rates in the tropical lower stratosphere are smaller in MERRA-2 than ERA-I or JRA-55, which is attributable to the larger long-wave cloud radiative effect and the unique assimilation process in MERRA-2. CLaMS-MRA produces a more realistic pattern of QBO-related signals in SWV (but a bit weaker than observation) than the other two runs, especially in the middle stratosphere. The “tape recorder” signal is 25 % weaker and $\sim \mathrm{1.5}$ month slower (between 450 and 600 K) in CLaMS-MRA relative to SWOOSH. By contrast, the tape-recorder signals are both stronger and faster than observed in CLaMS-ERA (20 % stronger and $\sim \mathrm{1}$ month faster than SWOOSH) and CLaMS-JRA (40 % stronger and $\sim \mathrm{1}$ month faster than SWOOSH). We find that differences in the rate of tropical upwelling among the reanalyses not only modulate the propagation of the simulated tape-recorder signal, but also interact with ${\mathrm{CH}}_{\mathrm{4}}$ photochemistry in ways that tend to amplify the differences due to propagation rates alone (Sect. ).

With respect to residual variability in SWV, the consistency among the model results and the SWOOSH analysis improves with time, including the responses of SWV to variations in ENSO and volcanic aerosol. All simulations indicate strong and consistent quasi-decadal variability. Trends in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ entry mixing ratios over the 1980–2013 period are highly sensitive to the reanalysis used to drive the model, leading to different magnitudes and even different signs. This sensitivity of simulated ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ trends to choice of reanalysis cannot be fully explained by discrepancies in the response of modeled ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ to ENSO variability or major volcanic eruptions. The ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ trends are particularly sensitive to the quality of the reanalyses during the first 5 years of the analysis period (1980–1985).

As best estimates of the true state of the atmosphere, meteorological reanalyses are widely used to drive tracer transport or constrain the meteorological state during model simulations. Our study indicates that the ERA-I, JRA-55 and MERRA-2 reanalyses are all capable of reproducing the seasonality and interannual variability in SWV related to QBO, ENSO and volcanic aerosol. However, particular attention should be paid to potential reanalysis-dependent biases when studying longer-term variability in SWV or other atmospheric tracers. Moreover, the available observations allow no conclusions about which reanalysis is most realistic. Therefore, multiple reanalyses should be used whenever possible to better characterize the effects of uncertainties in the atmospheric state.

Data availability

The CLaMS model data used for this paper may be requested from the corresponding author ([email protected])

Appendix A

SWOOSH period of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ data used for comparison with CLaMS

In Sect. , we have clarified the procedure used to determine the “SWOOSH period”. The start time for each latitude–$\mathit{\theta }$ grid location is the first month when SWOOSH ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ has more than 12 months of available data within the next 2 years. The end time is always December 2013. This procedure excludes some early months of the SWOOSH record (mainly intermittent observations from SAGE-II in the lower stratosphere), thereby improving the continuity of the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ time series. Figure  shows the distribution of the SWOOSH period start years identified via this procedure. Due to data filtering applied during the production of SWOOSH , the SWOOSH period is typically (1) shorter at lower levels than at higher levels, (2) shorter in the tropics and polar regions than in the midlatitudes and (3) shorter in the Southern Hemisphere than in the Northern Hemisphere. The shortest SWOOSH periods start from 2004 and are found in the tropical lower stratosphere and at Southern Hemisphere high latitudes. The longest SWOOSH period starts from 1986, in the midlatitudes between $\mathit{\theta }=\mathrm{550}$ K and $\mathit{\theta }=\mathrm{650}$ K.

Although SWOOSH ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ data have some gaps from the start month (shown in Fig. ) to the end of 2013, the full monthly and zonally mean modeled ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ are taken into account for comparison. This allows us to keep the modeled ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ time series as continuous as possible, which in turn allows us to extract more reliable variance estimates. The results are virtually unchanged when (SWOOSH) missing months are also excluded from the CLaMS runs.

Figure A1

Distribution of start years for the “SWOOSH period” used to compare against the model results. The start time for each latitude–$\mathit{\theta }$ grid point is the first month for which more than 12 months of SWOOSH ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ data are available over the following 2 years.

[Figure omitted. See PDF]

Appendix B Tropospheric methane in CLaMS

The boundary condition for ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratios in CLaMS is prescribed in the lowest model layer, corresponding to the hybrid vertical coordinate $\mathit{\zeta }$ equal to 0–100 K. For the simulation period between 1985 and 2012, ${\mathrm{CH}}_{\mathrm{4}}$ measurements from the NOAA/CMDL ground-based measurement network are used . For simulations in 2012 and later, observations from the AIRS instrument are used. The ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratios between 1980 and 1984 in Fig. B1 are just repetition of the year 1985 since the NOAA/CMDL data are not available before 1985. As ${\mathrm{CH}}_{\mathrm{4}}$ is a long-lived tracer, we use zonal-mean mixing ratios to prescribe the lower boundary condition. Figure  shows the time series of global-mean ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratios prescribed at the model lower boundary. Since the time series before 1985 is artificial, Fig.  shows a start of the methane from the mid-1980s before leveling off in the mid-1990s. The long-term trend in ${\mathrm{CH}}_{\mathrm{4}}$ at the model lower boundary (1980–2013) is 0.054 ppmv decade${}^{-\mathrm{1}}$. In fact, there is one option to extend the period of tropospheric methane before 1985 following by using the Antarctic ice core data . We will consider this option for the lower boundary condition in future applications of CLaMS.

Figure B1

Temporal variations from 1980–2013 in the global mean methane mixing ratio prescribed at the model lower boundary, i.e., the hybrid vertical coordinate $\mathit{\zeta }$ equal to 0–100 K.

[Figure omitted. See PDF]

Appendix C Semiannual variation

Figure  shows the spatial distribution of the amplitude and corresponding variance fraction of semiannual harmonic (SAH) signals in simulated SWV. The enhanced SAH signal in region 1 is associated with the semiannual oscillation (SAO) of the upper stratospheric circulation. Satellite observations of the SAH feature in the tropical upper stratosphere and the mechanisms behind this feature have been studied by and . The enhanced SAH signals in region 2 and region 3 are located between two strong AC signals in the polar regions (SH $>$ NH; see Fig. ). The mechanism behind the SAH signals in regions 2 and 3 is related to the combined effects of seasonality in the formation of polar stratospheric clouds (PCSs) and seasonality in vertical transport e.g.,.

Figure C1

The top row shows SAH amplitudes in the same way as in the second row of Fig. . The bottom row shows the fraction of SAH variance relative to total variance, similar to the bottom row of Fig. .

[Figure omitted. See PDF]

The SAH fraction shows that the SAH contributed more than 10 % of the total variance in region 1 and region 3 (both located in the upper stratosphere). The SAH amplitudes in region 1 simulated by the CLaMS runs are in good agreement (within 0.1–0.2 ppmv), and also agree well with HALOE and MLS observations . As with the AC and QBO features, the SAH signal in region 1 is most pronounced in CLaMS-JRA, followed by CLaMS-ERA and CLaMS-MRA. Similarly, the SAH fraction relative to the total variance is again higher in CLaMS-MRA (greater than 15 %–20 %) due to the low total variance simulated by CLaMS-MRA in the mid-upper stratosphere.

Appendix D MERRA-2 versus MERRA-2 AMIP diabatic heating rates

The MERRA-2 AMIP dataset is an ensemble of 10 Atmospheric Model Intercomparison Project (AMIP)-type simulations conducted using the atmospheric model used to produce MERRA-2. All simulations used identical boundary conditions and model configurations to MERRA-2 but did not assimilate observations . Therefore, a comparison of MERRA-2 against that of MERRA-2 AMIP gives a rough indication of the effect of the iterative predictor–corrector data assimilation procedure used in MERRA-2. The comparison of the diabatic heating rates is shown in Fig. . The climatological patterns in MERRA-2 and MERRA-2 AMIP are qualitatively similar; however, their quantitative differences are substantial (bottom panel). Most notably, diabatic heating rates are about 0.2–0.4 K d${}^{-\mathrm{1}}$ larger in MERRA-2 than in MERRA-2 AMIP at pressures greater than 200 hPa, while diabatic heating rates in the lower stratosphere (pressures less than 100 hPa) are 0.1–0.2 K d${}^{-\mathrm{1}}$ smaller in MERRA-2 than in MERRA-2 AMIP. The results suggest that the data assimilation procedure suppresses upwelling in the tropical lower stratosphere by as much as 20 %–40 %.

Figure D1

Comparison of time-mean zonal-mean diabatic heating rates in the tropical (30${}^{\circ }$ N–30${}^{\circ }$ S) lower stratosphere between the MERRA-2 reanalysis (a) and the MERRA-2 AMIP ensemble mean (b) over 1980–2017. Differences between MERRA-2 and MERRA-2 AMIP are shown in (c).

[Figure omitted. See PDF]

Author contributions

MT carried out the analysis on the simulations and on the reanalysis data. PK and FP contributed the code for analysis. JSW prepared the MERRA-2 reanalysis data and provided the attribution for slow upwelling rates in MERRA-2 for the discussion section. PK, FP and XY performed the simulations driven by ERA-I, JRA-55 and MERRA-2, respectively. MD contributed the regression on natural variability. SF and MR gave helpful suggestions, especially on the dehydration scheme in the model and the interpretation of trends. MT wrote the paper. All the co-authors provided helpful discussions and comments on the paper.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “The SPARC Reanalysis Intercomparison Project (S-RIP) (ACP/ESSD inter-journal SI)”. It is not associated with a conference.

Acknowledgements

This research was supported by a joint DFG–NSFC research project with DFG project number 392169209 and NSFC project number 20171352419. Mengchu Tao thanks the German Helmholtz-Gemeinschaft within the Helmholtz-CAS Joint Research Group (JRG) “Climatological impact of increasing anthropogenic emissions over Asia”, enabling her research position in Institute of Energy and Climate Research, Stratosphere (IEK-7), Forschungszentrum in Jülich during which this work was carried out. We thank the Helmholtz Association under grant number VH-NG-1128 (Helmholtz-Hochschul-Nachwuchsforschergruppe), providing the research funding for the young investigator group in IEK-7 including the co-authors Felix Ploeger and Mohamadou Diallo. We thank Krzysztof Wargan for assistance with the MERRA-2 AMIP data. The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association.

Review statement

This paper was edited by Farahnaz Khosrawi and reviewed by three anonymous referees.