Interpreting the nature of Northern and Southern Annular Mode variability in CMIP5 Models

Characteristic timescales for the Northern Annular Mode (NAM) and Southern Annular Mode (SAM) variability are diagnosed in historical simulations submitted to the Coupled Model Intercomparison Project Phase 5 (CMIP5) and are compared to the European Centre for Medium‐Range Weather Forecasts ERA‐Interim data. These timescales are calculated from geopotential height anomaly spectra using a recently developed method, where spectra are divided into low‐frequency (Lorentzian) and high‐frequency (exponential) parts to account for stochastic and chaotic behaviors, respectively. As found for reanalysis data, model spectra at high frequencies are consistent with low‐order chaotic behavior, in contrast to an AR1 process at low frequencies. This places the characterization of the annular mode timescales in a more dynamical rather than purely stochastic context. The characteristic high‐frequency timescales for the NAM and SAM derived from the model spectra at high frequencies are ∼5 days, independent of season, which is consistent with the timescales of ERA‐Interim. In the low‐frequency domain, however, models are slightly biased toward too long timescales, but within the error bars, a finding which is consistent with previous studies of CMIP3 models. For the SAM, low‐frequency timescales in November, December, January, and February are overestimated in the models compared to ERA‐Interim. In some models, the overestimation in the SAM austral summer timescale is partly due to interannual variability, which can inflate these timescales by up to ∼40% in the models but only accounts for about 5% in the ERA‐Interim reanalysis.


Introduction
The ability of global climate models to accurately represent the key aspects of climate variability enhances confidence in projections of future climate change. A particular focus of climate research has been the Northern and Southern Annular Modes (NAM and SAM, respectively), which are understood to represent dominant patterns of extratropical atmospheric variability on short and long timescales [Thompson and Wallace, 2000;Baldwin, 2001]. The annular modes are defined as the first empirical orthogonal function (EOF) of geopotential height or sea level pressure in the Northern and Southern Hemisphere. Annular mode variability has been linked with different weather and climate-related phenomena, such as precipitation [Field, 2010], sea ice distribution [Yuan and Li, 2008], and surface air temperature [Karpechko et al., 2009]. The annular modes are thought to be influenced by ocean sea surface temperature [Sen Gupta and England, 2007], stratospheric variability [Baldwin and Dunkerton, 2001;Simpson et al., 2011], changing greenhouse gas concentrations [Arblaster and Meehl, 2006], and ozone concentrations [Sexton, 2001], making them a useful identifier of future climate change.
In this study, Coupled Model Intercomparison Project Phase 5 (CMIP5) model simulations are tested for their ability to reproduce the timescales of the observed NAM and SAM following methods outlined in Osprey and Ambaum [2011] (henceforth OA11). OA11 suggest dividing the near-surface annular mode power spectrum into two parts, from which characteristic timescales, identifiable with high and low-frequency variability, ( low and high , respectively) can be calculated. For the high-frequency spectra, an exponential fit is made which is associated with chaotic processes [Sigeti, 1995;Ohtomo et al., 1995;Frisch and Morf , 1981]. For low frequencies, a first-order autoregressive process (AR1) is fit, which parameterizes possible stochastic variability [Keeley et al., 2009;Feldstein, 2000;Hasselmann, 1976]. These methods differ from previous studies, such as Baldwin et al. [2003] and Gerber et al. [2008b], who fit an exponential to the autocorrelation function, and thereby assume a purely stochastic description of annular mode variability across all timescales. 4.2 ± 0.8 5.7 ± 0.3 34 4.5 ± 0.6 7.6 ± 0.6 20 a The spectral cut frequency c is given in cycles per year. The upper values for each model are NDJF timescales calculated from the spectrum centered at 15 January/1 December for NH/SH, and the lower MJJA timescales are centered at 15 July/1 June for NH/SH. The errors are the ones from the fitting procedure Δ f . All timescales were calculated from spectra omitting power in frequencies below 1 yr −1 . A dynamical framing of annular mode timescales opens up the possibility of using dynamical systems approaches for better understanding and forecasting the future annular mode response to external forcing (e.g., increases in greenhouse gas concentrations).
Section 2 outlines those models and procedures used to construct the annular modes and characteristic timescales from the CMIP5 model ensemble and reanalyses. Section 3 identifies and describes summertime (May, June, July, and August (MJJA)) and wintertime (November, December, January, and February (NDJF)) timescales in the models and also quantifies the contribution of interannual variability to the low-frequency timescale. We conclude with a summary and discussion of our results in section 4.

Data and Methods
For this study, daily geopotential height data are taken from the historical r1i1p1 run of 15 different CMIP5 models, namely, BCC-CSM1.  Table 1 and institution information in Table 2). Additionally, the analysis was performed on three different reanalysis data sets, namely, ERA-40, ERA-Interim, and JRA55 (institution information in Table 3).
We define the NAM/SAM as leading order EOFs of the NH/SH extratropics. For reconstructing the time evolution of the annular modes, principal component time series are calculated of the leading empirical orthogonal function (EOF) of 850 hPa daily mean zonally averaged geopotential height. For the Northern Hemisphere (NH), the leading EOF is calculated between 20 ∘ N and 90 ∘ N, while in the Southern Hemisphere (SH) geopotential height data between 20 ∘ S and 85 ∘ S are used. The latter is due to lack of data close to the South Pole. Before performing the EOF analysis, data were deseasonalized and the linear trend removed. Following Baldwin and Dunkerton [2001], EOFs are constructed on the 90 day smoothed full time series, and annular mode time series are formed by the projection of the unsmoothed height anomalies onto the leading EOF.
Characteristic NAM/SAM timescales are derived from the respective power spectra using the methods outlined in OA11. Power spectra P of the different seasons are calculated by weighting the annular mode time series with a Gaussian (standard deviation = 45 days) centered at 15 January/July for the NAM and 1 December/June for the SAM and then applying a Fourier transform. We will refer to the spectra calculated from the weighted time series as NDJF and MJJA. A least squares fit P fit ( ) to the power spectrum P is performed, before using an inverse Fourier transform of P fit to obtain a fitted autocovariance function R fit ( ) =  −1 (P fit ( )), which is then normalized by its value at lag-zero, recovering the autocorrelation function (ACF).
For calculating P fit ( ) the power spectra are divided into two parts; the low-frequency part ( < c ) is fitted by a Lorentzian, the high-frequency part ( ≥ c ) by an exponential function, The parameters A and B are identified as relating to the integrated power across the piecewise spectra, and and are linked to the characteristic low and high-frequency timescales, as described below. The cutoff frequency c , which loosely represents the frequency separating weather variability and its integrated low-frequency effects, is adapted to optimize the fit R fit ( )∕R fit (0) to the actual autocorrelation function of P.
To obtain that, c was varied between 5 and 50 cycles per year and the value chosen for which the difference at 1∕e between fit and actual autocorrelation function was minimized.
Characteristic timescales, relating to the low and high-frequency parts, are then retrieved as low = 1 for < c and high = 4 for ≥ c as in OA11. Chi-square related uncertainty estimates on these timescales are calculated by propagating the errors in the fitted parameters, Δ and Δ (1 error, calculated from the covariance matrix) to the timescales as Δ low Additional uncertainty estimates of the method were obtained by varying the cut frequency c by ±20% as in OA11. Further, a subsampling test of the time series was carried out, using 100 samples of 30 randomly chosen years (bootstrapping with replacement), and the standard deviation in the random sample of the retrieved timescales, Δ s , was calculated. The CMIP5 annular mode timescales are compared with those derived for ERA-Interim.
Finally, the possible influence of external interannual variability (eIAV) on the characteristic timescales of the low-frequency power spectra is explored. Here "external" means that the power is not accounted for by an AR1 process calculated using frequencies between 1 yr −1 and c . This is achieved by omitting those frequencies lower than 1 cycle per year in the fitting procedure ( Figure 1). Therefore, the Lorentzian fit only retains interannual variability consistent with an AR1 process defined on intraannual timescales. Hemisphere. For comparison, spectra were scaled to 1 at f = 0. Single-model spectra are depicted in grey, a mean spectrum in black. Also shown are Lorentzian (blue) and exponential (red) fits to the mean spectrum that characterize AR1 and chaotic processes (solid for the according part of the spectrum, dashed otherwise). It can be clearly seen that the spectra cannot be characterized by either of these functions alone.    The reported CMIP5 models show similar behavior at high frequencies, leading to characteristic timescales between 4 and 6 days for both summer and winter in the NH and SH, a result consistent with similar features found in reanalysis data, discussed in OA11. As seen in OA11 for reanalysis data, for the CMIP5 models studied, the shape of the spectra over high frequencies is consistent with an exponential function ( Figure 2). However, the behavior at low frequencies differs a lot across the models, giving NAM timescales between 6.3 to 11.6 days in NDJF and 6.2 to 10.5 days in MJJA. These NAM timescales compare with 8.7 days in NDJF and 5.9 days in MJJA, respectively, in ERA-Interim. The spread in the timescales for the SAM are even higher: 7.5 to 18.0 days in MJJA, and 14.4 to 30.3 days in NDJF, with the multimodel mean being 10.3 and 19.7 days (ERA-Interim: 7.4 and 11.2 days, respectively). In Figure 3 the ensemble mean retrieved fit of the autocorrelation functions is shown as well as the fit to ERA-Interim data for NH and SH. The mismatch between the autocorrelation of the fits and the actual autocorrelation shows the overestimation of the timescales that are particularly large in SH NDJF. These model biases are consistent with previous studies using a different method for calculating the timescale, such as Gerber et al. [2008aGerber et al. [ , 2010. Taking various error estimates into account; however, the biases are unlikely to be outside the error bars for the NAM timescales.

Results
Tables 4 and 5 show the different errors for all the timescales. One possible uncertainty in the timescales arises from the choice of c . For high , changing c has only minor influence, whereas for low the error can be of the order of days. The uncertainties are greater for higher timescales, especially in SH NDJF, and for very low cutoffs (<10 cycles per year). The latter is likely to be due to fitting the Lorentzian to a reduced number of data points when c is small. Compared to the other errors, Δ is usually the biggest and the only one that exceeds 10 days for some of the models. However, this error is only of this order for SAM NDJF timescales; in this case, the spectra at low frequencies exhibit a very steep rise (Figure 2). The higher c , the less will this power influence the fit, which gives a rather steep decline of high with c and thus large errors. Errors from the fitting procedure Δ f are relatively small, usually below 1 day-the exception again being SH NDJF. The same is true for errors estimated by using a subsampling test Δ s .
Different behavior of the models can be assessed by looking at the power spectra ( Figure 2). The spectra were scaled to match 1 at = 0 for comparison. They show the two regime behaviors that were also noted by OA11 for reanalysis data, exhibiting a steep rise at low frequencies and dropping exponentially at high frequencies.
Although the slope at low frequencies can be different, the behavior at high frequencies is very similar across all models, which can also be seen by the agreement of the high-frequency timescale high (Figure 4). Also shown in Figure 2 is a mean of the model spectra and the Lorentzian and exponential fit to it. These fits apply well in their respective frequency domain but are both not well representing the other part. The Lorentzian fit does not represent the spectra at high frequencies, being too shallow and the exponential fit does not adequately represent the steep rise at low frequencies.
The power in the very low (≲1 cycle per year) frequencies defines interannual variability. The effect of eIAV, which is not characterized by an AR1 process, on the fit can be seen in Figure 1. Taking all frequencies into account, the fit shows a similar steep rise as the data but a poor agreement at higher frequencies and leads to a higher characteristic timescale. Neglecting frequencies below 1 cycle per year in the fitting procedure, in contrast, gives a better estimate of the data and less discontinuity to the exponential fit. The slower rise then accounts for a lower characteristic timescale low .
An overview of the influence of the eIAV on the low-frequency timescale low is given in Figure 5. The blue bars show the timescales when frequencies less than 1 cycle per year are omitted; the underlying red bars show the timescales obtained by fitting to the whole spectrum. In the Northern Hemisphere, eIAV only accounts for up to 10% in MJJA and 13% in NDJF, whereas in the Southern Hemisphere, the contribution can be as high as 19% 10.1002/2014JD022989 Figure 5. Characteristic low-frequency timescales low in days for the Northern and Southern Annular Modes in summer and winter. Blue bars show the timescale when interannual variability is neglected, underlying red the timescale when taking the spectrum at frequencies of less than 1 cycle per year into account as well. For NAM, the weighting Gaussian was centered at 15 January (NDJF) and 15 July (MJJA); for SAM the centers are at 1 December and 1 June, respectively.

Discussion and Conclusion
Northern and Southern Annular Mode timescales for 15 different CMIP5 models were investigated by looking at the power spectra and autocorrelation function at 850 hPa. In contrast to other studies, the power spectra are divided into two parts, fitted by an exponential function and a Lorentzian, which leads to two different timescales, high and low . As for reanalysis data in Osprey and Ambaum [2011], model spectra at high frequencies exhibit low-order chaotic rather than stochastic behavior, with timescales of 4-6 days in both hemispheres and seasons. As this coincides with the timescale for the growth and decay of baroclinic eddies, it might well be that the high-frequency part of the spectrum is dominated by the effect of individual eddy lifecycles on the mean flow. This contribution is not taken into account in simpler stochastic models that assume an AR1 process behind the annular mode dynamics. A similar behavior of the spectra was found in a study of a simplified atmosphere [Kunz and Greatbatch, 2014], who link the decaying spectra at high frequencies to interaction of a frictional boundary with mechanical forcing and thermal damping. As their model does not include any nonlinearity, spectra at even higher frequencies would-for algebraic reasons-not show exponential decay [Sigeti and Horsthemke, 1987]. So even if the spectra they show resemble those presented in this study, theirs cannot exhibit an exponential form over a larger frequency range; thus, Kunz and Greatbatch [2014] are not providing a counterexample for the chaotic interpretation of the spectra shown here, which are consistent with an exponential decay. Kunz and Greatbatch [2014] themselves do not rule out nonlinear contributions to the spectrum, specifically from baroclinic eddies.
Model output was compared to values retrieved from ERA-Interim data. At high frequencies, spectra show similar behavior for models and reanalysis data, for the Northern and Southern Hemisphere which is independent of time of year, with characteristic timescale high of about 4.5 days. At low frequencies, however, behavior is different, with models being biased too long in the SH NDJF timescales.
For the Northern Annular Mode, ERA-Interim data suggest characteristic timescales of about 8.7 days in winter and 5.9 days in summer, whereas models on average give 8.4/8.2 days, respectively. However, taking various uncertainty sources into account (Table 4), the disagreement in the Northern Hemisphere is unlikely to be significant, which has also been pointed out by Simpson et al. [2011].
For the Southern Annular Mode, low-frequency characteristic timescales are overestimated for both NDJF and MJJA, model average being 19.7/10.3 days, whereas ERA-Interim data give 11.2/7.4 days, respectively. The values retrieved from the three different reanalysis data sets are consistent but are not identical to the ones given in Osprey and Ambaum [2011], especially in the Northern Hemisphere winter. This can be attributed to a couple of differences: The cutoff frequency c used here is higher, which results in a smaller influence of power at very low frequencies, as there are more data points to fit to. Further, OA11 do not remove eIAV. They also use a time series starting in 1958 for calculation of their spectra and look at the 1000 hPa level. Calculation of the NAM NDJF timescale for ERA-40 at 1000 hPa following the procedure presented here gives low = 10.1 days, which is comparable to the other reanalyses data sets.
Another main difference between NAM and SAM is the influence of "external" interannual variability. Neglecting the spectrum below 1 cycle per year for the Lorentzian fit leads to a decrease in low-frequency timescale for 13 out of 15 models in the Southern Hemisphere NDJF, which can be as high as 39%, whereas the changes in the Northern Hemisphere are comparably small.
As the SAM is strongly influenced by ozone concentration [Sexton, 2001;Roscoe and Haigh, 2007;McLandress et al., 2011], this difference could be due to circulation effects following different ozone levels in both hemispheres. Another major influence on the SAM is El Niño-Southern Oscillation [e.g., Gong et al., 2010], which is stronger in December, January, and February, and may account for larger timescales as well as larger uncertainties. Further, a study looking at eddy feedbacks [Simpson et al., 2013] links the overestimation of the SAM NDJF timescale to a deficiency in planetary wave feedback in CMIP5 models. Eddies also seem important in establishing and maintaining a certain jet position [Kidston et al., 2010], which in turn influences the timescale of the annular modes. The connection between jet stream position and the annular modes is a robust feature in general circulation models (GCMs) [e.g., Barnes and Polvani, 2013]. This also means that dynamic behavior of the jet could be reflected in the annular modes; power in the low frequencies might be explained by the Hurst effect from regime behavior of the jet latitude, as described by Franzke et al. [2015].
The small changes in the NH seem to contradict Keeley et al. [2009] who claim that removing IAV reduces the NAM timescale substantially. However, their method differs regarding the treatment of interannual variability. Keeley et al. [2009] aim to remove all externally forced IAV, whereas this study only reduces its influence on the very low-frequency part of the power spectrum. Further, their timescale is derived differently by calculating it from the autocorrelation function, which consequently leads to different values. An advantage of our current method is that it will minimize the influence of numerical artifacts in retrieved ACFs from any uncritical removal of low-frequency signals in power spectra. Despite those differences, the two studies agree on the claim that NAM timescales are similar in both seasons after reduction of IAV.
As the annular modes have a major impact on surface climate, e.g., precipitation in European winter as indicated by Li and Wang [2013] (who also point out the correlation between sea ice cover and the NAM) it is crucial to represent them adequately in models in order to get a reliable prediction of more local phenomena. So far especially the SAM seems not to be well captured regarding the behavior at low frequencies (≲20 cycles per year), however, they are well represented at higher frequencies. This different behavior is consistent in all the CMIP5 models which justifies the use of two different functions to fit the annular mode power spectra and indicates different underlying processes. This distinction is not captured by studies that only define one timescale based on the autocorrelation function as it only shows in the power spectra. Therefore, the method presented here better accounts for the possible underlying dynamics of the annular modes. However, for simply identifying biases in the timescales, the findings for the low-frequency timescale in this study agree with previous simpler studies regarding seasonal and hemispheric differences as well as the overestimation of the