Systematic attribution of observed Southern Hemisphere circulation trends to external forcing and internal variability

A critical question in the global warming debate concerns the causes of the observed trends of the Southern Hemisphere (SH) atmospheric circulation over recent decades. Secular trends have been identified in the frequency of occurrence of circulation regimes, namely the positive phase of the Southern Annular Mode (SAM) and the hemispheric wave-3 pattern which is associated with blocking. Previous studies into the causes of these secular trends have either been purely model based, have not included observational forcing data or have mixed external forcing with indices of internal climate variability impeding a systematic and unbiased attribution of the causes of the secular trends. Most model studies also focused mainly on the austral summer season. However, the changes to the storm tracks have occurred in all seasons and particularly in the austral winter and early spring when midlatitude blocking is most active and stratospheric ozone should not play a role. Here we systematically attribute the secular trends over the recent decades using a non-stationary clustering method applied to both reanalysis and observational forcing data from all seasons. While most previous studies emphasized the importance of stratospheric ozone depletion in causing austral summer SH circulation trends, we show observational evidence that anthropogenic greenhouse gas concentrations have been the major driver of these secular trends in the SAM and blocking when all seasons are considered. Our results suggest that the recovery of the ozone hole might delay the signal of global warming less strongly than previously thought and that effects from all seasons are likely crucial in understanding the causes of the secular trends.


Introduction
The SH climate and atmospheric circulation has undergone significant changes over the last few decades.It is important to understand its causes and anthropogenic contributions because this will not only help to constrain future climate projections but is also essential to evaluate the ability of the current generation of climate models to accurately simulate these changes.Previous attempts at understanding the causes of SH climate change focused on changes in the mean climate and its variance (Arblaster and Meehl, 2006;Turner et al., 2009;Thompson et al., 2011) during austral summer.
Here we instead focus on changes in the frequency of occurrence of circulation regimes (O'Kane et al., 2013;Lee and Feldstein, 2013) in all seasons.Such circulation regimes have a significant impact on surface weather and climate.For instance, the frequency of occurrence of blocking strongly affects temperature and precipitation (Risbey et al., 2009;Pook et al., 2013).Furthermore, the secular trend towards the increased occurrence of the positive SAM phase, linked in the austral summer to stratospheric ozone (Son et al., 2010;Polvani et al., 2011;Thompson et al., 2011;Previdi and Polvani, 2014;Barnes et al., 2013), affects Antarctic temperatures and sea ice extent (Turner et al., 2009).It also has to be noted that the way SAM is defined can impact on the outcomes of any study using a SAM index (Ho et al., 2012).Hence, we use here a method which does not presuppose any particular spatial structure on the resultant regimes and the particular form of the SAM we find is determined from the data.
Over the recent decades (1980-present) large changes in the SH storm track modes have occurred in all seasons including the austral winter, when blocking is at its most active (Marques and Rao, 2000).The austral winter storm track changes manifest as reduced baroclinicity and a decrease in the July zonal winds of about ≈ 10 m s −1 in the subtropical jet relative to the earlier period   (Frederiksen and Frederiksen, 2007).The changes in the storm track modes provide a dynamical mechanism for the observed systematic linear downward trends in the annual number of SH blocking events (Weidenmann et al., 2002;O'Kane et al., 2013).Such events predominantly occur in preferred locations about the Australian (110-210 • E), eastern Pacific (260-315 • E) and Indian (20-80 • E) ocean sectors.These regions are associated with the ridges of the hemispheric wave-3 pattern (van Loon and Jenne, 2002;Trenberth and Mo, 1985).
A recent study (O'Kane et al., 2013) using non-stationary clustering has shown that, consistent with reduced blocking activity post the late 1970s, the wave-3 pattern has weakened while the corresponding zonal state (SAM) has strengthened and moved poleward (positive phase).While an early study hypothesized that global warming is likely to change the frequency of occurrence of these circulation regimes but not their spatial patterns (Corti et al., 1999;O'Kane et al., 2013).Furthermore, O' Kane et al. (2013) has shown that the spatial character of the SH persistent climate regimes has changed significantly over the reanalysis period 1948-2009.Here we extend our earlier study of SH circulation regimes to show observational evidence that when all seasons are considered these changes are mainly in response to radiative forcing trends of anthropogenic CO 2 emissions and to a much lesser degree to stratospheric ozone depletion.This is likely due to the fact that ozone depletion plays a minor role outside of the austral summer season.
The main method of investigating the role of different forcings of the SH secular trends has been through coupled climate models (Stocker et al., 2013).Modeling studies (Arblaster and Meehl, 2006;Turner et al., 2009;Son et al., 2010;Previdi and Polvani, 2014;Barnes et al., 2013) have identified stratospheric ozone depletion as an important driver of the observed austral summertime intensification of the SAM over the recent decades.However, stratospheric ozone depletion is a highly seasonal effect and can play no role in the austral winter-spring atmospheric circulation dynamics.The SH storm tracks are equally active all year around (Trenberth, 1991).However, the austral wintertime is the season when the observed changes to the storm track activity, namely reduced blocking and baroclinicity of the subtropical jet (Frederiksen and Frederiksen, 2007), have been particularly evident and cannot be solely attributed to the ozone mass deficit (OMD).
Modeling studies of the effect of various individual and combined radiative forcings on the SH circulation have largely compared trends in mean zonal indices, mainly the SAM index, without consideration of related systematic changes in the spatially coherent zonally asymmetric features of the circulation (O' Kane et al., 2013).Such studies often rely on ensemble averaging of repeated forcing experiments to enhance the forcing signal while simultaneously reducing intrinsic interannual-to decadal-scale variability.The problem with this approach is the prohibitive computational cost of coupled models, allowing only for relatively small ensemble sizes whose models often have rather coarse resolutions and contain many biases.
CMIP5 (Coupled Model Intercomparison Project Phase 5) models are known to poorly represent midlatitude blocking and to be limited in their ability to capture important SH circulation responses such as the response of the SAM to large volcanic eruptions (Karpechko, 2010;Gleixner, 2012;Charlton-Perez et al., 2013).Furthermore, recent studies showed that very high horizontal resolutions (grid spacing of about 16 km) are necessary for climate models to accurately simulate the geographical structure and probability distribution of blocking and regional weather regimes (Dawson et al., 2012;Dawson and Palmer, 2014).
Given that the circulation changes are a key signature of the forcing for attribution, the main aim of this paper is to complement model-based results with observationally based studies and to try to separate natural variability from the forced response.For example, the late 1970s climate shift occurred coincident with the shift in phase of the Interdecadal Pacific Oscillation (IPO), so separating the low frequency intrinsic ENSO (El Niño-Southern Oscillation) behavior from a response to the constituent components of the radiative forcing is an important problem.
Here we argue that most studies (see Thompson et al. (2011) for a review) attributing secular trends in the SH circulation have exclusively focused on changes in the mean zonal circulation and the trend toward positive SAM in the summer months.This has been largely attributed to ozone because model simulations with (without) ozone can (cannot) reproduce the magnitude of the trend and because there is only a small trend towards the positive SAM in the winter; thus, the conclusion is that ozone is the major driver of the secular SAM trend.The mechanism proposed (Son et al., 2010) is that stratospheric cooling/heating allows the high latitude tropopause to rise in the summer enabling a poleward movement of the westerlies and consequently the Hadley cell.Polvani et al. (2011) found that the only statistically significant relationship between the position of the Hadley cell and the midlatitude jet exists in the summer; however, in the winter no such relationship exists even though the wintertime is when one expects such a relationship to be most robust.Here we show results of a statistical analysis using zonally asymmetric fields which not only represent the SAM pattern but also zonal asymmetries, i.e., blocking, which has also undergone secular changes over the last few decades.
In Sect. 2 we present the data used in this study and in Sect. 3 we describe the statistical method used for the nonstationary clustering.In Sect. 4 we present the attribution re-Nonlin.Processes Geophys., 22, 513-525, 2015 www.nonlin-processes-geophys.net/22/513/2015/ sults and also describe in detail our sensitivity tests regarding the number of parameters to be estimated and demonstrate the robustness of our results.We provide our conclusions in Sect. 5.

Data
We use daily NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) reanalysis data (Kalnay et al., 1996) covering the period 1980-2007 for 500 hPa geopotential height and surface air temperature.We consider only anomalies with respect to the climatological mean where the mean seasonal cycle has been removed but not detrended.Note that there is still an annual cycle in higher moments and in the frequency of occurrence present in the time series.While there are still large biases in the Antarctic region in the various reanalysis products (Bromwich and Fogt, 2004;Bromwich et al., 2007), we have shown in a similar study (O'Kane et al., 2015) that by using the Japanese 55-year Reanalysis (JRA-55) conducted by the Japan Meteorological Agency (JMA) we found very similar results.Hence, our results do not depend strongly on the used reanalysis data set.
As forcing data we use the Cape Grim CO 2 measurements (Steele et al., 2007), sulfate aerosols (Skeie et al., 2011), stratospheric aerosol optical thickness (Bourassa et al., 2010) (available at http://data.giss.nasa.gov/modelforce/strataer/),stratospheric ozone mass deficit (Roscoe and Haigh, 2007), and the solar constant (Fröhlich, 2000).Most of the forcing data is in monthly mean resolution.Since we are using daily reanalysis data for the clustering, we expand the monthly forcing data to daily resolution by using the monthly mean values for each day of the respective month.Because stratospheric ozone depletion has a strong annual cycle we carried out sensitivity analysis by lagging the ozone mass deficit values by 1, 2 or 3 months and we also used a 365-day backward running mean.The forcing time series are displayed in Fig. 1.
As internal modes of climate variability, we use an ENSO 3.4 index, the Madden-Julian Oscillation (MJO) index, the Indian Ocean Dipole (IOD) and the eastern Indian Ocean Dipole mode indices and the annual cycle here defined as sin(2 • π /365 • t).These indices describe tropical sea surface temperature (SST) variability (ENSO, IOD) or an intrinsic mode of tropical variability (MJO).We consider these to be intrinsic drivers of midlatitude variability but recognize that they likely also respond to changes in external radiative forcing such that a clear separation between cause and effect is difficult.However, it is still important to elucidate which role they play in the frequency changes of the regime patterns.We do not consider Antarctic sea ice extent because of its marginal expansion and because this slight expansion in extent is largely wind-driven (Holland and Kwok, 2012) and likely a response to the changes in the large-scale circulation.Furthermore, the changes in sea ice extent and area have been spatially heterogeneous, with increases in some areas like in the Ross Sea and decreases in other areas like in the Bellingshausen/Amundsen seas (Parkinson and Cavalieri, 2012).This is despite the trend towards the positive SAM; thus, it is unlikely that sea ice extent would have a significant impact on the secular circulation trends.

Non-stationary clustering
We first give an intuitive description of the clustering method used before we explain it in much more detail in Sect.3.1.That section can be skipped by readers who are more interested in the clustering results.
Many studies have provided evidence that the atmospheric circulation can be efficiently described by a few persistent cluster states (Cheng and Wallace, 1993;Kimoto and Ghil, 1993;Corti et al., 1999;Horenko et al., 2008;Majda et al., 2006;Franzke et al., 2008Franzke et al., , 2009Franzke et al., , 2011;;O'Kane et al., 2013;Risbey et al., 2015).Conventional clustering methods such as k-means partition phase space using heuristic algorithms, for example using empirical orthogonal functions (EOFs), into an a priori arbitrarily determined k set of cluster centroids whose points within each cluster are close but where each centroid is in some sense far apart from each other (Dawson and Palmer, 2014).Similarly, self organizing maps (SOMs; Johnson et al., 2008) are typically based on minimizing the geometric (Euclidean) distance between the observational data and some specified set of recurrent patterns but without consideration of the persistency of those states and mostly without considering the dynamics and differences in dynamics within these states (Michelangeli et al., 1995).Fur-thermore, classical clustering methods do not consider differences in the dynamics of the cluster states (Christiansen, 2007).Recently, Lee and Feldstein (2013) applied SOMs to reanalyzed daily zonal-mean zonal wind data for the austral summer period employing four SOM patterns with the DJF (December-January-February) global mean temperature taken as an indicator of the response to GHG (greenhouse gas) forcing.They attribute the main response to ozone by correlating the third SOM DJF zonal-mean zonal wind pattern with the November Antarctic ozone index.One of the central problems in applying these methods to historical geophysical data is related to the robustness, meaning that increasing the k away from k = 1 increases the total number of parameters -thereby increasing the risk of overfitting.Concepts from information theory like the Akaike information criterion (AIC) (Burnham and Anderson, 2002) are usually deployed to find the optimal number of clusters k that allows avoiding overfitting.
Our approach considers the hemispheric response of not only the zonal annular mode but also systematic changes in wave 3 and blocking.We also use reanalyzed data but consider all possible combinations of the observed radiative forcings and relevant indices of internal variability.Moreover, we make no a priori assumptions on the number of states and employ an approach that considers persistency, changes in the dynamics and that is able to ascertain causation between the time series and the external forcings.
Specifically, we use the non-stationary clustering method FEM-BV-VARX (finite element method of time series analysis with bounded variation and vector autoregressive factor of model parameters; Horenko, 2010;O'Kane et al., 2013;Risbey et al., 2015;Gerber and Horenko, 2014) to systematically attribute circulation trends to observed external forcings.In this study we assume that the large-scale circulation can be effectively decomposed into a small set of distinct patterns or regime states.On each day the atmosphere is only in one of these patterns, where it might stay for some time before it switches to one of the other regime states.If we order the regime states from 1 to n, we get a daily index denoting in which state the atmosphere is in on that particular day; this index is referred to in the literature as the Viterbi path (Viterbi, 1967;Franzke et al., 2008).We then construct a statistical model which simultaneously estimates the geographical structure of these patterns and the evolution of the switching between the patterns.We do this by minimizing the distance between the observed atmospheric circulation and the regime states (see next subsection for a more detailed description of this procedure).Furthermore, we also allow external factors, like CO 2 , ozone, or ENSO, to influence the evolution of the regime states.This is again done by minimizing the distance.This also enables us to evaluate different forcing combinations.By using different forcing combinations and looking for the statistical model with the minimum distance, we can systematically find the forcings which are most likely responsible for the observed evolution of the regime states.

Overview of the FEM-BV methodology
The FEM-BV-VARX approach is a general variational framework that is reduced to the well-known methods of linear regression, autoregressive models, k means and hidden Markov approaches when more restrictive assumptions are made on the nature of the underlying data-generating process (Lean and Rind, 2008;Bromwich et al., 2013;Roscoe and Haigh, 2007;Metzner et al., 2012).Furthermore, as VARX is a tool for inferring the Granger causality (Granger, 1988) (causation between time series variables in terms of predictability and not correlation), FEM-BV-VARX is a more general approach which allows going beyond the standard stationarity assumption of the usual methods currently used for inferring cause-response relationships, e.g., in ecology (Sugihara et al., 2012), economics (Granger, 1988) and climate science (Mosedale et al., 2006;Wang et al., 2004).The non-stationary FEM-BV-VARX framework contains the standard stationary VARX and the concept of standard Granger causality as a particular special case and allows for a systematic comparison of causality relations inferred with and without the stationarity assumption for a given set of observational data.See Gerber and Horenko (2014) for more details explaining the relation of the FEM-BV framework to standard stationary approaches of data-driven causality inference.
The approach we are following here is that we perform FEM-BV-VARX fits to all possible combinations of the external forcings.Then we apply a standard information theoretic criterion: the AIC (Akaike, 1998;Horenko, 2010;O'Kane et al., 2013;Risbey et al., 2015), which is a measure of the relative goodness of fit of a model to data.That means that AIC can be used to assess the goodness of fit relative to the number of fitted parameters or external factors used, allowing us to find the model that is least overfitting and best fitting the analyzed data (Burnham and Anderson, 2002).This standard procedure can tell us which forcing combination is best able to explain the secular circulation trends given by the FEM-BV-VARX model.
The presence of unresolved external covariates (which are not statistically independent or identically distributed) may result in the non-stationarity and non-homogeneity of the resulting data-driven statistical models and may be manifested in the presence of secular trends and/or in regime-transition behavior.By covariate we not only mean external forcings but also unresolved physical processes and scales (e.g., due to EOF truncation).This may then introduce problems when applying the standard stationary approaches common to machine learning and statistics (de Wiljes et al., 2014).In the context of this paper, this issue plays a very important role when analyzing atmospheric data since many of the potentially relevant covariates might not be available explicitly in the set of covariates that we have chosen for testing.Therefore, when deploying statistical time series analysis methods, they should be capable of dealing with non-stationarity and non-homogeneity issues that emerge in the models as a result of these systematically missing (and potentially important) external influences.
Combining the concepts and ideas from pure and applied mathematics (such as the FEM from numeric of partial differential equations, regularization in infinitely dimensional spaces from the theory of ill-posed problems, stochastic calculus and theory of stochastic processes, information criteria from information theory, and embedding theorems from the theory of dynamical systems), Horenko and colleagues developed a family of non-stationary, non-homogeneous and non-parametric time series analysis methods.This family of time series analysis techniques, which is reviewed concisely by Metzner et al. (2012), allows for systematic timedependent model identification when assumptions of temporal stationarity or spatial homogeneity of some underlying statistics are not justifiable.The main idea is based on regularized variational minimization of a scalar-valued functional describing the error g (x(t), u(t), θ (t)) of some model for a given observation x(t) subject to available external impacts/covariates u(t) and characterized by the timedependent set of model parameters θ (t) subject to the constraints on (t) = (γ 1 (t), . .., γ K (t)): where θ (t) = K k=1 γ k (t) k , thereby allowing the problem of statistical data analysis beyond the usual stationarity assumption to be reformulated and algorithmically solved as a clustering problem with K clusters, distance kernel g (x(t), u(t), θ (t)) and a regularization parameter N C , controlling the number of transition between different clusters in time (special case N C ≡ 0 and K ≡ 1 corresponds to the homogeneous and stationary statistical problem).Imposed regularization confines the bounded variation (BV) of the regime-switching process γ k (t) in time, thereby making the temporal change of inferred model parameters θ (t) more or less persistent.Changing the constraining variable N C , one can test the whole range of possible statistical models, going from N C ≡ 0 (stationary/well-posed) to N C ≡ N T (completely non-stationary/very ill-posed problems).The above variational problem is non-convex since the parameters k are not known a priori and have to be inferred simultaneously with the γ k (t).
Many classical methods of data analysis and machine learning (e.g., multilinear statistical regression, k-means clustering, Gaussian mixture models (GMMs) and hidden Markov models (HMMs) Majda et al., 2006;Franzke et al., 2008) can be derived as special cases of this FEM-BV methodology.For example, HMMs are obtained additionally assuming that is an output of a homogeneous Markov chain and setting N C ≡ N T .For more detailed information please see Sect.2h "Relation to classical methods of unsupervised learning" in Metzner et al. (2012).It is an important feature of the deployed methodology, since it allows us to test different standards or more advanced methods in the context of the same theoretical and algorithmic FEM-BV framework.
In the FEM-BV methodology, finite element methods are employed in the numerical representation of indicator functions γ k (t) for the time domain of applicability of different models from a common model class.Model class is defined by the choice of the particular analytical form of the error function g (x(t), u(t), θ (t, j )); the explicit VARX-form deployed in this paper is given below.As shown above, these indicators are regularized using a bounded variation criterion, hence the acronym FEM-BV.The choice of the model class depends on the type of data and a specific form of error function g considered.Here we have implemented vector autoregressive models with external influences (VARX), which are defined as x t = µ(t)+ MEM i=1 A i (t)x t−iτ +B(t)φ (u(t))+ t , parameters being θ (t) = (µ(t), A 1 (t), . .., A MEM (t), B(t)), and the model error functional defined as g (x(t), u(t), k ) 2 (for more details please see Metzner et al., 2012;Horenko, 2010).
The algebraic structure of the above problem allows us to deploy efficient numerical algorithms based on the iterative application of linear and/or quadratic programming problems (optimizing for fixed ) followed by stationary/convex inference of for a fixed .This procedure is repeated iteratively, until the change in L( , ) is less than some small a prioridetermined threshold, resulting in monotonic minimization of L.Here we have used the adaptive finite element method from the numerics of partial differential equations (PDEs) (Horenko, 2010;de Wiljes et al., 2014) for numerical minimization of the FEM-BV problems (Eqs.1, 2).
The number of different spatiotemporal K regimes/clusters -the model parameters to be chosen within these regimes, such as memory depth and number of EOFs, and the indicator functions γ k (•) signaling activation of the respective models -are all determined simultaneously in a global optimization procedure.This yields a judicious compromise between low residuals, reproducing the data of a training set on the one hand, and the demand for the smallest-possible overall number of free parameters of the complete model on the other.The optimization is based on a new non-parametric modified AIC and may thus be interpreted as a constructive implementation of "Occam's razor" for data analysis problems.Thereby, the resulting FEM-BV-framework is essentially free of parameters that should be defined and tuned by the user.The only parameter that needs to be set externally is the overall number of optimization repetitions with different randomly chosen initial values of or for parameter optimization (in the following referred to as the number of annealing steps).Increasing this number reduces the probability of getting trapped in one of the local minima of L (for N C > 0), simultaneously linearly increasing the amount of computations.Therefore, the number of annealing steps should be chosen carefully, depending on the available computational resources and the size of the data to be analyzed.

Application of FEM-BV-VARX to atmospheric reanalysis and observed forcing data
We have chosen to examine a time series of 500 hPa geopotential height anomalies (seasonal cycle subtracted) projected on the 20 leading EOFs and with u(t) taken as 32 combinations of forcings (to be described in more detail below).
Deploying the FEM-BV-VARX method, comprehensive sensitivity tests were carried out for the cluster parameters involving memory depth MEM and the number of annealing steps.Although the resulting optimal choice of forcing combinations has converged after 16 annealing steps, the model affiliation sequence already remained unchanged after 4 annealing steps.The results presented in the manuscript are for 64 annealing steps.The optimal memory depth τ was 2 days; however, the degree of memory tested ranged from 0 to 5 days.Daily forcing agents was spline interpolated with no lag apart from OMD.Every possible combination of forcing agents was considered including the observed OMD lagged by 0, 30, 60 and 90 days as well as a variant with a lag average of 365 days span (here we artificially introduce persistency into the OMD time series).The optimal external radiative forcing agent was found ultimately to correspond to the Cape Grim CO 2 time series.
In order to choose the optimal model M we first used the AIC: (3) Application of the AIC is equivalent to assuming that the scalar-valued squared model errors are χ 2 -distributed and that the vector-valued FEM-BV-VARX model errors are Gaussian, i.e., dependent on the residuals having a lognormal distribution.We tested this assumption using a nonparametric information-theoretic algorithm from Metzner et al. (2012) and found that for all of the model errors t the most optimal parametric family was indeed the log-normal distribution.
Additionally, the log-normal distribution was fitted to the model errors, the respective log-likelihoods computed and used to calculate the AIC for the non-stationary models.
The most informative non-stationary model that emerged using the AIC criteria (with posterior probability almost equal to one in each case) was the model with Cape Grim CO 2 and memory of 3 days.Comparable results were found using Akaike information criteria corrected (AICc) and the Bayesian information criteria (BIC).The two latter criteria also take the size of the statistics into account and are derived under very different mathematical assumptions than the AIC.This further confirmed our results, demonstrating that they are not induced by the implicit assumptions necessary for the information criteria applicability.
In our study we have included only a selected number of several external forcings.One might argue now that we have neglected some additional forcing (e.g., sea ice extent or the strength of the Antarctic Circumpolar Current) which might be responsible for the observed secular trend.However, if we include something like sea ice extent into the set of forcings and get a result showing that sea ice extent is more statistically significant, then this would not contradict this study simply because the variable γ (describing the switching process and describing the bias coming from all the unresolved covariates in our study) would be different from the one that we obtained in this study.From our study we can guarantee that in a given set of explicit covariates we found the one covariate that is most important (in the Granger-causality sense) and -taking into account the presence of the unresolved covariates -this covariate is CO 2 (Metzner et al., 2012;de Wiljes et al., 2014).All of the other eventually important covariates are in the regime-switching process that we have also identified but were not explicitly represented in our results.

Attribution results
A previous study with the FEM-BV-VARX method (O'Kane et al., 2013) revealed the existence of statistically significant persistent circulation regimes corresponding to the positive phase of the SAM and a hemispheric wave-3 blocking pattern (Fig. 2).That study also found evidence for significant secular trends in all seasons and, when the full reanalysis period 1948-2010 was considered, that a distinct regime transition occurred around 1980 towards a preference for the positive SAM phase; here shown in Fig. 3 in terms of time of residence in either a wave-3 blocked or zonal (positive SAM) state.The increasingly frequent and persistent positive SAM pattern is accompanied by a corresponding decrease in the frequency and persistence of the wave-3 hemispheric pattern, which features blocking in the Australian sector in particular.The increasing frequency of occurrence of the zonal state and the corresponding decreasing frequency of occurrence of the wave-3 blocking state are consistent with the study of O' Kane et al. (2013) and are, thus, statistically significant.Furthermore, this behavior occurs throughout all four seasons and, thus, is very robust.However, the changes oc- cur slower in the summer season.While in all other seasons the wave-3 blocking state occurred more often than the zonal state up to about 1980; this transition has been much slower in the summer season and starting in 2000 both states occur about equally frequent in the summer season.There has been a much more distinct regime transition during the other three seasons (Fig. 3).
Our FEM-BV-VARX analysis finds strong evidence that anthropogenic greenhouse gas concentrations have caused the secular trends in the SAM and hemispheric wave-3 pattern.The clustering analysis with CO 2 forcing has the smallest AIC value, AIC min = −63 053 (which corresponds to a Akaike weight of about 1).This denotes the most parsimonious explanation of the observational data among all of the other fitted explanatory statistical models with all possible combinations of considered forcings.Hence, providing the best compromise between the quality of fit to the data and a low number of parameters.The absolute value of the AIC is less meaningful, only its relative size compared with the other tested models is useful.The next best forcing combinations relative to the optimal choice are the solar constant (Akaike weight of 2.3550 × 10 −5 ) and stratospheric aerosol (Akaike weight of 1.8217 × 10 −5 ).
The Akaike weight value w i = exp((AIC min − AIC(i))/2) provides a measure for how much better the best FEM-BV-VARX fit explains the data relative to the other FEM- BV-VARX models; this quantity can also be interpreted as the posterior model probability w i and, thus, as how less likely the FEM-BV-VARX model i is to explain the data relative to the best-fit model.The Akaike weight value (or posterior model probability) w i = exp((AIC min − AIC(OMD))/2) = 3.1083 × 10 −18 reveals statistically overwhelming support of CO 2 compared to stratospheric OMD in explaining the secular trend in the cluster frequency of occurrence.The corresponding Akaike weights reveal that the CO 2 forcing is significantly better than all other possible used combinations in explaining the observed trends.
We tested the sensitivity of our results using different information criteria like the BIC (Burnham and Anderson, 2002) and the AIC corrected for finite sample size (AICc) (Burnham and Anderson, 2002) and obtained very similar results in that CO 2 was always the forcing which best explained the secular trends.The next best fit varies according to the FEM-BV-VARX setup but it typically includes stratospheric ozone depletion and stratospheric sulfate aerosols.While our results are unclear as to the relative importance of stratospheric ozone depletion over the other leading forcings, we find strong evidence for the role of CO 2 in contributing to the secular trends.
Figure 2 also shows the geographical structure of the two cluster states in terms of 500 hPa geopotential height and surface air temperature (SAT) (note that the upper atmospheric regime states are very similar to the ones in O' Kane et al. (2013) which did not include external forcings in their analysis.).The wave-3 blocked state (state 1) is associated at the surface with a cold anomaly over the Antarctic Peninsula and a warm anomaly over the Ross Ice Shelf and the coast of Antarctica's Victoria land as well as with SAT over South America in SAT.The zonal positive SAM state (state 2) is associated at the surface with a warm anomaly over East Antarctica and Australia and a cold anomaly along the coast of Antarctica's Wilkes Land.The cluster states have very different surface temperature signatures.With the strong trend towards the SAM state in recent decades, it is interesting to ask whether the pattern of surface temperature trends over the Antarctic region reflect the SAM state surface signature.Figure 4 shows the trend in surface air temperature over the same period (1979-2010) calculated from yearly averaged Had4Krig version 2.0.0 (Cowtan and Way, 2014) data.Because large trends are evident in both blocking (wave 3) and SAM in all seasons (Fig. 3), we have used annual mean data to calculate the SAT trends.The remarkable agreement between the SAM state surface temperature anomaly pattern (Fig. 2d) and the Had4Krig SAT trend pattern over Antarctica is further evidence of the weakening of the wave-3 state and the shift towards the positive SAM state.
Our results are in contrast to earlier studies which found that ozone depletion is up to 9 times more important than anthropogenic CO 2 concentrations (Roscoe and Haigh, 2007;Son et al., 2010;Lee and Feldstein, 2013).One possible explanation for this difference is that the earlier studies mainly examined the austral summer season focusing on the response of the linear trend in the zonal mean circulation towards the positive SAM phase and the poleward shift of the Hadley cell (Son et al., 2010;Polvani et al., 2011;Previdi and Polvani, 2014).Our study focuses on attributing systematic changes in the circulation over the latter half of the NCEP reanalysis period employing a data-driven methodology that can infer causation (as explained above, FEM-BV-VARX is a non-stationary extension of the Granger causality inference and can describe the standard Granger causality as a particular stationary case).Moreover, we are not simply considering changes to the zonal SAM index in the austral summer in isolation but are explicitly attributing changes to the entire SH circulation including coherent features to all possible combinations of the relevant radiative forcings.Previous studies mainly analyzed changes in the mean state and not in the frequency of occurrence or changes in structure.This might also partly explain why our findings differ from previous studies.
To increase the confidence in our results, we systematically examined the sensitivity of our results to the treatment of the ozone data.Considering a 365-day-averaged and timelagged seasonally varying OMD leads to more robust results because we account for the strong annual cycle of stratospheric ozone and its delayed impact on the tropospheric circulation.Ozone has a strong seasonal component with OMD known to impact the tropospheric circulation (from the observational record) in December-January.Thus, we repeated our analysis using lagged (by 0, 1, 2 and 3 months), seasonally varying and 365-running-mean OMD data.While we did find some sensitivity to lag interval, our results were qualitatively unchanged.

Sensitivity to the number of model parameters
O' Kane et al. (2013) showed that, in order to accurately capture the dynamics and amplitude of Southern Hemisphere midlatitude atmospheric blocking regime transitions, any dimension reduction of the 500 hPa atmospheric reanalysis data required including a minimum of 9 but more generally the leading 20 PCs (principle components).Retaining these dominant modes makes the VARX model 20dimensional with the resulting FEM-BV-VARX, where each dimension communicates with every other dimension, leading to quadratic growth in the number of VARX parameters in the matrix A and of dimension n.This also applies to the stationary VARX model of the form where x t has the dimension n, u t has the dimension m and t is independent and identically distributed.Here the total number of model parameters will be N param = n+MEM•n 2 + n•m.With n = 20 and MEM ≥ 2 resulting in many thousands of free model parameters.Since the length of the available observational data is limited, a classical ill-posed problem is manifested in the overfitting of the data.This is the case even for the standard stationary data sets.BIC and AIC implicitly see this problem, attributing these models to much higher values of the information criteria than they would for more informative and well-posed models.A more complete discussion can be found in Metzner et al. (2012).
In order to try and address this problem, and given that we have ascertained that a minimum of 20 PCs should be retained, we have considered additional sensitivity experiments in which we diagonalize the matrix A i .Since the only coupling between the different dimensions of x t is induced by the off-diagonal elements of A i , diagonalizing is equivalent to a separate identification of the n following problems with respect to µ(j ), A i (j, j ), B(j, :); i.e., ∀j = 1, n : x t (j ) = µ(j ) (5) It is straightforward to see that the total number of parameters in such a case will be growing linearly with n as N param = n • (1 + MEM + m) (multiplied by K if considering non-stationary models where K > 1).In sensitivity experiments with A i (j, j ) and for memory depths ≥ 2, Cape Grim CO 2 was again found to be optimal and the results were insensitive to annealing steps ≥ 4. In descending order, and for the diagonalized experiments, the leading five combinations were found to be (1) Cape Grim CO 2 , (2) Cape Grim CO 2 and ozone, (3) optical thickness, (4) sulfate aerosols, (5) stratospheric ozone, and where a memory depth of 2 days was found to be optimal.We further note that even such a diagonally restricted VARX model is much more general than when applying multilinear regression: which is a special case of VARX for K = 1 and A j = 0∀j .Thus, even with the diagonality restriction, this analysis is more general than the stationary multilinear regression approaches, e.g., employed by Roscoe and Haigh (2007).One might further seek to reduce the number of parameters such that the problem is not ill posed.While reducing the number of EOFs is a possible approach even a reduction from 20 to 9 EOFs (the absolute minimum number of modes required to capture the Southern Hemisphere wave-3 blocking state) will not sufficiently reduce the number of parameters such that the FEM-BV-VARX with full A i matrix is well posed.Also in the Southern Hemisphere blocking is quite transient and, as shown in Fig. 11 of Zidikheri et al. (2007), the difference between equilibrium zonal and blocked states in terms of hemispheric zonal wind speeds averaged over a midlatitude zonal band is only 3 ms −1 at 500 mb as compared to the Northern Hemisphere which is 30 ms −1 an order of magnitude larger.
One strategy that was comprehensively tested was the sensitivity to persistency over a large range of values as N C from Eqs. (1) and ( 2) is changed from 0 to N T .However, due to the transient nature of the Southern Hemisphere atmospheric circulation, the residuals (model errors) for the multiple-state model are not significantly smaller than for the one-state model.By including the number of transitions in the definition of N param it is immediately obvious that the one-state model is always preferred in such cases and taking N C → N T we simply converge on the one-state solution.
Clearly time-series analysis where persistency of the respective metastable states is weak represents a serious challenge.One approach we explored assumes that -for fixed time-series x t , number of metastable states K, and persistency N C , as well as for fixed local VARX parameters θ i = (µ i , A i 1 , . .., A i MEM , B i ) -we can straightforwardly compute the respective Viterbi path solving the linear programming step of the FEM-BV procedure (see Step 2 of the FEM-BV algorithm description on p. 23 of Horenko, 2010).Since it is a linear minimization problem with convex constraints, it has a unique solution.That is, one can uniquely recover the distinct Viterbi path, (t) = (γ 1 (t), . .., γ k (t)), knowing only the full data series x t and preserving only + n • m) + 1) parameters while preserving the value of persistency C. Having computed the Viterbi path, one can also compute the distinct values of the model errors The sensitivity experiments we describe effectively bound the problem of overfitting inherent in analyzing atmospheric observational data.Importantly, we achieve the same results for diagonalization (well-posed, with no cross terms) and for the full FEM-BV-VARX (ill conditioned, with all cross terms).

Internal climate variability
Next we examined whether intrinsic climate variability statistically significantly affected the regime behavior.We find that the combination of ENSO and the first component of the multivariate MJO index (MJO1) provide the major intrinsic driver of the observed atmospheric regime behavior.The next best combinations are MJO1, MJO1 together with the eastern Indian Ocean Dipole index and the annual cycle together with MJO1.
The low frequency variability of ENSO is highly correlated to the IPO and, as pointed out earlier, the IPO shifted phase in the late 1970s coinciding with the transition to reduced blocking.Thus, it is natural to expect ENSO to be a major component of internal variability driving changes in the wave 3. The first component of the MJO index corresponds to enhanced convection over the maritime continent (Indonesia, Philippines and Papua New Guinea) close to the tropical warm pool (Wheeler and Hendon, 2004).This provides evidence for a tropical origin of SH mid-and high-latitude climate variability on interannual to decadal timescales.Since SST changes occur on longer timescales, this might open the opportunity for making skillful longrange predictions on seasonal to decadal timescales.However, the external CO 2 forcing still explains the observed secular trends best; thus, the intrinsic climate indices taken alone are not able to statistically explain the secular trends.We also find evidence that the Mt.Pinatubo volcanic eruption in 1991, as measured by stratospheric aerosol optical thickness, could have triggered a dramatic sudden increase in the regime frequency of occurrence in its immediate aftermath.Figure 5 shows that the eruption of Mt.Pinatubo is followed by a sudden drop in the frequency of occurrence (1-year running mean) of the positive phase of the SAM and a corresponding increase in blocking.In the long term this only delayed the secular increase in the SAM.From Fig. 5 we infer also that the response timescale to the eruption is about 3-4 years.In contrast, the 1982 eruption of El Chichón did not cause a drop of the frequency of occurrence of the positive SAM phase.This result is consistent with the EOF analysis of ERA-40 reanalysis data where a significant shift to negative stratospheric and surface SAM was observed only after the Fuego and Mt.Pinatubo eruptions (Karpechko, 2010;Gleixner, 2012).Importantly, CMIP5 models have been shown to be unable to reproduce a realistic dynamical response by the annular mode to even large intermittent volcanic signals like Mt. Pinatubo, suggesting that the extratropical circulations of current CMIP5 models are not able to simulate the response to short-lived abrupt perturbations in stratospheric forcing (Karpechko, 2010;Gleixner, 2012;Charlton-Perez et al., 2013).

Conclusions
Our examination of reanalysis data together with observed forcing data reveals that greenhouse gas emissions are an important driver of SH circulation changes over the last few decades.Recent studies have suggested that stratospheric ozone depletion is many times more dominant than CO 2 in driving systematic changes in the SH midlatitude circulation.However, our results highlight that for understanding the anthropogenic impact on SH circulation changes the delayed influence of stratospheric ozone depletion is relevant but not the dominant mechanism.Previous studies mainly focused on the austral summer and the zonal response (SAM).Our results from analyzing data from the whole year suggest that other seasons and the changes in coherent features (wave-3 blocking), including implications for anomalous surface warming of the high latitudes, need to be incorporated in order to more completely understand and attribute SH circulation changes.Such studies are particularly needed for evaluating and improving climate models.Temperature trends in the Antarctic region are spatially heterogeneous.The pattern of these trends reflects the surface signature of the shift toward a more zonal SAM circulation regime at the expense of the blocking regime.
Our observationally based results are also confirmed by numerical modeling studies (Miller et al., 2006;Freitas et al., 2013Freitas et al., , 2015)).Atmospheric general circulation model simulations forced by observed SST fields (Rayner et al., 2003) and CO 2 concentrations but with climatological O 3 have been able to reproduce the observed SH circulation changes although the magnitude is underestimated by about 40-50 % (Freitas et al., 2013(Freitas et al., , 2015)).This provides support from a numerical climate model for the observational data analysis presented in this manuscript.
Our finding that anthropogenic CO 2 is the dominant driver whereas stratospheric ozone depletion makes a somewhat lesser contribution to the SH circulation changes over recent decades has important implications for future SH climate change.In particular, it may be that the recovery of ozone has less relevance to changes in Southern Hemisphere extratropical circulation than projected in many modeling studies (Barnes et al., 2013;Shindell et al., 2004).Our findings regarding the sudden but short-lived increase in blocking and negative SAM after the Mt.Pinatubo eruption highlights the potential of the climate system to abruptly change in response to large transient perturbations in stratospheric forcing while emphasizing the dominant role of systematic changes in anthropogenic CO 2 on the climate.

Figure 1 .
Figure1.Forcing time series: Cape Grim CO 2 (dark blue), sulfate aerosols (green), stratospheric aerosol optical thickness (red), lagged OMD (blue), OMD (magenta) and solar constant (khaki).Time series are normalized by subtracting the respective mean and dividing by the respective standard deviation.

Figure 3 .
Figure 3. Percentage of time in either the hemispheric wave-3 state (black dashed) or the zonal state (blue dashed) for the NCEP reanalysis 500 hPa geopotential height field for all seasons and annually.The dashed lines are a LOESS fit to the time-averaged data where the solid lines indicate the values and averaging periods of the data.

Figure 4 .
Figure 4. Surface air temperature trend (K decade −1 ) over the period 1979-2010 calculated from the yearly averaged Had4Krig version 2.0.0 data set.