The 11-year solar cycle in current reanalyses : A ( non ) linear attribution study of the middle atmosphere

This study focusses on the variability of temperature, ozone and circulation characteristics in the stratosphere and lower mesosphere with regard to the influence of the 11-year solar cycle. It is based on attribution analysis using multiple nonlinear techniques (Support Vector Re5 gression, Neural Networks) besides the multiple linear regression approach. The analysis was applied to several current reanalysis datasets for the 1979-2013 period, including MERRA, ERA-Interim and JRA-55, with the aim to compare how this type of data resolves especially the double10 peaked solar response in temperature and ozone variables and the consequent changes induced by these anomalies. Equatorial temperature signals in the lower and upper stratosphere were found to be sufficiently robust and in qualitative agreement with previous attribution studies. The anal15 ysis also pointed to the solar signal in the ozone datasets (i.e. MERRA and ERA-Interim) not being consistent with the observed double-peaked ozone anomaly extracted from satellite measurements. Consequently the results obtained by linear regression were confirmed by the nonlinear approach 20 through all datasets, suggesting that linear regression is a relevant tool to sufficiently resolve the solar signal in the middle atmosphere. Furthermore, the seasonal evolution of the solar response was also discussed in terms of dynamical causalities in the winter hemispheres. The hypothetical mechanism 25 of a weaker Brewer Dobson circulation at solar maxima was reviewed together with discussion of polar vortex behaviour.


Introduction
The Sun is a prime driver of various processes in the climate system. From observations of the Sun's variability on decadal or centennial timescales, it is possible to identify temporal patterns and trends in solar activity, and consequently to derive the related mechanisms of the solar influence on the Earth's climate (e.g. Gray et al., 2010). Of the semi-regular solar cycles, the most prominent is the approximate 11-year periodicity which manifests in the solar magnetic field or through fluctuations of sunspot number, but also in the total solar irradiance (TSI) or solar wind properties. For the dynamics of the middle atmosphere, where most of the ozone production and destruction occur, the changes in the spectral solar irradiance (SSI) are the most influential, since the TSI as the integral over all wavelengths exhibits variations of orders lower than the ultraviolet part of the spectrum (Lean, 2001). This fact was supported by original studies (e.g. Labitzke, 1987;Haigh, 1994) that suggested the solar cycle (SC) influence on the variability of the stratosphere. Gray et al. (2009) have shown, with the fixed dynamical heating model, that the response of temperature in the photochemically controlled region of the upper tropical stratosphere is due to both direct solar heating and an indirect effect caused by the ozone changes.
Numerous studies have identified temperature and ozone changes linked to the 11-year cycle by multiple linear regression. The use of ERA-40 reanalysis (Frame and Gray, 2010) pointed to a manifestation of annually averaged solar signal in temperature, exhibited predominantly around the Equator with amplitudes up to 2 K around the stratopause and with a secondary amplitude maximum of up to 1 K in the lower stratosphere. Soukharev and Hood (2006), Hood et al. (2010) A. Kuchar et al.: Solar cycle in current reanalyses and Randel and Wu (2007) have used satellite ozone data sets to characterise statistically significant responses in the upper and lower stratosphere. The observed double-peaked ozone response in the vertical profile around the Equator was reproduced in some chemistry climate models, although concerns about the physical mechanism of the lower stratospheric response were expressed (Austin et al., 2008).
The ozone and temperature perturbations associated with the SC have an impact on the middle atmospheric circulation. They produce a zonal wind anomaly around the stratopause (faster subtropical jet) during solar maxima through the enhanced meridional temperature gradient. Since planetary wave propagation is affected by the zonal mean flow (Andrews and McIntyre, 1987), we can suppose that a stronger subtropical jet can deflect planetary waves propagating from higher latitudes. Reduced wave forcing can lead to decreasing/increasing or upwelling/downwelling motions in the equatorial or higher latitudes respectively (Kodera and Kuroda, 2002). The Brewer-Dobson circulation (BDC) is weaker during solar maxima (Kuroda and Kodera, 2001), although this appears to be sensitive to the state of the polar winter. Observational studies, together with model experiments (e.g. Matthes et al., 2006), suggest a so-called "topdown" mechanism where the solar signal is transferred from the upper to lower stratosphere, and even to tropospheric altitudes.
Statistical studies (e.g. Labitzke et al., 2006;Camp and Tung, 2007) have also focused on the lower stratospheric solar signal in the polar regions and have revealed modulation by the Quasi-Biennial Oscillation (QBO), or the well known Holton-Tan relationship (Holton and Tan, 1980) modulated by the SC. Proposed mechanisms by Matthes et al. (2004Matthes et al. ( , 2010 suggested that the solar signal induced during early winter in the upper equatorial stratosphere propagates poleward and downward when the stratosphere transits from a radiatively controlled state to a dynamically controlled state involving planetary wave propagation (Kodera and Kuroda, 2002). The mechanism of the SC and QBO interaction, which stems from reinforcing each other or canceling each other out , has been verified by WACCM3.1 model simulations . These proved the independence of the solar response in the tropical upper stratosphere from the response dependent on the presence of the QBO at lower altitudes. However, fully coupled WACCM-4 model simulations by Kren et al. (2014) raised the possibility of occurrence by chance of the observed solar-QBO response in the polar region. The internally generated QBO was not fully realistic though. In particular, the simulated internal QBO descended down to only about 50 hPa.
It has been shown that difficulties in the state-of-the-art climate models arise when reproducing the solar signal influence on winter polar circulation, especially in less active sun periods (Ineson et al., 2011). The hypothesis is that solar UV forcing is too weak in the models. Satellite measurements indicate that variations in the solar UV irradiance may be larger than previously thought (Harder et al., 2009). However, the measurements by Harder et al. (2009) from the SORCE satellite may have been affected by instrument degradation with time and so may be overestimating the UV variability (see the review by Ermolli et al., 2013). The latter authors have also concluded that the SORCE measurements probably represent an upper limit on the magnitude of the SSI variation. Consequent results of general circulation models, forced with the SSI from the SORCE measurements, have shown a larger stratospheric response than for the NRL SSI data set. Thus, coordinated work is needed to have reliable SSI input data for GCM and CCM simulations (Ermolli et al., 2013), and also to propose robust conclusions concerning SC influence on climate (Ball et al., 2014).
At the Earth's surface, the detection of the SC influence is problematic since there are other significant forcing factors, e.g. greenhouse gases, volcanoes and aerosol changes (e.g. Chiodo et al., 2012), as well as substantial variability attributable to internal climate dynamics. However, several studies (van Loon et al., 2007;van Loon and Meehl, 2008;Hood and Soukharev, 2012;Hood et al., 2013;Gray et al., 2013;Scaife et al., 2013) detected the solar signal in sea level pressure and sea surface temperature, which supports the hypothesis of a troposphere-ocean response to the SC. Some studies (e.g. Hood and Soukharev, 2012) suggest a so-called "bottom-up" solar forcing mechanism that contributes to the lower stratospheric ozone and temperature anomaly in connection with the lower stratosphere deceleration of the BDC.
The observed double-peaked ozone anomaly in the vertical profile around the Equator was supported by the simulations of coupled chemistry climate models (Austin et al., 2008). However, the results presented by Chiodo et al. (2014) suggest the contribution of SC variability could be smaller since two major volcanic eruptions are aligned with solar maximum periods and also given the shortness of the analysed time series (in our case, 35 years). These concerns related to the lower stratospheric response of ozone and temperature derived from observations have already been raised (e.g. Solomon et al., 1996;Lee and Smith, 2003). However, another issue is whether or not the lower stratospheric response could depend on the model employed in the simulations (Mitchell et al., 2015b).
Several past studies (e.g. Soukharev and Hood, 2006;Frame and Gray, 2010;Gray et al., 2013;Mitchell et al., 2014) used multiple linear regression to extract the solar signal and separate other climate phenomena like the QBO, the effect of aerosols, North Atlantic Oscillation (NAO), El Niño-Southern Oscillation (ENSO) or trend variability. Apart from this conventional method, it is possible to use alternative approaches to isolate and examine particular signal components, such as wavelet analysis (Pisoft et al., 2012(Pisoft et al., , 2013 or empirical mode decomposition (Coughlin and Tung, 2004). The nonlinear character of the climate system also suggests potential benefits from the application of fully nonlinear attribution techniques to study the properties A. Kuchar et al.: Solar cycle in current reanalyses 6881 and interactions in the atmosphere. However, such nonlinear methods have been used rather sporadically in the atmospheric sciences (e.g. Walter and Schönwiese, 2003;Pasini et al., 2006;Blume and Matthes, 2012), mainly due to their several disadvantages such as the lack of explanatory power (Olden and Jackson, 2002).
To examine middle atmospheric conditions, it is necessary to study reliable and sufficiently vertically resolved data. Systematic and global observations of the middle atmosphere only began during the International Geophysical Year (1957Year ( -1958 and were later expanded through the development of satellite measurements (Andrews and McIntyre, 1987). Supplementary data come from balloon and rocket soundings, though these are limited by their vertical range (only the lower stratosphere in the case of radiosondes) and the fact that the in situ observations measure local profiles only. By assimilation of these irregularly distributed data and discontinuous measurements of particular satellite missions into an atmospheric/climatic model, we have modern basic data sets available for climate research, so-called reanalyses. These types of data are relatively long, globally gridded with a vertical range extending to the upper stratosphere or the lower mesosphere and thus suitable for 11-year SC research. In spite of their known limitations (such as discontinuities in ERA reanalysis -McLandress et al., 2014), they are considered an extremely valuable research tool (Rienecker et al., 2011).
Coordinated intercomparison has been initiated by the SPARC (Stratospheric Processes and their Role in Climate) community to understand them, and to contribute to future reanalysis improvements (Fujiwara et al., 2012). Under this framework, Mitchell et al. (2014) have examined nine reanalysis data sets in terms of 11-year SC, volcanic, ENSO and QBO variability. Complementing their study, we provide here a comparison with nonlinear regression techniques, assessing robustness of the results obtained by multiple linear regression (MLR). Furthermore, EP flux diagnostics are used to examine solar-induced response during the winter season in both hemispheres, and solar-related variations of assimilated ozone are investigated.
The paper is arranged as follows. In Sect. 2 the used data sets are described. In Sect. 3 the analysis methods are presented along with regressor terms employed in the regression model. Section 4 is dedicated to the description of the annual response results. In Sect. 4.1 solar response in MERRA reanalysis is presented. Next, in Sect. 4.1.1 other reanalyses are compared in terms of SC. Comparison of linear and nonlinear approaches is presented in Sect. 4.1.2. Section 4.2 describes monthly evolution of SC response in the state variables. Section 5 is aimed at dynamical consequences of the SC analysed using the EP flux diagnostics.

Data sets
Our analysis was applied to the most recent generation of three reanalysed data sets: MERRA (Modern Era Reanalysis for Research and Applications, developed by NASA) (Rienecker et al., 2011), ERA-Interim (ECMWF Interim Reanalysis) (Dee et al., 2011) and JRA-55 (Japanese 55-year Reanalysis) (Ebita et al., 2011). We have studied the series for the period 1979-2013. All of the data sets were analysed on a monthly basis. The Eliassen-Palm (EP) flux diagnostics (described below) was computed on a 3-hourly basis from MERRA reanalysis and subsequently monthly means were produced. A similar approach has already been used by Seviour et al. (2012) and Mitchell et al. (2015a). The former study proposed that even 6-hourly data are sufficient to diagnose tropical upwelling in the lower stratosphere. The vertical range extends to the lower mesosphere (0.1 hPa) for MERRA, and to 1 hPa for the remaining reanalyses. The horizontal resolution of the gridded data sets was 1.25 • × 1.25 • for MERRA and JRA-55 and 1.5 • × 1.5 • for ERA-Interim respectively.
In comparison with previous generations of reanalyses, it is possible to observe a better representation of stratospheric conditions. This improvement is considered to be connected with increasing the height of the upper boundary of the model domain (Rienecker et al., 2011). For example, the Brewer-Dobson circulation was markedly overestimated by ERA-40; an improvement was achieved in ERA-Interim, but the upward transport remains faster than observations indicate (Dee et al., 2011). Interim results of JRA-55 suggest a less biased reanalysed temperature in the lower stratosphere relative to JRA-25 (Ebita et al., 2011).
In addition to the standard variables provided in reanalysis, i.e. air temperature, ozone mixing ratio and circulation characteristics -zonal, meridional or omega velocity -we have also analysed other dynamical variables. Of particular interest were the EP flux diagnostics -a theoretical framework to study interactions between planetary waves and the zonal mean flow (Andrews and McIntyre, 1987). Furthermore, this framework allows the study of the wave propagation characteristics in the zonal wind and the induced (large-scale) meridional circulation as well. For this purpose the quasi-geostrophic approximation of transformed Eulerian mean (TEM) equations were used in the form employed by Edmon et al. (1980), i.e. using their formula (3.1) for EP flux vectors, (3.2) for EP flux divergence and (3.4) for residual circulation. These variables were then interpolated to a regular vertical grid. For the visualisation purposes, the EP flux arrows were scaled by the inverse of the pressure. The script was publicly released (Kuchar, 2015). To detect variability and changes due to climate-forming factors, such as the 11-year SC, we have applied an attribution analysis based on multiple linear regression (MLR) and two nonlinear techniques. The regression model separates the effects of climate phenomena that are supposed to have an impact on middle atmospheric conditions. Our regression model of a particular variable X as a function of time t, pressure level p, latitude ϕ and longitude λ is described by the following equation: (1) After deseasonalising, which can be represented by the α index for every month in a year, the individual terms represent a trend regressor TREND(t) either in linear form or including the equivalent effective stratospheric chlorine (EESC) index (this should be employed due to the ozone turnover trend around the middle of the 90s), a SOLAR(t) represented by the 10.7 cm radio flux as a proxy for solar ultraviolet variations at wavelengths 200-300 nm that are important for ozone production and radiative heating in the stratosphere, and which correlates well with sunspot number variation (the data were acquired from Dominion Radio Astrophysical Observatory (DRAO) in Penticton, Canada).
We have also included the quasi-biennial proxies QBO 1,2,3 (t) as another stratosphere-related predictor. Similar studies have represented the QBO in multiple regression methods in several ways. Our approach involves three separate QBO indices extracted from each reanalysis. These three indices are the first three principal components of the residuals of our linear regression model (1) excluding QBO predictors applied to the equatorial zonal wind. The approach follows the paper by Frame and Gray (2010) or the study by Crooks and Gray (2005) to avoid contamination of the QBO regressors by the solar signal or other regressors. The three principal components explain 49, 47 and 3 % of the total variance for the MERRA; 60, 38 and 2 % for the JRA-55; and 59, 37 and 3 % for the ERA-Interim. The extraction of the first two components reveals a 28-month periodicity and an out-of phase relationship between the upper and lower stratospheres. The out-of phase relationship or orthogonality manifests approximately in a quarter period shift of these components. The deviation from the QBO quasi-regular period represented by the first two dominant components is contained in the residual variance. Linear regression analysis of the zonal wind with the inclusion of the first two principal components reveals a statistically significant linkage between the third principal component and the residuals of this analysis. Furthermore, the regression coefficient of this QBO proxy was statistically significant for all variables p value < 0.05 (see below for details about significance testing techniques). Wavelet analysis for the MERRA demonstrates three statistically significant but non-stationary periods exceeding the level of the white noise wavelet spectrum (not shown): an approximate annual cycle (a peak period of 1 year and 2 months), a cycle with a peak period of 3 years and 3 months and a long-period cycle (a peak period between 10 and 15 years). Those interferences can be attributed to the possible nonlinear interactions between the QBO itself and other signals like the annual cycle or long-period cycle such as the 11-year SC at the equatorial stratosphere.
The El Niño-Southern Oscillation is represented by the multivariate ENSO index ENSO(t) which is computed as the first principal component of the six main observed variables over the Pacific Ocean: sea level pressure, zonal and meridional wind, sea surface temperature, surface air temperature and total cloudiness fraction of the sky (NCAR, 2013). The effect of volcanic eruptions is represented by the stratospheric aerosol optical depth SAOD(t). The time series was derived from the optical extinction data (Sato et al., 1993). We have used globally averaged time series in our regression model. The North Atlantic Oscillation has also been included through its index NAO(t) derived by rotated principal component analysis applied to the monthly standardised 500 hPa height anomalies obtained from the Climate Data Assimilation System (CDAS) in the Atlantic region between 20 and 90 • N (NOAA, 2013).
The robustness of the solar regression coefficient has been tested in terms of including or excluding particular regressors in the regression model; e.g. the NAO term was removed from the model and the resulting solar regression coefficient was compared with the solar regression coefficient from the original regression set-up. The solar regression coefficient seems to be highly robust since neither the amplitude nor the statistical significance field was changed significantly when NAO or QBO 3 or both of them were removed. However, cross-correlation analysis reveals that the correlation between NAO and TREND, SOLAR and SAOD regressors is statistically significant, but small (not shown).
The multiple regression model of Eq.
(1) has been used for the attribution analysis, and supplemented by two nonlinear techniques. The MLR coefficients were estimated by the least squares method. To avoid the effect of autocorrelation of residuals and to obtain the best linear unbiased estimate (BLUE) according to the Gauss-Markov theorem (Thejll, 2005), we have used an iterative algorithm to model the residuals as a second-order autoregressive process. A Durbin-Watson test (Durbin and Watson, 1950) confirmed that the regression model was sufficient to account for most of the residual autocorrelations in the data.
As a result of the uncorrelated residuals, we can suppose the standard deviations of the estimated regression coefficients not to be diminished (Neter et al., 2004). The statistical significance of the regression coefficients was computed with a t test.
The nonlinear approach, in our case, consisted of a multilayer perceptron (MLP) and the relatively novel epsilon support vector regression (ε-SVR) technique with the threshold parameter ε = 0.1. The MLP as a technique inspired by the human brain is capable of capturing nonlinear interactions between inputs (regressors) and output (modelled data) (e.g. Haykin, 2009). The nonlinear approach is achieved by transferring the input signals through a sigmoid function in a particular neuron and within a hidden layer propagating to the output (a so-called feed-forward propagation). The standard error back-propagation iterative algorithm to minimise the global error has been used.
The support vector regression technique belongs to the category of kernel methods. Input variables were nonlinearly transformed to a high-dimensional space by a radial basis (Gaussian) kernel, where a linear classification (regression) can be constructed (Cortes and Vapnik, 1995). However, cross-validation must be used to establish a kernel parameter and cost function searched in the logarithmic grid from 10 −5 to 10 1 and from 10 −2 to 10 5 respectively. We have used 5fold cross-validation to optimise the SVR model selection for every point in the data set as a trade-off between the recommended number of folds (Kohavi, 1995) and computational time. The MLP model was validated by the holdout crossvalidation method since this method is more expensive in order of magnitude in terms of computational time. The data sets were separated into a training set (75 % of the whole data set) and a testing set (25 % of the whole data set). The neural network model was restricted to only one hidden layer with the maximum number of neurons set up to 20.
The earlier mentioned lack of explanatory power of the nonlinear techniques in terms of complicated interpretation of statistical models (Olden and Jackson, 2002) mainly comes from nonlinear interactions during signal propagation and the impossibility to directly monitor the influence of the input variables. In contrast to the linear regression approach, the understanding of relationships between variables is quite problematic. For this reason, the responses of our variables have been modelled by a technique originating from sensitivity analysis studies and also used by e.g. Blume and Matthes (2012). The relative impact RI of each variable was computed as where I k = σ (ŷ −ŷ k ). σ (ŷ −ŷ k ) is the variance of the difference between the original model outputŷ and the model outputŷ k when the k-input variable was held at its constant level. There are many possibilities with regard to which constant level to choose. It is possible to choose several levels and then to observe the sensitivity of model outputs varying for example on minimum, median and maximum levels. Our sensitivity measure (relative impact) was based on the median level. The primary reason comes from purely practical considerations -to compute our results fast enough as an-other weakness of the nonlinear techniques lies in the larger requirement of computational capacity. In general, this approach was chosen because of their relative simplicity for comparing all techniques to each other and to be able to interpret them too. The contribution of variables in neural network models has already been studied and Gevrey et al. (2003) produced a review and comparison of these methods. Figure 1a, d, g, j shows the annually averaged solar signal in the zonal means of temperature, zonal wind, geopotential height and ozone mixing ratio. The signal is expressed as the average difference between the solar maxima and minima in the period 1979-2013, i.e. normalised by 126.6 solar radio flux units. Statistically significant responses detected by the linear regression in the temperature series (see Fig. 1a) are positive and are located around the Equator in the lower stratosphere with values of about 0.5 K. The temperature response increases to 1 K in the upper stratosphere at the Equator and up to 2 K at the poles. The significant solar signal anomalies are more variable around the stratopause and not limited to the equatorial regions. Hemispheric asymmetry of the statistical significance can be observed in the lower mesosphere. From a relative impact point of view (in Fig. 2a-c marked as RI), it is difficult to detect a signal with an impact larger than 20 % in the lower stratosphere where the volcanic and QBO impacts dominate. In the upper layers (where the solar signal expressed by the regression coefficient is continuous across the Equator) we have detected relatively isolated signals (over 20 %) around ±15 • using the relative impact method. The hemispheric asymmetry also manifests in the relative impact field, especially in the SVR field in the mesosphere.

Annual response (MERRA)
The annually averaged solar signal in the zonal mean of zonal wind (Figs. 1d and 2d-f) dominates around the stratopause as an enhanced subtropical westerly jet. The zonal wind variability due to the SC corresponds to the temperature variability due to the change in the meridional temperature gradient and via the thermal wind equation. The largest positive anomaly in the Northern Hemisphere reaches 4 m s −1 around 60 km (Fig. 1d). In the Southern Hemisphere, the anomaly is smaller and not statistically significant. There is a significant negative signal in the southern polar region. The negative anomalies correspond to a weakening of the westerlies or an amplification of the easterlies. The relative impact of the SC is similarly located zonally even for both nonlinear techniques (Fig. 2d-f). The equatorial region across all the stratospheric layers is dominantly influenced by the QBO (expressed by all three QBO regressors) and for this reason the solar impact is minimised around the Equator. The pattern of the solar response in geopotential height  shows positive values in the upper stratosphere and lower mesosphere. This is also consistent with the zonal wind field through thermal wind balance. In the geopotential field, the SC influences the most extensive area among all regressors. The impact area includes almost the whole mesosphere and the upper stratosphere. Figure 1j also shows the annual mean solar signal in the zonal mean of the ozone mixing ratio (expressed as a percent change per annual mean). By including an EESC regressor term in the regression model instead of a linear trend over the whole period (for a detailed description see the methodology Sect. 3), we tried to capture the ozone trend change around the year 1996. Another possibility was to use our model over two individual periods, e.g. 1979-1995 and 1996-2013, but the results were quantitatively similar. The main common feature of the MERRA solar ozone response in Fig. 1j with observational results is the positive ozone response in the lower stratosphere, ranging from a 1 to 3 percent change.
In the equatorial upper stratosphere, no solar signal was detected that is comparable to that estimated from satellite measurement (Soukharev and Hood, 2006). By the relative impact method (Fig. 2j-l), we have obtained results comparable with linear regression coefficients, but especially around Atmos. Chem. Phys., 15, 6879-6895

Annual response -comparison with JRA-55, ERA-Interim
Comparison of the results for the MERRA, ERA-Interim and JRA-55 temperature, zonal wind and geopotential height shows that the annual responses to the solar signal are in qualitative agreement (compare individual plots in Fig. 1). The zonal wind and geopotential response seem to be consistent in all presented methods and data sets. The largest discrepancies can be seen in the upper stratosphere and especially in the temperature field (the first row in these figures). The upper stratospheric equatorial anomaly was not detected by any of the regression techniques in the case of the JRA-55 reanalysis although the JRA-25 showed a statistically significant signal with structure and amplitude of 1-1.25 K comparable with ERA-Interim in the equatorial stratopause (Mitchell et al., 2014). Although the anomaly in the MERRA temperature in Fig. 1a in the upper stratosphere is comparable to that in the ERA Interim temperature in Fig. 1b, the former signal is situated lower down at around 4 hPa (see also Mitchell et al., 2014). However, upper stratospheric temperature response could be less than accurate due to the existence of discontinuities in 1979, 1985and 1998(McLandress et al., 2014 coinciding with major changes in instrumentation or analysis procedure. Therefore, the temperature response to solar variation may be influenced by these discontinuities in the upper stratosphere. The revised analysis with the adjustments of ERA Interim temperature data from McLandress et al. (2014) showed in comparison with the original analysis without any adjustment that the most pronounced differences are apparent in higher latitudes and especially in 1 hPa. The regression coefficients decreased by about 50 % when using the adjusted data set, but the differences are not statistically significant in terms of 95 % confidence interval. The difference in tropical latitudes is about 0.2 K/(S max − S min ). The trend regressor t from Eq. 1 reveals a large turnaround from positive trend to negative in the adjusted levels, i.e. 1, 2, 3 and 5 hPa. Other regressors do not reveal any remarkable difference. The results in Figs. 1b, e, h, k and 3 from the raw data set were kept in order to refer and discuss the accordance and differences between our results and results from Mitchell et al. (2014), where no adjustments have been considered either.
The variability of the solar signal in the MERRA stratospheric ozone series was compared with the ERA-Interim results. The analysis points to large differences in the ozone response to the SC between the reanalyses and in comparison with satellite measurements by Soukharev and Hood (2006). In comparison with the satellite measurements, no relevant solar signal was detected in the upper stratosphere in the MERRA series. The signal seems to be shifted above Figure 3. The annually averaged response of the solar signal in the ERA-Interim zonal-mean temperature t (a-c), unit: K; zonal wind u (d-f), unit: m s −1 ; geopotential height h (g-i), unit: gpm; and ozone mixing ratio o3 (j-l), unit: percentage change per annual mean. The response is expressed as a relative impact RI approach. The relative impact was modelled by MLR, SVR and MLP techniques. The black contour levels in the RI plots are 0.2, 0.4, 0.8 and 1.0. the stratopause (confirmed by all techniques, shown in Figs. 2 and 3j-l). Regarding the ERA-Interim, there is a statistically significant ozone response to the SC in the upper stratosphere, but it is negative in sign, with values reaching up to 2 % above the Equator and up to 5 % in the polar regions of both hemispheres. However, a negative ozone and a positive temperature response in the upper stratosphere to a positive UV flux change from solar minimum to maximum is not physically reasonable. It must reflect an artifact of the assimilation model scheme and/or internal variability of the model rather than an effect of solar forcing (for more details about ozone as a prognostic variable in ERA-Interim, see Dee et al., 2011). There is a clear inverse correlation between the ERA-Interim temperature response in Fig. 1b and the ozone response in Fig. 1k. This does probably imply that the temperature response is producing the negative ozone response in the assimilation model. However, it is not physically reasonable because both the ozone and the temperature in the upper stratosphere respond positively to an increase in solar UV (e.g. Hood et al., 2015). In the case of MERRA, while SBUV ozone profiles are assimilated with SC passed to the forecast model (as the ozone analysis tendency contribution), no SC was passed to the radiative part of the model. The same is also true for ERA-Interim and JRA-55 (see the descriptive table of the reanalysis product on SC in irradiance and ozone in Mitchell et al. (2014). Despite the fact that the analysed ozone should contain a solar signal, the signal is not physically reasonable and is dominated by internal model variability in terms of dynamics and chemistry. Since the SBUV ozone profiles have very low vertical resolution, this may also affect the ozone response to the SC in the MERRA reanalysis. These facts should also be taken into account in case of monthly response discussion of particular variables in Sect. 4.2.
The lower stratospheric ozone response in the ERAinterim is not limited to the equatorial belt ±30 • up to 20 hPa, as in the case of the MERRA reanalysis, and the statistical significance of this signal is rather reduced. The solar signal is detected higher and extends from the subtropical areas to the polar regions. The results suggest that the solar response in the MERRA series is more similar to the results from satellite measurements (Soukharev and Hood, 2006). Nevertheless, further comparison with independent data sets is needed to assess the data quality in detail.

Comparison of the linear and nonlinear approaches (MLR vs. SVR and MLP)
In this paper, we have applied and compared one linear (MLR) and two nonlinear attribution (SVR and MLP) techniques. The response of the studied variables to the solar signal and other forcings was studied using the sensitivity analysis approach in terms of averaged response deviation from the equilibrium represented by the original model outputŷ ( Blume and Matthes, 2012). This approach does not recognise a positive or negative response as the linear regression does. For this reason, the relative impact results are compared to the regression's coefficients. Using linear regression, it would be possible to assess the statistical significance of the regression's coefficients and a particular level of the relative impact since they are linearly proportional. A comparison between the linear and nonlinear approaches by the relative impact fields shows qualitative and in most regions also quantitative agreement. The most pronounced agreement is observed in the zonal wind (Figs. 2, 3 and 4d-f) and geopotential height fields (Figs. 2, 3 and 4g-i). On the other hand worse agreement is captured in the ozone and temperature field. In the temperature field the upper stratospheric solar signal reaches values over 20 %, some individual signals in the Southern Hemisphere even reach 40 %. However, using the relative impact approach, the lower stratospheric solar signal in the temperature field (which is well established by the regression coefficient) does not even reach 20 % because of the dominance of the QBO and volcanic effects. These facts emphasise that nonlinear techniques contribute to the robustness of attribution analysis since the linear regression results were plausibly confirmed by the SVR and MLP techniques.
In conclusion, the comparison of various statistical approaches (MLR, SVR and MLP) should actually contribute to the robustness of the attribution analysis including the statistically assessed uncertainties. These uncertainties could partially stem from the fact that the SVR and neural network techniques are dependent on an optimal model setting which is based on a rigorous cross-validation process, which places a high demand on computing time.
The major differences between the techniques can be seen in how much of the temporal variability of the original time series is explained, i.e. in the coefficient of determination. For instance, the differences of the explained variance reach up to 10 % between linear and nonlinear techniques, although the zonal structure of the coefficient of determination is almost the same. To conclude, nonlinear techniques show an ability to simulate the middle atmosphere variability with a higher accuracy than cross-validated linear regression.

Monthly response (MERRA)
As was pointed out by Frame and Gray (2010), it is necessary to examine the solar signal in individual months because of a solar impact on polar-night jet oscillation (Kuroda and Kodera, 2001). For example, the amplitude of the lower stratospheric solar signal in the northern polar latitudes in February exceeds the annual response since the SC influence on vortex stability is most pronounced in February. Besides the radiative influences of the SC, we discuss the dynamical response throughout the polar winter (Kodera and Kuroda, 2002). ±0.25, ±0.5, ±1, ±2, ±5, ±10, ±15, ±30; zonal wind u (e-h), unit: m s −1 , contour levels: 0, ±1, ±2, ±5, ±10, ±15, ±30; geopotential height h (j-l), unit: gpm, contour levels: 0, ±10, ±20, ±50, ±100, ±150, ±300; EP flux divergence EPfD (m-p), unit: m s −1 day −1 ; together with EP flux vectors scaled by the inverse of the pressure, unit: m 2 s −2 ; and ozone mixing ratio, unit: percentage change per monthly mean; with residual circulation o3 + rc (q-t), units: m s −1 ; −10 −3 Pa s −1 during northern hemispheric winter. The response is expressed as a regression coefficient (corresponding units per S max minus S min ). The statistical significance of the scalar fields was computed by a t test. Red and yellow areas in Panels (a-l) and grey contours in Panels (m-t) indicate p values of < 0.05 and 0.01 respectively. Statistically significant upper stratospheric equatorial anomalies in the temperature series (winter months in Figs. 5 and 6a-d) are expressed in almost all months. Their amplitude and statistical significance vary throughout the year. The variation between the solar maxima and minima could be up to 1 K in some months. Outside the equatorial regions, the fluctuation could reach several Kelvin. The lower stratospheric equatorial anomaly strengthens during winter. This could be an indication of dynamical changes, i.e. alteration of the residual circulation between the equatorial and polar regions (for details, please see Sect. 5). Aside from the radiative forcing by direct or ozone heating, other factors are linked to the anomalies in the upper levels of the mid-dle atmosphere (Haigh, 1994;Gray et al., 2009). It is necessary to take into consideration the dynamical coupling with the mesosphere through changes of the residual circulation (see the dynamical effects discussion below). That can be illustrated by the positive anomaly around the stratopause in February (up to 4 K around 0.5 hPa). This anomaly extends further down and, together with spring radiative forcing, affects the stability of the equatorial stratopause. Hemispheric asymmetry in the temperature response above the stratopause probably originates from the hemispheric differences, i.e. different wave activity (Kuroda and Kodera, 2001). These statistically significant and positive temperature anomalies across the subtropical stratopause begin to descend and move Figure 6. The monthly averaged response of the solar signal in the MERRA zonal-mean temperature t (a-d); unit: K; contour levels: 0, ±0.25, ±0.5, ±1, ±2, ±5, ±10, ±15, ±30; zonal wind u (e-h), unit: m s −1 ; contour levels: 0, ±1, ±2, ±5, ±10, ±15, ±30; geopotential height h (j-l); unit: gpm; contour levels: 0, ±10, ±20, ±50, ±100, ±150, ±300; EP flux divergence EPfD (m-p), unit: m s −1 day −1 ; together with EP flux vectors scaled by the inverse of the pressure; unit: m 2 s −2 ; and ozone mixing ratio, unit: percentage change per monthly mean; with residual circulation o3+ rc (q-t); units: m s −1 , −10 −3 Pa s −1 during southern hemispheric winter. The response is expressed as a regression coefficient (corresponding units per S max minus S min ). The statistical significance of the scalar fields was computed by a t test. Red and yellow areas in Panels (a-l) and grey contours in Panels (m-t) indicate p values of < 0.05 and 0.01 respectively. to higher latitudes in the beginning of the northern winter. The anomalies manifest fully in February in the region between 60 and 90 • N and reach tropospheric levels -contrary to the results for the Southern Hemisphere (see Fig. 10 in Mitchell et al., 2014). The southern hemispheric temperature anomaly is persistent above the stratopause and the SC influence on the vortex stability differs from those in the Northern Hemisphere.
The above described monthly anomalies of temperature correspond to the zonal wind anomalies throughout the year (Figs. 5 and 6e-h). The strengthening of the subtropical jets around the stratopause is most apparent during the winter in both hemispheres. This positive zonal wind anomaly gradu-ally descends and moves poleward, similar to the Frame and Gray (2010) analysis based on ERA-40 data. In February, the intensive stratospheric warming and mesospheric cooling is associated with a more pronounced transition from winter to summer circulation attributed to the SC (in relative impact methodology up to 30 %). However, GCMs have not yet successfully simulated the strong polar warming in February (e.g. Schmidt et al., 2010;Mitchell et al., 2015b). Due to the short (35-year) time series, it is possible that this pattern is not really solar in origin but is instead a consequence of internal climate variability or aliasing from the effects of the two major volcanic eruptions aligned to solar maximum periods.
In the Southern Hemisphere, this poleward motion of the positive zonal wind anomaly halts approximately at 60 • S. For example, in August, we can observe a well-marked latitudinal zonal wind gradient (Fig. 6h). Positive anomalies in the geopotential height field correspond to the easterly zonal wind anomalies. The polar circulation reversal is associated with intrusion of ozone from the lower latitudes, as is apparent e.g. in August in the Southern Hemisphere and in February in the Northern Hemisphere (last rows of Figs. 5 and 6).
When comparing the results from the MERRA and ERA-40 series studied by Frame and Gray (2010), distinct differences were found (Fig. 5e, f) in the equatorial region of the lower mesosphere in October and November. While in the MERRA reanalysis we have detected an easterly anomaly above 1 hPa in both months (only November shown), a westerly anomaly was identified in the ERA-40 series. Further distinct differences in the zonal mean temperature and zonal wind anomalies were not found.

Dynamical effects discussion
In this section, we discuss the dynamical impact of the SC and its influence on middle atmospheric winter conditions. Linear regression was applied to the EP diagnostics. Kodera and Kuroda (2002) suggested that the solar signal produced in the upper stratosphere region is transmitted to the lower stratosphere through the modulation of the internal mode of variation in the polar-night jet and through a change in the Brewer-Dobson circulation (prominent in the equatorial region in the lower stratosphere). In our analysis, we discussed the evolution of the winter circulation with an emphasis on the vortex itself rather than the behaviour of the jets. Furthermore, we try to describe the possible processes leading to the observed differences in the quantities of state between the solar maximum and minimum period. Because the superposition principle only holds for linear processes, it is impossible to deduce the dynamics merely from the fields of differences. As noted by Kodera and Kuroda (2002), the dynamical response of the winter stratosphere includes highly nonlinear processes, e.g. wave-mean flow interactions. Thus, both the anomaly and the total fields, including climatology, must be taken into account.
We start the analysis of solar maximum dynamics with the period of the northern hemispheric winter circulation formation. The anomalies of the ozone, temperature, geopotential in the lower stratosphere only and Eliassen-Palm flux divergence mostly in the upper stratosphere support the hypothesis of weaker BDC during the solar maximum due to the less intensive wave pumping. This is possible through the "downward control" principle when modification of wavemean flow interaction in the upper levels governs changes in residual circulation below (Haynes et al., 1991). The finding about weaker BDC during the solar maximum is consistent with previous studies (Kodera and Kuroda, 2002;Matthes et al., 2006). The causality is unclear, but the effect is visible in both branches of BDC as is illustrated by Fig. 5 and summarised schematically in Fig. 7.
During the early Northern Hemisphere (NH) winter (including November) when westerlies develop in the stratosphere, we can observe a deeper polar vortex and consequent stronger westerly winds both inside and outside the vortex. However, only the westerly anomaly outside the polar region and around 30 • N from 10 hPa to the lower mesosphere is statistically significant (see the evolution of zonal wind anomalies in Fig. 5e-h). The slightly different wind field has a direct influence on the vertical propagation of planetary waves. From the Eliassen-Palm flux anomalies and climatology we can see that the waves propagate vertically with increasing poleward instead of equatorward meridional direction with height. This is then reflected in the EP flux divergence field, where the region of maximal convergence is shifted poleward and the anomalous convergence region emerges inside the vortex above approximately 50 hPa ( Fig. 5m-p).
The poleward shift of the maximum convergence area further contributes to the reduced BDC. This is again confirmed by the temperature and ozone anomalies. The anomalous convergence inside the vortex induces anomalous residual circulation, the manifestation of which is clearly seen in the quadrupole-like temperature structure (positive and negative anomalies are depicted schematically in Fig. 7 using red and blue boxes respectively). This pattern emerges in November and even more clearly in December. In December, the induced residual circulation leads to an intrusion of the ozonerich air into the vortex at about the 1 hPa level (Fig. 5r). The inhomogeneity in the vertical structure of the vortex is then also pronounced in the geopotential height differences. This corresponds to the temperature analysis in the sense that above and in the region of the colder anomaly there is a negative geopotential anomaly and vice versa. The geopotential height difference has a direct influence on the zonal wind field (via the thermal wind balance). The result is a deceleration of the upper vortex parts and consequent broadening of the upper parts (due to the conservation of angular momentum).
Considering the zonal wind field, the vortex enters January approximately with its average climatological extent. The wind speeds in its upper parts are slightly higher. This is because of the smaller geopotential values corresponding to the negative temperature anomalies above approximately 1 hPa. This probably results from the absence of adiabatic heating due to the suppressed BDC, although the differences in the quantities of state (temperature and geopotential height) are small and insignificant (see the temperature anomalies in Fig. 5c). It is important to note that these differences change sign around an altitude of 40 km inside the vortex further accentuating the vertical inhomogeneity of the vortex. This might start balancing processes inside the vortex, which is confirmed by analysis of the dynamical quantities, i.e. EP flux and its divergence (Fig. 5o) Figure 7. Solar cycle modulation of the winter circulation: schema of the related mechanisms. The upper and lower figure show early and later winter respectively. The heating and cooling anomalies are drawn with red and blue boxes. The EP flux divergence and convergence are drawn with green and yellow boxes. The wave propagation anomaly is expressed as a wavy red arrow in contrast to the climatological average drawn by a wavy grey arrow. The induced residual circulation according to the quasi-geostrophic approximation is highlighted by the bold black lines.
Significant anomalies of the EP flux indicate anomalous vertical wave propagation resulting in the strong anomalous EP flux convergence being significantly pronounced in a horizontally broad region and confined to upper levels (convergence (negative values) drawn by green or blue shades in Fig. 5m-p). This leads to the induction of an anomalous residual circulation starting to gain intensity in January. The situation then results in the disruption of the polar vortex visible in significant anomalies in the quantities of state in February -in contrast to January. Further strong mixing of air is suggested by the ozone fields. The quadrupole-like structure of the temperature is visible across the whole NH middle atmosphere in February (indicated in the lower diagram of Fig. 7), especially in the higher latitudes. This is very significant and well pronounced by the stratospheric warming and mesospheric cooling.
The hemispheric asymmetry of the SC influence can be especially documented in winter conditions, as was already suggested in Sect. 4.2. Since the positive zonal wind anomaly halts at approximately 60 • S and intensifies over 10 m s −1 , one would expect the poleward deflection of the planetary wave propagation to be according to NH winter mechanisms discussed above. This is actually observed from June to August when the highest negative anomalies of the latitudinal component of EP flux are located in the upper stratosphere and in the lower mesosphere ( Fig. 6m-p). The anomalous divergence of EP flux develops around the stratopause between 30 and 60 • S. Like the hypothetical mechanism of weaker BDC described above, we can observe less wave pumping in the stratosphere and consequently less upwelling in the equatorial region. In line with that, we can see in the lower stratosphere of equatorial region (Figs. 5b and 6b) a more pronounced temperature response in August (above 1 K) than in December (around 0.5 K) as already mentioned in previous observational (van Loon and Labitzke, 2000) or reanalysis (Mitchell et al., 2014) studies. Although this can point to a weaker BDC, the residual circulation ( Fig. 6q-t) as a proxy for BDC (Butchart, 2014) does not reveal this signature. Hypothetically, this could be due to a higher role of unresolved wave processes in reanalysis (small-scale gravity waves) or due to the worse performance of residual circulation as a proxy for the large-scale transport in SH (e.g. larger departure from steady waves approximation comparing to NH), or because of the other processes than BDC leading to the temperature anomaly, e.g. aliasing with volcanic signal.
Overall, the lower stratospheric temperature anomaly is more coherent for the SH winter than for the NH winter, where the solar signal is not so apparent or statistically significant in particular months and reanalysis data sets.

Conclusions
We have analysed the changes in air temperature, ozone and circulation characteristics driven by the variability of the 11year solar cycle's influence on the stratosphere and lower mesosphere. Attribution analysis was performed on the three reanalysed data sets, MERRA, ERA-Interim and JRA-55, and aimed to compare how these types of data sets resolve the solar variability throughout the levels where the "top-down" mechanism is assumed. Furthermore, the results that originated in linear attribution using MLR were compared with other relevant attribution studies and supported by nonlinear attribution analysis using SVR and MLP techniques.
The nonlinear approach to attribution analysis, represented by the application of the SVR and MLP, largely confirmed the solar response computed by linear regression. Consewww.atmos-chem-phys.net/15/6879/2015/ quently, these results can be considered quite robust regarding the statistical modelling of the solar variability in the middle atmosphere. This finding indicates that linear regression is a sufficient technique to resolve the basic shape of the solar signal through the middle atmosphere. However, some uncertainties could partially stem from the fact that the SVR and MLP techniques are highly dependent on an optimal model setting that requires a rigorous cross-validation process (which places a high demand on computing time). As a benefit, nonlinear techniques show an ability to simulate the middle atmosphere variability with higher accuracy than linear regression.
The solar signal extracted from the temperature field from MERRA and ERA-Interim reanalysis using linear regression has the amplitudes around 1 and 0.5 K, in the upper stratospheric and in the lower stratospheric equatorial region respectively. However, the peak amplitudes of the temperature response in the equatorial upper stratosphere occur at different levels (about 4 and 2 hPa respectively). These signals, statistically significant at a p value < 0.01, are in qualitative agreement with previous attribution studies (e.g. Frame and Gray, 2010;Mitchell et al., 2014). A statistically significant signal was only observed in the lower part of the stratosphere in the JRA-55 reanalysis, however with similar amplitudes as the other data sets.
Similar to the temperature response, the double-peaked solar response in ozone was detected in satellite measurements (e.g. Soukharev and Hood, 2006), although concerns were expressed about the physical mechanism of the lower stratospheric response (e.g. Austin et al., 2008). However, the exact position and amplitude of both ozone anomalies remain a point of disagreement between models and observations. The results of our attribution analysis point to large differences in the upper stratospheric ozone response to the SC in comparison with the studies mentioned above and even between reanalyses themselves. The upper stratospheric ozone anomaly reaches 2 % in the SBUV(/2) satellite measurements (e.g. Soukharev and Hood, 2006, Fig. 5) which were assimilated as the only source of ozone profiles in MERRA reanalysis. This fact is remarkable since the same signal was not detected in the upper stratosphere in the MERRA results. However, the solar signal in the ozone field seems to be shifted above the stratopause where similar and statistically significant solar variability was attributed. Concerning the solar signal in the ERA-Interim, there is a negative ozone response via a regression coefficient in the upper stratosphere, although the solar variability expressed as relative impact appears to be in agreement with satellite measurements. The negative ozone response in the tropical upper stratosphere is not consistent with physical expectations for a nominal positive change in solar UV irradiance (e.g. Hood et al., 2015).
Furthermore, the lower stratospheric solar response in the ERA-Interim's ozone around the Equator is reduced in this data set and shifted to higher latitudes. Another difference was detected in the monthly response of the zonal wind in October and November in the equatorial region of the lower mesosphere between the results for the MERRA series and ERA-40 data studied by Frame and Gray (2010). While in the MERRA reanalysis we have detected an easterly anomaly, a westerly anomaly was identified in the ERA-40 series.
A similar problem with the correct resolving of the doublepeaked ozone anomaly was registered in the study of Dhomse et al. (2011) which investigated the solar response in the tropical stratospheric ozone using a 3-D chemical transport model. The upper stratospheric solar signal observed in SBUV/SAGE and SAGE-based data could only be reproduced in model runs with unrealistic dynamics, i.e. with no inter-annual meteorological changes.
The reanalyses have proven to be extremely valuable scientific tools (Rienecker et al., 2011). On the other hand, they have to be used with caution, for example, due to the existence of large discontinuities occurring in 1979, 1985and 1998(McLandress et al., 2014) that translated into errors in the derived solar coefficients. Our revised analysis with the adjustments from McLandress et al. (2014) resulted in an 0.2 K/(S max − S min ) reduction in the temperature solar regression coefficients in tropical latitudes of the upper stratosphere.
In the dynamical effects discussion, we described the dynamical impact of the SC on middle atmospheric winter conditions. The relevant dynamical effects are summarised in schematic diagrams (Fig. 7). Both diagrams depict average conditions and anomalies induced by the SC. The first one summarises how equatorward wave propagation is influenced by the westerly anomaly around the subtropical stratopause. The quadrupole-like temperature structure is explained by anomalous residual circulation in the higher latitudes together with the anomalous branch heading towards the equatorial region already hypothesised by Kodera and Kuroda (2002). The second diagram concludes the transition time to vortex disruption during February. Again, a very apparent quadrupole-like temperature structure is even more pronounced, especially in the polar region, and seems to be more extended to lower latitudes.
Fields of residual circulation and EP flux divergence in February are opposite to what would be expected from the suppressed BDC in the SC max. There is an enhanced downwelling in the polar and an enhanced upwelling in the equatorial region below 1 hPa. This suggests a need to diagnose the influence of SC on transport at least on a monthly scale because the changes in the underlying dynamics (compare the upper and lower diagrams in Fig. 7) would make the transport pathways more complicated. The negative zonal wind response in late northern winter may be caused by an increased likelihood of major stratospheric warmings later in the winter under solar maximum conditions when the polar vortex in early winter is stronger, on average, and therefore less susceptible to disruption (e.g. Gray et al., 2004). Since GCMs have not yet successfully simulated this pattern (e.g. Schmidt et al., 2010;Mitchell et al., 2015b) and due to the short (35-year) time series, it is possible that this pattern is not really solar in origin but is instead a consequence of internal climate variability or aliasing from effects of the two major volcanic eruptions aligned to solar maximum periods.
However, we can strongly assume that the dynamical effects are not zonally uniform, as is shown here using two-dimensional (2-D) EP diagnostics and TEM equations. Hence, it would be interesting to extend the discussion of dynamical effects for other relevant characteristics, for example, for the analysis of wave propagation and wave-mean flow interaction using the 3-D formulation (Kinoshita and Sato, 2013).
This paper is fully focused on the SC influence, i.e. on decadal changes in the stratosphere and lower mesosphere, although a huge number of results concerning other forcings was generated by attribution analysis. The QBO phenomenon in particular could be one of the points of future interest since the solar-QBO interaction and the modulation of the Holton-Tan relationship by the SC are regarded as highly challenging, especially in global climate simulations .