Interannual-to-decadal variability of the stratosphere during the 20th century: ensemble simulations with a chemistry-climate model

. Interannual-to-decadal variability in stratospheric ozone and climate have a number of common sources, such as variations in solar irradiance, stratospheric aerosol load-ing due to volcanic eruptions, El Ni˜no Southern Oscillation variability and the quasi-biennial oscillation (QBO). Cur-rently available data records as well as model simulations addressing stratospheric chemical climate variability mostly cover only the past few decades, which is often insufﬁcient to address natural interannual-to-decadal variability. Here we make use of recently reconstructed and re-evaluated data products to force and validate transient ensemble model simulations (nine members) across the twentieth century com-puted by means of the chemistry-climate model SOCOL (SOlar Climate Ozone Links). The forcings include sea surface temperatures, sea ice, solar irradiance, stratospheric aerosols, QBO, changes in land properties, greenhouse gases, ozone depleting substances, and emissions of carbon monox-ides, and nitrogen oxides. The transient simulations are in good agreement with observations, reconstructions and reanalyses and allow quantiﬁcation of interannual-to-decadal variability during the 20th century. All ensemble members are able to capture the low-frequency variability in tropical and mid-latitude total ozone as well as in the strength of the subtropical jet, suggesting a realistic response to external forcings in this area. The region of the northern polar vortex exhibits a large internal variability that is found in the frequency, seasonality, and strength of major warmings as well as in the strength of the modeled polar vortex. Results from process-oriented analysis, such as correlation between the vertical Eliassen Palm ﬂux (EP ﬂux) component and polar as ozone in order


Introduction
The stratosphere exhibits chemical and dynamical variability on different time scales, ranging from day-to-day variability to interdecadal and longer variability, as well as exhibiting secular trends (Solomon, 1999). Through stratosphere-troposphere coupling, interannual-to-decadal variability originating from the troposphere can affect stratospheric climate (e.g., Hadjinicolaou et al., 2002) and vice versa (see e.g., Baldwin et al., 2003). For instance, volcanic eruptions (e.g., Robock, 2000;Rozanov et al., 2002), solar variability (e.g., Hood, 1999;Egorova et al., 2004), ozone depletion (Gillett and Thompson, 2003), or El Niño Southern Oscillation (ENSO) (e.g., Sassi et al., 2004) may all affect both the troposphere and the stratosphere, and it is an open question to what extent climate changes are mediated by the stratosphere. Understanding the mechanisms underlying stratospheric interannual-to-decadal variability is important for detecting and attributing a possible ozone recovery as well as for assessing current and future climate change.
To study low-frequency variability in the stratosphere we have applied the new version of the chemistry-climate model (CCM) SOCOL (SOlar Climate Ozone Links)  to simulate the 20th century in transient mode with a set of nine ensemble members. Most CCM-based studies A. M. Fischer et al.: 20th Century Ensemble Simulations with the CCM SOCOL focus on the past 30 years, analyzing the increase in anthropogenic ozone depleting substances (ODSs), the effects of volcanic eruptions, and of solar variability (e.g. Shindell et al., 1998;Rozanov et al., 2002;Austin et al., 2003;Schnadt et al., 2002;Langematz et al., 2003;Rozanov et al., 2005). This is also one of the goals of the "REF1-simulations" carried out in the CCM Validation Activity (CCMVal) under the auspices of SPARC (Stratospheric Processes And their Role in Climate) (see Eyring et al., 2005;Eyring et al., 2006). The past 30 years, however, is too short to be representative of the full decadal variability of the stratosphere (see e.g., Brönnimann et al., 2004). It is also the period with largest anthropogenic influences due to emissions of ODSs and greenhouse gases (GHGs) and thus less suitable for studying natural (solar, volcanic, ENSO) variability. There are a few transient CCM simulations that go back to the 1960s or 1950s (Shindell et al., 1998;Dameris et al., 2005;Garcia et al., 2007;Austin and Wilson, 2006;Tegtmeier and Shepherd, 2007). Yet, none of these simulations were carried out with a large set of ensemble members.
One aspect of scientific interest concerns changes in stratospheric wave driving over the past decades and future climate. As a response to increasing GHG concentrations, middle atmosphere models, irrespective of coupling with interactive chemistry, predict a positive trend in the total amount of wave energy reaching the stratosphere . Conversely, Austin et al. (2003) report on CCM simulations over 1979-2001 suggesting a reduction in stratospheric wave drag. However, the exact mechanisms for the circulation change suggested by the models remain unclear. The spread amongst models is large and the trend over 1960-2000 is less clear than over future periods underscoring the need for ensemble member simulations to enhance the robustness of model results. Our simulations contribute to this discussion by allowing analysis of more simulations over a longer period. First results are shown in this paper.
While simulations covering the 20th century are standard for atmospheric general circulation models (AGCMs) (see also The Climate of the Twentieth Century (C20C) project, Folland et al., 2002;Scaife et al., 2008), this presents a challenge for CCMs as it requires additional data on boundary conditions as well as additional validation data. The necessary data, including information on the quasi-biennial oscillation (QBO) (Brönnimann et al., 2007), upper-level geopotential height and temperature fields (Griesser et al., 2008), as well as historical total ozone (TOZ) data (Vogler et al., 2006Staehelin et al., 1998), are only now becoming available. To simulate the 20th century with SOCOL in transient mode we have specified solar irradiance, stratospheric aerosols, sea surface temperatures (SST) and sea ice (SI), emission of GHGs, ODSs, nitrogen oxides (NO x ) and carbon monoxide (CO), nudged QBO and time evolving land surface changes at the model's boundaries. The focus of the paper is on boundary conditions and analysis of long-term behavior in stratospheric ozone and dynamics.
The paper is organized as follows: Section 2 briefly describes the updated version of SOCOL, the compilation of boundary conditions, experimental model setup and observational data. Results of the long-term model performance in the stratosphere are shown in Sect. 3. Our conclusions and outlook can be found in the last section.

Descriptions of new model version, boundary conditions and observational data for comparison
A transient simulation over the 20th century  was carried out in ensemble mode (9 ensemble members) using the CCM SOCOL .

Model description
The version of SOCOL used in the present work is essentially SOCOL vs2 described in detail by Schraner et al. (2008). In brief, SOCOL vs1 was first described and validated by Egorova et al. (2005). It is a combination of the Middle Atmosphere version of the European Centre/HAMburg Model 4 (MA-ECHAM4) spectral AGCM (Manzini and McFarlane, 1998) and the chemistry-transport model (CTM) MEZON (Model for Evaluation of oZONe trends) (Egorova et al., 2003). MA-ECHAM4 is a spectral model with T30 horizontal truncation and 39 vertical levels on a hybrid sigmapressure coordinate system. The model spans the atmosphere from the surface to 0.01 hPa. Chemical substances are transported by the hybrid numerical advection scheme (Zubov et al., 1999), which consists of the Prather scheme in the vertical direction (Prather, 1986) and of the Semi-Lagrangian scheme in the horizontal direction (Williamson and Rasch, 1989). Through its flux-form, the Prather scheme is mass conservative enabling the maintenance of strong vertical gradients, whereas the Semi-Lagrangian scheme is not mass conserving. SOCOL was used for the CCMVal REF1 simulation from 1975-2000 Eyring et al., 2006).
The CCM intercomparison of Eyring et al. (2006) revealed several shortcomings in SOCOL compared to observations and other CCMs. A number of modifications of the transport and chemistry/micro-physics scheme in SOCOL have been applied to overcome these problems leading to SOCOL vs2 . This includes an improved parameterization of stratospheric water vapor condensation, extension of chemical transport to all species and a more sophisticated heterogeneous chemistry scheme. A major deficiency of the original model version concerned artificial mass redistribution of several chemical species by a mass fixer after each horizontal transport step. The mass fixer guarantees mass conservation on a given layer by adjusting the whole field and penalizing those regions with the steepest horizontal gradient. However, Schraner et al. (2008) showed that this procedure led to artificial mass loss or accumulation in particular regions (i.e. when the transport error is regionally confined). To overcome this, a family-based mass fixing procedure for species of the nitrogen, chlorine, and bromine families was introduced. In case of ozone, where a family approach could not be pursued, sensitivity runs by Schraner et al. (2008) revealed that the transport error in ozone is mainly produced in a region between 40 • S to 40 • N, but predominantly corrected at the high latitudes. To avoid artificial mass transport by the mass fixer the fixer correction is only applied to 40 • S to 40 • N in the latest model version. By following this methodology global mass conservation on the model layers is still ensured.
These modifications led to considerable improvements in key chemistry-climate parameters, most notably in highlatitude stratospheric chlorine (Fig. 1) and ozone. Compared to observations as well as to other CCMs, SOCOL has significantly improved in the region of the southern polar vortex. For a detailed description of the modifications to and results from the latest model version we refer to Schraner et al. (2008). In this study, in addition to the nine ensemble member simulations performed with SOCOL vs2, we carried out one model simulation without the areal restriction (40 • S to 40 • N) of the mass fixing-procedure for ozone in order to analyze the effects of this implementation (see below).

Boundary conditions
Monthly and annually changing SST and SI distributions were obtained from the Hadley Centre HadISST data set (Rayner et al., 2003). For the description of solar variability we used spectral solar irradiance data compiled by Lean (2000) (Fig. 2a), which we applied to calculate the time evolution of solar heating and photolysis rates in the model.

Emissions
GHGs (Fig. 2b) and organic chlorine and bromine containing source gases (ODSs) are prescribed in SO-COL in the 5 lowermost model layers (corresponding to the lowermost ∼1.5 km in the model). Figure 2c shows the time development of equivalent effective stratospheric chlorine (EESC) derived from ODS mixing ratios.
Sixteen organic chlorine and bromine source gases are taken into account for the representation of ODSs. Values based on the CCMVal website (see above) for 1959-1999 were complemented with CH 2 Br 2 and CHBr 3 , two short-  Schraner et al., 2008).
lived bromocarbons emitted naturally from the ocean, for which we assumed a constant contribution of 1.63 pptv and 1.21 pptv, respectively, throughout the century (see also Warwick et al., 2006). Prior to 1959, we extended the ODSs time series with data from WMO (2003) where available. For methychloride and methylbromide a linear atmospheric trend was applied during the 20th century inferred from Antarctic firn air measurements (Butler et al., 1999;Sturges et al., 2001). Surface emissions of CO and NO x are represented by fluxes similar to those used for the CTM MOZART-2, representing values of the year 1990 (Horowitz et al., 2003). Natural sources were assumed to remain constant during the simulated period, while the time dependence of anthropogenic emissions was taken into account using data from EDGAR-HYDE 1.3 (van Aardenne et al., 2001), given on a 1 • ×1 • grid at 10 year intervals. These data include emissions from fossil fuel consumption, biofuel combustion, industrial processes, agricultural land, savannah burning, deforestation, and agricultural waste burning. The seasonal cycle of MOZART-2 was superimposed on the total emissions. Time-dependent aircraft NO x emissions are identical to those used for the CCM E39/C (Dameris et al., 2005), regridded to the SOCOL grid. They are based on a data set by Schmitt and Brunner (1997).

Aerosols
The time evolution of stratospheric aerosols is prescribed by the GISS data set described in Sato et al. (1993) with later updates available on http://data.giss.nasa.gov/modelforce/ strataer/ (see also Fig. 2d). For the period 1850-1999 these data contain monthly resolved and zonally averaged optical depths at 550nm for four altitude bands from 15 to 35 km   Lean, 2000), (b) Greenhouse gas concentrations for CO 2 (solid, in ppmv), N 2 O (dashed, in ppbv) and CH 4 (dotted, in ppbv), (c) Ozone Depleting Substances shown in equivalent effective stratospheric chlorine (EESC, in pptv), (d) Global mean extinction by stratospheric aerosols at 550 nm and an altitude of 22.5 km (Ext. in 1/km) (see Sato et al., 1993), (e) Vegetation classes of HYDE (see  for the year 1900 (above) and 1990 (below). height and a latitude grid spacing of 8 • . Reference to the GISS data set is less accurate than using directly the gapfilled SAGE aerosol record of Thomason and Peter (2006), as was done by Schraner et al. (2008). However, that record dates back only to the late 1970s, which is why we used the GISS data set. Aerosol optical depths are estimated from optical extinction data for 1882-1990, whose quality increases with time. Since 1979 data are complemented with satellite observations of the Stratospheric Aerosol and Gas Experi-ment (SAGE I/II) and Stratospheric Aerosol Measurement (SAM II) (Sato et al., 1993).
The monthly stratospheric aerosol properties of SOCOL (as zonal cross sections) are prescribed with surface area densities as well as with extinction coefficients, asymmetry factors and single-scattering albedos for all 8 spectral bands (2 shortwave and 6 longwave bands). The computation of these parameters was performed offline using optical depths and effective radii of GISS as input parameters. Effective radii during major volcanic eruptions before 1980 were reconstructed based on observations during the eruptions of El Chichon and Mt. Pinatubo. These data were also obtained from the GISS website. The computation of aerosol properties is based on Mie theory assuming an H 2 SO 4 concentration in the particle of 70 wt%, a width of a log-normal distribution of 1.8 and using temperatures from ERA40.
After the model runs discussed here were performed, an analysis of the data set revealed an error in the singlescattering albedo, which was set to 1.0 (scattering only) instead of around 0.995 in the near infrared band for the Mie calculation. This results in a slight underestimation in heating due to the presence of stratospheric aerosols. However, outside of the volcanically most perturbed periods this error will not have any significant effect on the model results.
Tropospheric aerosols are prescribed by a climatology described by Lohmann et al. (1999), used for the calculation of local heating rates and shortwave backscatter.

Quasi-biennial oscillation
The oscillation of stratospheric equatorial winds is nudged in SOCOL according to Giorgetta (1996). The prescribed zonal mean zonal wind field includes 19 pressure levels (from 90 hPa up to 3 hPa). From 1957 to 1999 data were obtained from ERA40 reanalyses and back to mid-1953 from the Free University of Berlin. Prior to that, we made use of a QBO-reconstruction that is based on historical pilot balloon wind data and on the QBO signature in the semidiurnal surface pressure oscillation which was extracted from historical sea-level pressure data. The reconstruction was validated with an independent reconstruction based on historical TOZ data, and it was found that after around 1910, the phase of the QBO in the middle stratosphere during peak phases should be correct around 80% of the time. Details on the final reconstruction and validation are given by Brönnimann et al. (2007).

Land surface changes
Effects of anthropogenic vegetation changes have been taken into account using the HYDE "A" data set , which is compiled on a 0.5 • by 0.5 • grid. It includes pasture and cropland areas from 1700 to 1950 in 50 year intervals and additionally for 1970 and 1990. Wherever possible Goldewijk (2001) organized the allocation of cropland and pasture on a country level using historical population density maps. Figure 2e illustrates the spread of pasture and cropland areas over the globe from 1900 to 1990 (violet colors). In order to assign the historical vegetation maps to ECHAM4, we mapped the 18 vegetation classes to the ECHAM vegetation classes (Claussen et al., 1994;Hagemann et al., 1999;Hagemann, 2002) and interpolated them to the model grid.
Vegetation parameters in ECHAM4 (Roeckner et al., 1996) are implemented following estimates by Claussen et al. (1994), prescribing land surface characteristics as a global field for present-day conditions. For our purpose, we follow later updates by Hagemann et al. (1999) and Hagemann (2002), who made use of improved estimates of the dependence of vegetation parameters on different vegetation classes as well as of a much higher resolved global distribution of major ecosystem types (Olson, 1994).
Using the mapped and regridded historical vegetation maps we assigned tabulated estimates of surface background albedo, surface roughness length, vegetation ratio, forest coverage, leaf area index, soil water capacity and volumetric wilting point to each vegetation class (Hagemann, 2002). The aggregation of surface roughness length to the coarser grid was not performed by linear upscaling but using an inverse logarithmic formulation (Claussen et al., 1994). As recommended by Hagemann et al. (1999), we replaced the constant value of volumetric wilting point (as in Roeckner et al., 1996) with a global spatially varying data set (Patterson, 1990). For years not included in the data set we interpolated the fields linearly in time.

Experimental design
For spin-up we used the off-line CTM version driven by daily winds, temperature and tropospheric water vapor representing the last year of the 25-year long time-slice experiment with CCM SOCOL described by Egorova et al. (2005). We carried out a 10 year run with this model using boundary conditions (NO x -and CO-emissions as well as ODSs and GHG concentrations) for 1890-1899. All mixing ratios, except tropospheric water vapor, were set to zero for the initialization. The coupled simulation started in 1900 with a one-year spin-up of the GCM part. In the second year of the coupled simulation, we perturbed the CO 2 concentration during one month in the range of 0.01% to obtain the nine different ensemble members, which after 1901 were driven by identical boundary conditions. The simulations were run until December 1999.

Comparison to observational data
Our analysis in the present study focuses on ozone and dynamics for which a wealth of observational data are available. The chemical performance of SOCOL vs2 was extensively validated by Schraner et al. (2008). Key tropospheric climate parameters for the 20th century have been analyzed in SOCOL within the framework of the CLIVAR C20C Project. This intercomparison revealed a good match with observations as well as with other atmospheric general circulation models on interannual and decadal time scales.
Modeled TOZ is compared against three observational data sets, the National Institute of Water and Atmospheric Research (NIWA) (Bodeker et al., 2005) combined database, data from the Backscatter Ultraviolet (BUV) instrument on the Nimbus 4 satellite in the early 1970s (see also Stolarski et al., 1997), and ground-based TOZ measurements prior to 1970, obtained from the World Ozone and Ultraviolet Radiation Data Centre (WOUDC). We contrast ozone distributions to the newly reprocessed Solar Backscatter Ultraviolet (SBUV) data set (version 8) (Rosenfield et al., 2005) and to a combined data set of SAGE-I/SAGE-II (55 • S-55 • N) supplemented with high latitude ozonesonde measurements at Syowa (69 • S) and Resolute (75 • N) providing global coverage from 10 to 50 km altitude (referred to as SAGE data set; for details see Randel and Wu, 2007). Prior to 1970 we compare SOCOL results to ozonesonde measurements from several European and North American stations, obtained from WOUDC. In addition, to estimate stratospheric ozone trends we use the Candidoz Assimilated Three-dimensional Ozone (CATO) data set (Brunner et al., 2006a). This data set, providing global coverage, was reconstructed by assimilating TOZ satellite observations, corrected for offsets and drifts against the Dobson network, into an equivalent latitude -potential temperature framework.
To evaluate the model's dynamics, modeled zonal wind and zonal temperature are compared with the ERA40 reanalysis results (Uppala et al., 2005). The northern subtropical jet index (referred to as "SJI" and measured as the maximum zonal mean zonal wind at 200 hPa in the 0 • N-50 • N belt) and polar vortex strength (referred to as "PVS" and measured as the difference between zonal mean geopotential height of the 100 hPa level averaged over 75 • N-90 • N and 40 • N-55 • N, respectively) were reconstructed as in Brönnimann et al. (2006a), i.e. using the same predictor data and reconstruction method. The reconstructions cover the period 1922 to 1947 and were merged with NCEP reanalyses (Kistler et al., 2001) thereafter. Eliassen-Palm (EP) fluxes were calculated from ERA40 and NCEP data by I. Wohltmann and are based on formulations by Andrews et al. (1987). Note that in the following we refer to reanalysis results of ERA40 and NCEP as "observations", acknowledging that they are the output from a model based assimilation of observations. The increasing concentrations of ODSs in the atmosphere (see Fig. 2c) initially led to depleted TOZ only in the region south of 70 • S for August to October (climatology over 1961-1980). Compared to BUV data from 1970-1976 (in a period where ODS concentrations were still relatively low), SOCOL underestimates the TOZ peak in the Northern Hemisphere (NH) and the Southern Hemisphere (SH). For the more recent time period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999), modeled TOZ amounts dropped globally and hence altered the seasonal cy-cle at all latitudes in line with observations. However, the simulated ozone hole is somewhat deeper (values drop to 100 DU) and persists longer than in the NIWA climatology. In terms of latitudinal extent of the ozone hole, model and observations are in good agreement. In autumn over midlatitudes TOZ is generally too high compared to NIWA, leading to an underestimation of the seasonal cycle. In general, the simulated TOZ climatologies over the past two decades are comparable to those simulated by other CCMs described in Eyring et al. (2006).

Analysis of ozone
To address the evolution of TOZ in the nine different ensemble members and hence to qualitatively judge the internal variability of the simulations, we have analyzed TOZ anomalies back to 1970 with respect to a climatology over 1979 to 1999 in five latitude bands (see Fig. 4). This is compared with a special run with the ozone mass fixer applied to the global field (blue line) instead of the standard restriction to the latitude band 40 • S-40 • N as described above (for details of the mass fixing procedure see Schraner et al., 2008). For mid-latitudes and the tropics the low-frequency variability in all SOCOL ensemble members is in good agreement with the NIWA collection of observations. The TOZ reduction due to the volcanic eruptions of El Chichon (1982) and Mt. Pinatubo (1991) are visible in the tropical belt in the model and the observations. However, the modeled decay seems to be too fast. In polar regions the amplitude of interannual variations is larger than at lower latitudes. Consistent with observations, a negative trend due to increasing ODS concentrations is visible at all latitude belts.
A shortcoming of SOCOL is the reduced amplitude of the seasonal cycle poleward of 30 • in both hemispheres, which is independent of the detailed configuration of the mass fixer (compare blue and grey curves in climatologies of Fig. 4). Sensitvity tests with the mass fixer completely switched off for ozone revealed that this underestimation is at least partly caused by the mass fixing procedure . As can be seen in Fig. 4, the major effect of the mass fixer restriction for ozone to the 40 • S-40 • N band is to shift the whole seasonal cycle of TOZ downward (upward) in the tropics (at mid-to high latitudes). This is because the global layer transport error and its correction, respectively, is only applied to a confined area (i.e. where most of the Semi-Lagrangian transport error occurs, Schraner et al., 2008) instead of the global field leading to a more realistic seasonal amplitude especially in the tropics. We have compared other fields of the ensemble mean and the run with globally applied mass fixer (e.g. zonal wind, temperature, ozone), but did not find any irregularities that would disprove an areal restriction of the mass fixing procedure for ozone. Therefore, in the following we address the results of the nine ensemble simulations with the mass fixer restricted to 40 • S to 40 • N.
To validate modeled TOZ before 1970 we have used NH ground-based stations of Svalbard, Tromsø, Uppsala, Oxford and Arosa that allow longer-term analysis (starting between 1924 and 1955) and cover a wide latitude range (46 • N to  1901-1930 (upper left), 1931-1960 (upper right), 1961-1980 (middle left) and 1981-1999 (middle right). Measurement onboard satellites of BUV (1970-1976) and NIWA (1981 are shown in lower left and lower right panel, respectively. 77 • N, see Fig. 5). The series of Svalbard and Oxford have only recently been reevaluated (Vogler et al., 2006;Vogler et al., 2007) and provide an excellent data source to validate decadal scale transient simulations with a CCM. The ensemble members of SOCOL are in reasonable agreement with observations. This is especially true for Svalbard, but less so for the autumn months in Uppsala, Oxford, and Arosa. In general,the aforementioned model failures to capture the magnitude of the seasonal cycle is clearly visible at the different stations.
Similar to TOZ, zonally averaged vertical ozone also reveals little climatological change before the late 1970s. Conversely, over the past two decades ozone mixing ratios dropped significantly, e.g. for the month of January, by around 1 ppmv at 10 hPa for the latitudes north of 80 • N (see Fig. 6). Compared to observations of SBUV for the period of 1970-1976 and 1979-1999, the location and magnitude of the maximum ozone concentration in the tropics is well reproduced. However, model and observations diverge in the upper stratosphere (above the 10 hPa level), with largest discrepancies (up to 2 ppmv below observations) over the NH (SH) in boreal winter (summer). Below the 10 hPa level modeled ozone is too high (more than 1 ppmv), especially over the tropics and the SH (NH). This bias was already documented in Eyring et al. (2006) and Schraner et al. (2008) and is present in all months. The comparison of ozone number densities against observations by SAGE confirms these findings in the lower and upper stratosphere. It additionally shows an underestimation of ozone throughout the stratosphere in winter polar regions, leading to smaller than observed TOZ values (compare Fig. 3).
Prior to 1970 we validated the model runs with ozonesonde observations during the 1960s from Hohenpeissenberg, Uccle, Payerne, Boulder, and Resolute, e.g. Fig. 6e for the January intercomparison. Note, that a high quality for these old measurements is not assured and a validation to model results can therefore only be done in a qualitative way. For this comparison, only Brewer Mast measurements  1925 1930 1935 1940 1945 1950 1955 1960 1965 1970 Tromso (18.9E, 69.7N) [DU] 250 350 450 1925 1930 1935 1940 1945 1950 1955 1960 1965 1970 Uppsala ( with a correction factor between 0.9 and 1.35 (0.9-1.2 in case of Hohenpeissenberg) were considered (similar to Logan, 1994). Brewer Mast data are more reliable than Regener type of measurements which were frequently applied at that time (S. Oltmans, personal communication) but excluded in the present study. The European stations in Fig. 6e are at similar latitudinal location and shown together to obtain a more robust result.
In general, SOCOL is in qualitatively good agreement with these old station measurements. This is especially the case compared to European stations and Boulder throughout the year. In contrast to the Canadian station Resolute, which is located further north, modeled ozone mixing ratios in the lower stratosphere are comparable in winter but higher and outside the observed interannual standard deviation during summer (not shown). The model internal variability as well as the year-to-year fluctuation in station data is largest during winter time and reduces in summer.
To analyze the model runs in a process-oriented framework, we have tested whether the sensitivity of ozone upon  increasing concentrations of ODSs is correctly reproduced. We estimated ozone trends over all months from 1979-1999 by applying a multiple linear regression analysis with solar variability (see Fig. 2a), two orthogonal QBO time series at 10 and 30 hPa, (method based on Wallace et al., 1993) and extinction from stratospheric aerosols (at 550 nm, see Fig. 2d) as explanatory variables and zonally averaged ozone anomalies as target value.
O 3 (t) = b · EESC(t) + c · solar cylce + d 1 · QBO 10 +d 2 · QBO 30 + e · strat aer + ε(t), The trend coefficient b is regressed upon EESC time series (taken at the boundary, see Fig. 2c) and given in percentage per unit EESC. (Note that in a strict sense, the coefficient b is not itself a trend, but the sensitivity of ozone with respect to the trend prescribed by EESC, see Maeder et al. (2007). However, in the following for simplicity we refer to it as a trend.) Also note that percentage trends are calculated with respect to the average over the trend period. Autocorrelations within the residual time series ε(t), which affects the significance level of the trend estimate (p-value), were eliminated by applying a third-order autoregressive process. The linear regression model was applied to individual ensemble members (not shown) and ensemble mean simulation of SOCOL as well as to observational data sets of SBUV, CATO and SAGE (see Fig. 7). The general pattern of estimated trends in zonal mean ozone is comparable in all ensemble member simulations. The upper stratospheric ozone trends in the ensemble mean and individual members of SOCOL show two significant local minima around 40 km altitude of up to −10%/ppbv EESC. This is in excellent agreement with SBUV and SAGE data. Upper stratospheric negative trends in both hemispheres are also appearing in CATO but the magnitude is much lower and probably an artifact of the reconstruction . In the Antarctic lower stratosphere a downward trend of more than −20%/ppbv EESC can be observed in SAGE and CATO which is generally underestimated in SOCOL (−14%/ppbv EESC averaged over all members). However, the variability among the ensemble members in this region is larger and three ensemble members reach the magnitude of observed trends. Garcia et al. (2007) also reported a too small downward trend in the ensemble mean simulation with the Whole Atmosphere Community Climate Model (WACCM). Note, that the lower than observed TOZ values over southern high latitudes mentioned above is a result of a negative bias (due to remaining artificial problems of the transport scheme), present throughout the century, rather than an effect of a negative trend. In the Arctic lower stratosphere the ensemble mean of SOCOL shows a significant negative trend that is only somewhat smaller than in SAGE and CATO.

Analysis of dynamics
A realistic model representation of dynamical features is a prerequisite for a reliable simulation of ozone and other trace gases as well as for understanding the interaction between dynamics and chemistry. To assess the stratospheric model dynamics we first examined the well-known basic processes for the winter and summer hemispheres. Figure 8 illustrates the ensemble mean of zonal mean zonal wind for different climatic periods and in relation to ERA40 for the month of January and July. In SOCOL climatological zonal wind has hardly changed from 1901-1957 to 1958-1970. It simulates very well the main characteristics in both hemispheres and seasons. The northern winter polar vortex in SOCOL exceeds the observed wind speed by more than 7 m/s. This bias is less significant in the period 1958-1970 than in the more recent period . In the winter SH zonal winds in the middle to upper stratosphere are weaker than observed for 1958-1970 but stronger for 1971-1999, especially in a latitude band of 55 • S-75 • S. Upper stratospheric easterlies are simulated in both summer hemispheres in agreement with ERA40. However, the areal extent is confined to the midlatitudes in contrast to ERA40 where strong easterlies prevail over the whole hemisphere. A further discrepancy to ERA40 is a slight shift of the modeled subtropical jet towards the tropics in all periods (see also Brönnimann et al., 2006b). Its magnitude, however, looks reasonable compared to ERA40.
To further assess the model's performance in the stratosphere, we have analyzed zonal mean temperature (Fig. 9). Similar to zonal wind, no big alterations in modeled zonal temperature can be observed for the first seven decades. In line with a stronger northern polar vortex, modeled temperatures are colder at northern high latitudes between 100 and 10 hPa suggesting too little downwelling over high latitudes. In the Antarctic stratosphere, above 100 hPa, higher temperatures are simulated than in the reanalyses in the period 1958-1970, but to a lesser degree for 1971-1999 reflecting again the sensitivity to the polar vortex. In general, agreement between ERA40 and SOCOL is better for July than January where a cold bias is present in the lower and upper stratosphere over almost all latitudes. A cold bias in the upper stratosphere might be connected to a rather simplified absorption scheme of MA-ECHAM4 for solar UV radiation . Note that the upper stratospheric temperatures are subject to a large variability among different CCMs, as well as observational data (Eyring et al., 2006). Negative model biases near the tropopause can be explained by the coarse vertical resolution in the upper tropospherelower stratosphere . The modeled cold tropical tropopause temperature is well reproduced regarding location and magnitude and only slightly underestimated. In general, SOCOL is well able to capture the main dynamical characteristics in zonal mean wind and temperature back to 1901.
To address interannual variability in all nine ensemble members back to 1920, we have analyzed boreal winter time series of SJI and PVS (Fig. 10). Model results are compared to indices for the month of January to March (when reconstructions show sufficient skill, i.e. a reduction of error (RE) greater than 0.36) (see Brönnimann et al., 2006a). Note that a validation back to 1920 is only possible for the NH, where enough historical upper air station data is available. For the SH we compared the modeled PVS to ERA40 (not shown). In the case of SJI the correlation between ensemble mean simulation and the observed index is high (0.74) and internal variability within SOCOL rather low. SJI is strongly related to tropical SST and hence ENSO via changes in the Hadley Circulation. SOCOL simulates a very strong ENSO response and thus agrees extremely well with the reconstructed index. SOCOL captures comparably little of the year-to-year variability in the PVS reconstruction, although observations mostly lie within the large ensemble spread. The correlation between the ensemble mean and reconstruction is far lower (0.25) than for SJI and does not increase for individual ensemble members. Yet, it is interesting to note that two individual deviations from the mean are modeled fairly well: winter 1976 and 1987, which are possibly linked to anomalous ENSO events (Brönnimann et al., 2006b;Fischer et    2008). Similar to the NH, the ensemble variability in PVS of the SH is almost as large as interannual variability of the observations (not shown). Yet, the correlation between the ensemble mean and ERA40 is lower and drops to near zero.
To further elucidate the large model spread in the north polar region in TOZ and PVS, we have analyzed major warmings (MW), being a major source of interannual variability during boreal winter (see e.g., Matsuno, 1971). We have applied a standard WMO index with two criteria to diagnose MWs in the model and ERA40 for the months of November to March (see also Chaffey and Fyfe, 2001): 1) 10 hPa temperature difference between 60 • N and 90 • N is positive for four consecutive days (also referred to as "minor warming", note that due to horizontal model resolution we have extracted the fields at latitudes 61 • N and 85 • N) and 2) zonal mean zonal wind at 10 hPa and 60 • N changes from westerly to easterly flow within the minor warming period (also referred to as the "central date"). The central dates were counted separately in each ensemble member and analyzed upon frequency and seasonality (Fig. 11). Changes in zonal wind from westerly to easterly were counted only once within a minor warming period and possible final warmings in late March were excluded.
The total number of MWs for ERA40 over 1958-1999 amounts to 28, which is by a factor of 3-4 larger than in SOCOL averaged over all members for the same period. This is in agreement with the above mentioned findings of a too strong modeled northern polar vortex. On average, in 17% (66.6%) of analyzed years a MW is found in SOCOL (ERA40). For the period 1901-1957 the percentage is only slightly higher (22%). Dameris et al. (2005) also reported a too small number of modeled MWs in the E39/C model. Their model top, however, is at 10 hPa which hinders to fully capture the propagation of gravity waves and its interaction with the polar vortex.
A large fraction of MWs is normally observed in January/early February and in late February and mid-March. Five MWs occurred in late November to mid-December. Ensemble members diverge largely in the occurrence of MWs in early winter but show an increased frequency towards the end of winter on average in line with ERA40. In January the chance for a simulated MW, however, is generally low (similar to Dameris et al., 2005).
We defined four different parameters addressing the strength of individual MWs (Fig. 12a-d) in both periods, 1901-1957 and 1958-1999: a) Maximum zonal mean temperature at 85 • N and 10 hPa, b) Maximum zonal mean temperature at 61 • N and 10 hPa, c) Minimum zonal mean wind at 61 • N and 10 hPa, and d) Maximum temperature increase at 10 hPa and 85 • N within 5 days. In general the strength parameters in the ensemble members are in reasonable agreement with ERA40. The ensemble mean warming lies within the range of the first and third quartile of ERA40 (see thin horizontal lines in Fig. 12) and is hardly altered between the two time periods. For temperature at 85 • N (Fig. 12a), SO-COL shows a slight positive bias for 1958-1999 compared to ERA40, whereas agreement is better for the latitude 61 • N. The magnitude of average minimum easterly winds is somewhat lower in SOCOL than in ERA40. Good agreement is found for the fourth strength parameter (Fig. 12d), where temperature increases in the range of 10 to more than 60 K per 5 day window are simulated and observed.
Perturbations on the polar vortex and resulting MWs are caused by interaction of propagating waves, which is analyzed in the following. We have first evaluated the vertical component of the 100 hPa EP flux (EPz), averaged over 40 • N-80 • N and 40 • S-80 • S respectively, over the century to see whether the model captures the total midwinter wave activity that propagates from the troposphere into the stratosphere on an interannual-to-decadal scale. Time series of EPz (for SOCOL calculated from 12-hourly frequency output) during winter months (January to February and July to August, respectively) are shown in Fig. 13a. The overall magnitude of simulated EPz looks reasonable compared to observations of NCEP and ERA40 in both hemispheres. In the NH some individual developments in the averaged SO-COL time series appear as in observations: e.g., 1966, 1970-1975, 1983-1987. At least part of these anomalies are related to warm ENSO events which increase the vertical wave propagation in SOCOL    1901-1957 (red) and 1958-1999 (black). The green (red) shading in the upper left panel denotes the averaged number of model occurrences for 1958-1999 (1901-1957).
could be an artifact of the data assimilation, since more observational (satellite) data could be assimilated in later years which may influence the representation of waves. To overcome this, we have compared climatologies of EP flux diagnostics between the ensemble mean and ERA40 over the period 1990-1999 (not shown). In the NH (SH) modeled EPz is lower than ERA40 from 30 • N-60 • N throughout the stratosphere (from 50 • S-80 • S up to 20 hPa). Note that this negative bias is caused by averaging out some of the ensemble member variability, but is still consistent with Fig. 13. In addition, wave mean flow interaction measured by the di-vergence of the EP flux is attenuated in SOCOL (less negative) north of 40 • N (in the SH the signal is also positive but less dominant). This may explain the stronger than observed polar vortex and the reduced number of MWs seen above. However, the residual meridional circulation (RMC), analyzed by the Transformed Eulerian Mean equations (Andrews et al., 1987) and the tropical tape recorder reveal a less clear signal. The RMC is only decelerated in a region from mid-to high latitudes up to 10 hPa and not uniformly over the whole hemisphere as would be expected from reduced wave mean flow interaction. In the SH the RMC is  weaker in the upper stratosphere but stronger in the middle and lower stratosphere. In addition, the tropical tape recorder exhibits too strong upwelling compared to observations (not shown). Yet, the estimated upward velocities are not constant in altitude but rather show stronger windspeed from 100 hPa to 40 hPa and reduced windspeed from 40 hPa and 10 hPa, similar to Steil et al. (2003) for simulation with MA-ECHAM4/CHEM model.
Trends of EP flux and the RMC are presently a widely discussed topic of the scientific community (see e.g., Butchart et al., 2006). We have calculated trends of EPz (of Fig. 13a) within the modeled and observed time series for different time windows to account for decadal variability (Fig. 13b,  Tables 1 and 2). Reanalyses show a positive NH EPz trend of around 2%/decade (not significant) over 1958-1979 and a negative trend over 1980-1999. Model realizations over the same period simulate trends of opposite signs (Table 1) 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 1901-1929 1930-1959 1901-1959 1958-1979 1980-1999 1958-1999 1901-1929 1930-1959 1901-1959 1958-1979 1980-1999 1958-1999 a) b) which are mostly not significant. Over 1980-1999 the model realizations do not capture the negative trend as in the reanalyses. Note that most CCMs show a smaller than observed negative trend for the period 1979-2001 (Austin et al., 2003) which is not significant in most models. Model results before 1958 are somewhat more conclusive: seven out of nine members show a negative trend over 1930-1959 and a positive trend over 1901-1929 for the whole set of ensemble members. Yet, statistical significance on the 5% level is again almost absent in all ensemble members.
In the SH, for 1958-1999, SOCOL features a large model spread with trend estimates of both signs compared to very strong positive trends in the reanalyses (Table 2), which could be affected by data quality (see above). Note that ERA40 and NCEP show rather large dissimilar trends in the SH. One realization shows a significant positive trend (on the 5% level) of 7.5%/decade over 1958-1999. This is somewhat smaller than ERA40. Before 1958, all SOCOL mem-bers (two of them significantly on the 5% level) exhibit negative EPz trends for the trend period 1901-1959 similar to the NH.
Propagating waves affect the state of the winter polar vortex and related variables, such as temperature and TOZ. High correlation between meridional heat flux (proportional to EPz in the stratosphere) and 50 hPa polar temperatures has been documented by Newman et al. (2001) and addressed in several CCM-comparison studies (Eyring et al., 2006;Austin et al., 2003). Here we evaluate this process by showing scatter plots of January to February (July to August) 40 • −N-80 • N (40 • S-80 • S) EPz and February to March (August to September) 50 hPa zonal temperature averaged over 60 • N-90 • N (60 • S-90 • S) (Fig. 14a). In the NH the linear fit between modeled polar temperatures and EPz agrees extremely well with NCEP and ERA40 for the period 1958-1999 and also for 1901-1957. This is also true for the SH compared to NCEP. Differences in the intercept reflect deviations in Table 1. NH trends (Beta, in %/decade) and standard error (se) in the 100 hPa vertical component of EP flux, averaged over 40 • N-80 • N and over January to February for different periods. Light grey (dark grey) shading marks statistically significant estimates with p-value <0.05 (0.1). 1901-1929 1930-1959 1901-1959 1958-1979 1980-1999 1958-1999 1901-1929 1930-1959 1901-1959 1958-1979 1980-1999 1958-1999

Summary and discussion
In the first seven decades of the 20th century with still modest anthropogenic influence, modeled climatologies remain almost unaltered regarding TOZ, vertical ozone distribution, zonal mean zonal wind, and zonal temperature. SOCOL is able to capture TOZ development well compared to historic ground-based stations back to 1924, but the amplitude of the seasonal cycle is underestimated which is a remaining shortcoming of SOCOL vs2. Compared to satellite observations back to 1970 the low-frequency variability in TOZ is well reproduced over the tropics and mid-latitudes. In the same region back to 1920, SJI in SOCOL is highly correlated (0.74) with the reconstructed index. Its interannual variability reflects to a large degree ENSO variability and corresponding changes in the Hadley and Walker circulations. The A. M. Fischer et al.: 20th Century Ensemble Simulations with the CCM SOCOL SJI allows validating model results in an insightful, processoriented manner over long time-scales. It could potentially be further explored by other middle-atmosphere studies to test whether the response to oceanic forcing is correctly reproduced. Zonal mean vertical ozone distributions are in reasonable agreement with SAGE and SBUV. Yet, they show a positive bias in the lower stratosphere and a negative bias in the upper stratosphere. Modeled stratospheric trends over 1979-1999 perform well compared to observations for all ensemble members. Trend estimates among the ensemble members vary most in the southern lower stratosphere where on average the observed negative trend is underestimated.
In general, the polar regions exhibit a large internal model variability. This is revealed by the TOZ and PVS patterns suggesting that a sufficiently large set of ensemble members will be required if further progress in model development is sought for in these particular geographic regions. For PVS in the NH and SH the ensemble spread is almost as high as interannual variability in the observed index (except for some ENSO related events) and hence not predictable. MWs in the NH being a major source of interannual variability were analyzed with respect to their frequency, seasonality and intensity. A much smaller number (factor 3-4 smaller) of MWs occurred in the model than in observations, though the seasonal distribution of the occurrence frequency looks comparable, except for January. No discrepancies to ERA40 were found for analyzed strength parameters of individual MWs. This means that SOCOL underestimates the total number of MWs but it evolves in a comparable way in case a MW occurs. The underestimation of MWs indicates a too strong northern polar vortex which is further endorsed by higher zonal wind speed and colder temperatures over the Arctic stratosphere compared to ERA40. EPz and EP flux divergence reveal less wave activity and less wave mean flow interaction during winter months for 1990-1999 in the northern high latitudes and hence less perturbations of the polar vortex. In the SH the signal is less clear.
Time series of modeled EPz in the NH and SH exhibit a large internal variability as well as interannual variability. Modeled trends of EPz are generally not significant and the spread of estimates between different ensemble members is high. More sophisticated techniques are required to further explore the robustness of the trend estimates. The correlation between EPz and polar temperature as well as EPz and polar spring TOZ was examined, which shows good agreement with observations for both hemispheres and reveal comparable relationships for the pre-ERA40 and pre-satellite period.

Conclusion and outlook
We have applied the CCM SOCOL vs2  with a set of 9 ensemble members covering the 20th century in transient mode using prescribed SSTs, SI distri-bution, NO x -and CO-emissions, GHG and ODSs concentrations, QBO, and changes in land properties.
Despite its rather coarse horizontal and vertical resolution (3.75 • ×3.75 • (T30), 39 layers, 1000-0.01 hPa), simulations with the CCM SOCOL perform rather well compared to observations and allow addressing interannual-to-decadal variability during the twentieth century in relation to its forcing parameters. The model well reproduces observed ozone anomalies and reconstructed dynamics in the (sub-) tropical belt. That means that in this region the response of SOCOL to forcings of QBO, ENSO, solar variability and volcanic eruptions is very realistic and predictable even with less than nine ensemble members. Also, process-oriented analysis of model output (i.e. positive correlation of polar TOZ with EPz, positive correlation of polar temperature with EPz and negative ozone trend estimates upon EESC) reveal only little changes between different ensemble members and hence represent robust model results. On the contrary, our analysis of the polar vortex systems (i.e. PVS) show that in these areas the ensemble spread is very large which hinders to separate between internal and externally forced variability. The large ensemble range found for estimated trends in stratospheric wave drag, including estimates of opposite signs, reveals that the internal variability in models cannot be neglected when predicting the response in wave energy upon external forcings.
The presented model simulations provide a unique data set for further exploration of interannual-to-decadal variability. In particular, it is planned to quantify the impact of solar variability, ENSO variability and volcanic eruptions upon stratospheric ozone and climate individually over the whole twentieth century. The interaction between low-frequency variability in the stratosphere and tropospheric climate modes will be addressed in addition.