Nonlinear Processes in Geophysics ENSO ’ s non-stationary and non-Gaussian character : the role of climate shifts

El Niño Southern Oscillation (ENSO) is the dominant mode of climate variability in the Pacific, having socio-economic impacts on surrounding regions. ENSO exhibits significant modulation on decadal to inter-decadal time scales which is related to changes in its characteristics (onset, amplitude, frequency, propagation, and predictability). Some of these characteristics tend to be overlooked in ENSO studies, such as its asymmetry (the number and amplitude of warm and cold events are not equal) and the deviation of its statistics from those of the Gaussian distribution. These properties could be related to the ability of the current generation of coupled models to predict ENSO and its modulation. Here, ENSO’s non-Gaussian nature and asymmetry are diagnosed from in situ data and a variety of models (from intermediate complexity models to full-physics coupled general circulation models (CGCMs)) using robust statistical tools initially designed for financial mathematics studies. In particular α-stable laws are used as theoretical background material to measure (and quantify) the non-Gaussian character of ENSO time series and to estimate the skill of “na ı̈ve” statistical models in producing deviation from Gaussian laws and asymmetry. The former are based on non-stationary processes dominated by abrupt changes in mean state and empirical variance. It is shown that the α-stable character of ENSO may result from the presence of climate shifts in the time series. Also, cool (warm) periods are associated with ENSO statistics having a stronger (weaker) tendency towards Gaussianity and lower (greater) asymmetry. This supports the hypothesis of ENSO being rectified by changes in mean state through nonlinear processes. The relationship between Correspondence to: J. Boucharel (julien.boucharel@legos.obs-mip.fr) changes in mean state and nonlinearity (skewness) is further investigated both in the Zebiak and Cane (1987)’s model and the models of the Intergovernmental Panel for Climate Change (IPCC). Whereas there is a clear relationship in all models between ENSO asymmetry (as measured by skewness or nonlinear advection) and changes in mean state, they exhibit a variety of behaviour with regard to α-stability. This suggests that the dynamics associated with climate shifts and the occurrence of extreme events involve higher-order statistical moments that cannot be accounted for solely by nonlinear advection.


Introduction
El Niño Southern Oscillation (ENSO, see the glossary for the acronyms list) is the dominant mode of climate variability in the Pacific (MacPhaden et al., 1998).It impacts many surrounding regions and has major socio-economic consequences.Although our knowledge of the phenomenon has increased considerably in the last two decades, ENSO remains difficult to predict and its characteristics change in ways that are not yet understood by the scientific community.In particular, ENSO's characteristics (frequency, amplitude, propagating features and predictability) vary with changes in the mean state of the tropical Pacific (Fedorov and Philander 2000;Moon et al., 2004;An, 2004).The difficulty in predicting ENSO and its evolution lies partly in the limited ability of Gaussian statistics to account for Extreme Events (EEs).In fact most studies of ENSO implic- .A Gaussian curve that corresponds to the best fit to the PDF is overlapped in dashed thick line.The bins for SST anomalies larger (lower) than 2°C (-2°C) are duplicated (red shading) and represented on a scale with a 1/10 ratio (right scale) to highlight the asymmetry between the 'negative' and 'positive' tails of the PDF.
Fig. 1.K98 Dataset histogram constructed with a bin equal to max(X N ino3 )−min(X N ino3 )

19
= 0.3063 • C. A Gaussian curve that corresponds to the best fit to the PDF is overlapped in dashed thick line.The bins for SST anomalies larger (lower) than 2 • C (−2 • C) are represented on a scale with a 1/10 ratio (right scale) to highlight the asymmetry between the "negative" and "positive" tails of the PDF (red shading).
underestimated.In addition, recent studies have pointed out that Sea Surface Temperature Anomalies (SSTAs) over the eastern Pacific are positively skewed due to the nonlinearities of the tropical Pacific ocean-atmosphere system (Burgers and Stephenson, 1999;Hannachi et al., 2004;An and Jin, 2004).Thus ENSO has been depicted as a non-stationary and asymmetrical phenomenon (An and Jin, 2004;An et al., 2005) that can be rectified by changes in mean state (Rodgers et al., 2004;Dewitte et al., 2007a).The latter vary within decadal to inter-decadal time scales, partly reflecting the occurrence of abrupt transitions, named "climate shifts" (Trenberth and Hurrel, 1994;Zhang et al., 1997;Guilderson and Schrag, 1998;Urban et al., 2000).The source of these climate shifts remains unclear.Whereas some authors argue that extra-tropical variability can produce changes in tropical mean state through atmospheric teleconnections (Pierce et al., 2000) or oceanic "tunnels" (Gu and Philander, 1997), others suggest the importance of nonlinear processes within the tropics, in producing decadal variability and ENSO modulation (Timmermann and Jin, 2002;Timmermann, 2003;Timmerman et al., 2003;An and Jin, 2004;Dewitte et al., 2007a).
In this study, the focus was on ENSO statistics and their relationship with changes in mean state.However, unlike the aforementioned studies, the non-Gaussian nature of ENSO was explicitly taken into account.This property was diagnosed using relevant mathematical tools.Figure 1  Gaussian PDF overlapped (K98 data fitted).c) Empirical variance of K98 (solid line) and empirical variance of the Gaussian fitted distribution (dash line).d) Asymptotical test (g(T) as in section 2.2.1) for K98 data (solid line) and for the Gaussian fitted distribution (dash line) (see their description in Section 2.2.1.).  ) as in Sect.2.2.1) for K98 data (solid line) and for the Gaussian fitted distribution (dash line) (see their description in Sect.2.2.1.).plan et al., 1998).Figure 1 shows that the PDF had a small but significant deviation from Gaussianity.As an indication, the PDF of the Gaussian law fitted to the data is also plotted in Fig. 1.Note that in terms of information theory, the less Gaussian the PDF the more information it contains.The underlying question is then: within a simple theoretical framework, what causes ENSO's non-Gaussian and nonstationary character?More specifically, to what extent can climate shifts account for this particular ENSO property and are they part of the process of rectification of ENSO variability through the slowly varying mean state identified in earlier studies?
In a recent study, Hannachi et al. (2004) addressed a similar issue.Their approach was based on a nonlinear stochastic model to derive the nonlinearity associated with the NINO3 index.In their Fig. 15, the authors compared different "Lmoments" (equivalent to normalized statistical moments) of the NINO3 index.The different scatter plots displayed in the figure demonstrate that they found no significant relationship between the mean state and the interannual spread of the NINO3 SSTAs in the 24 models of ENSIP (the El Niño Simulation Intercomparison Project).Furthermore, it turned out that the majority of the models tended to concentrate in a cluster around the Normal distribution.Unlike recent physical studies (such as Rodgers et al., 2004 andDewitte et al., 2007a for instance), these diagnostics do not show any evidence of ENSO variability rectification through changes in mean state.It is however interesting to note that the observed values (from NCEP reanalyses) distance themselves from the simulated ones and display strong nonlinearity.In addition, most of the ENSIP coupled models (used in Hannachi et al., 2004) tend to underestimate the nonlinearity seen in the NINO3 index, which could be due to significant biases in the simulated mean state and to the limited skills of this first generation of coupled models (Latif et al., 2002).
In the light of Hannachi et al. (2004), this study aimed to examine the role of climate shifts and EEs (with the hypothesis that they emerge from nonlinear processes within the tropics) in controlling ENSO variability.It took advantage of newly designed statistical tools that diagnose the characteristics of the specific distribution law introduced below.As in Hannachi et al. (2004), we made use of CGCM simulations that provided long-term time series of ENSO variability, namely the simulations provided by the World Climate Research Programme Coupled Model Intercomparison Project phase 3 (CMIP3) multi-model data set that was collected for the needs of the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC-AR4).To complement Hannachi et al. (2004)'s approach, indices of the nonlinearity of the tropical Pacific system were used, namely the nonlinear advection within the mixed-layer (also called Nonlinear Dynamical Heating, cf.Timmerman and Jin, 2002) and the skewness of the NINO3 SST index according to An and Jin (2004).Lastly, the focus was on the role of climate shifts (decadal to multi-decadal variability) and EEs in producing departures from normality and the asymmetry of the observed ENSO indices, which extended the study by Hannachi et al. (2004).In summary, the main objective of this paper is to document the statistical properties of ENSO indices from data and different model outputs, in particular those from the new generation of IPCC models, and to corroborate (from a statistical point of view) the recent modelling studies (mentioned above) that emphasize the role of nonlinearity in modulating ENSO properties through variability time scale interactions.
Two main features of the PDF were then explicitly considered in this study: (1) the asymmetry and (2) the "weight" of the tail associated with warm events.In this context, we proposed the use of a specific parametric law as an alternative to Gaussian statistics (a more general framework including Normal distribution) to investigate these features.In particular, the α-stable law was proposed, as it better represents the processes exhibiting the ENSO properties of interest in this paper.In brief, non-Gaussian α-stable laws, also known as "heavy tailed laws" or "infinite variance laws", which were first introduced by Lévy (1924) and then generalized by Mandelbrot (1960), are characterized by four main parameters.The main ones are α and β.The parameter 0<α≤2 allows the "non-Gaussian degree" of the set to be measured.The parameter −1≤β≤1 represents the asymmetry of the law which matches the skewness of Gaussian statistics.Such a law has been used in previous studies to address issues related to financial time series analysis.Most of the time, the reason given for switching from Gaussian to αstable statistics is the desire to take into account "outliers" or EEs whose presence in the series leads to empirical variance bursts and weighs the distribution tails.Moreover, a salient feature of this particular law is that among infinite variance distributions, only stable distributions can be the limiting distribution of sums of independent identically distributed (iid) random variables.In other words, this characteristic of stable distribution can be regarded as the equivalent of the central limit theorem in the Gaussian framework.In addition, estimating α will document and reliably quantify the presence of EEs.This parameter can be viewed as a "proxy" for high order statistical moments (higher than the 4th order momentkurtosis-studied in Hannachi et al., 2004).
In this paper, based on robust statistical tests, we begin by demonstrating that ENSO can be accounted for by non-Gaussian statistics and non-stationary processes dominated by time-mean state and empirical variance shifts.We then hypothesise that these ruptures manifesting as abrupt switches from a cool to warm (warm to cool) ocean background tend to enhance (diminish) feedback processes allowing the burst of EEs.From models of various complexities, the nonlinearity is diagnosed and analysed along with the αstable character of relevant parameters of climate variability in the tropical Pacific (SST and thermocline depth anomalies).The results indicate that the models having the most consistent relationship between changes in mean state and nonlinearity are generally the ones exhibiting the largest deviation from Gaussianity in concordance with greater skill in accounting for EEs.
The paper is organized as follows: Sect. 2 presents observations and model outputs that were used.The α-stable law and the statistical methods used to diagnose the deviation from a Gaussian distribution of the series, as well as the so-called "naïve" statistical ENSO models that are proposed for interpreting the results are also presented in this section.Section 3 presents the results of the statistical analyses on the data and the models.In the light of the results from the dynamic model simulations, Sect. 4 proposes a definition for a model's skill in simulating EEs based on the comparison with the "naïve" statistical models.Section 5 is a discussion followed by concluding remarks.

Observations
In situ and reconstructed data were first used in order to validate our statistical tests and analyse the SSTA patterns with regard to the characteristics of the tropical Pacific mean state and statistics.The monthly SSTAs (referenced to the mean seasonal cycle) in the tropical Pacific region (29   E-60 • W) from the Kaplan optimal analysis of the MOHSST5 data set were calculated for the period January 1870-November 2007 (Kaplan et al., 1998).Note that because of the limited reliability of the reconstructed data (essentially due to a lack of data, see the time series in Fig. 2a), the first fifteen years (1855-1869) were not taken into account.Hereafter, we will refer to this data set as K98.K98 has been used extensively over the past years, for example to assess El Niño modelling forecasts (Chen et al., 2004) or to validate other long reconstructed fields (Rayner et al., 2003).

Model outputs
Analyses were also performed on model outputs.Models of different complexity were used.First of all, we used two intermediate ocean-atmosphere coupled models of the tropical Pacific; the so-called Zebiak and Cane model (Zebiak and Cane, 1987) and the model by Dewitte (2000), hereafter respectively referred to as the "ZC" and "LODCA" models.They are based on similar physics, namely shallow waters for both components, the ocean component including either one baroclinic mode (for ZC) or three baroclinic modes (for LODCA).Whereas the first one was used as a reference for the comparison of the full-physics models, the second one was only used as tool to diagnose the nonlinearities in the full-physics models, as will be explained later.ZC was integrated for 1200 years, with only the 1000 years after the 200year spin-up being analysed.This model has been used extensively for ENSO studies (see Karspeck et al., 2004, among others) because it comprehends the basic dynamics of ENSO.It also simulates an irregular ENSO cycle with chaotic behaviour (Tziperman et al., 1994), which resembles the observations (Karspeck et al., 2004).On the other hand, in a coupled mode, LODCA simulates a quasi-regular ENSO cycle but is more realistic in simulating ocean surface variability due to the consideration of the higher-order baroclinic modes (see Dewitte, 2000 for details).For this reason, LODCA was used in a forced configuration to assess to what extent the variability of the CGCMs (see below) could be explained by equatorial wave dynamics and to infer Nonlinear Dynamical Heating (NDH, see An and Jin (2004) for a definition), which is difficult to infer from direct model outputs.LODCA was therefore forced by monthly wind stress anomalies as derived from the CGCMs after being tuned with the climatology and wave parameters as derived from the CGCM outputs.Such methodology was used in Dewitte et al. (2007a) and the reader is invited to refer to this study for more details.For all the CGCMs used in this study, LODCA was able to simulate a NINO3 index that correlated at the 75% level at least with the NINO3 index inferred from the direct CGCM outputs.This indicates that the ENSO variability in the CGCMs can to a large extent be accounted for by equatorial processes.NDH was then diagnosed from the LODCA outputs for all the CGCMs.
The CGCMs used in the study came from the so-called IPCC data base (see Table 1).The pre-industrial control experiment, in which the concentration of greenhouse gases is fixed at estimates from 1850, was chosen in order to evaluate the performance of the models under past/present climate conditions for two main reasons: (1) this experiment is the one that provides the longest time series and thus has the greatest statistical confidence, (2) the fixed external forcing for long time series makes the analyses at interannual to decadal scales much easier to conduct (no need to remove the trend as in the "20c3m" climate of the 20th century experiment for instance).Monthly outputs were used since we were focusing on low-frequency mechanisms.
As a reference and for comparison with the CMIP3 models, the SODA ocean Reanalysis (Carton and Giese, 2008) was used (1.4.2.version) although it only spans the period 1958-2001.
For both observations and model outputs, the average SSTAs in the NINO3 area, the NINO1.2area (90 • W-70 • W; 0 • -10 • S) and the average SSTAs over the tropical Pacific (120 • E-60 • W; 29 • S-29 • N for K98 and 130 • E-80 • W; 10 • N-10 • S for ZC/LODCA) were retained as ENSO proxies to perform our tests; monthly climatology was removed from the SST at each grid point to derive anomalies.

Methods
A number of statistical tests were used in this study, in order to: (1) characterize and quantify the deviation from a Gaussian law and show some evidence of the α-stable character of the ENSO indices; (2) detect climate shifts in the time series.These tests are presented below.
When assessing model performance, it is common to think in terms of relative skill, or skill compared to some reference run.We have chosen four empirical strategies, which we called our "naïve models".They are also presented in this section.

Statistical tests
We shall begin by giving a short mathematical description of α-stable processes.As previously mentioned, stable distributions were first characterized by Lévy (1924) in a study of normalized sums of independent and identically distributed terms.The distribution F X of a random variable X is said to be stable if the distribution F S n of the random variable, where the X j are independent copies of X, is such that there is a n >0 and b n such that for every real number x, F S n (x)=F X (a n x + b n ).Standard references for the theory of stable distributions are Gnedenko and Kolmogorov (1954;Chapter VII), Feller (1966;Chapters VI, IX and XVII), Samorodnitsky and Taqqu (1994) and Zolotarev (1986).The coefficients a n above are necessarily of the form a n =n 1/α , 0<α≤2, with α being called the characteristic exponent or index of stability.A parameterization of all stable distributions in term of their characteristic functions ϕ is well known, see Gnedenko and Kolmogorov (1954, Sect. 34).It may be written: (2) where, and In Eq. ( 2) the parameter α is the index of stability introduced above, −1≤β≤1 is a measure of skewness, γ >0 is a scale parameter, and the real δ is a location parameter.When β=δ=0, X is said to be symmetric α-stable (SαS which means that X and −X have the same distribution) and its characteristic function takes the particular simple form: Remark: When α=2, the characteristic function (2) becomes ϕ(t)= exp −γ 2 t 2 + iδ t .This is the characteristic function of a Gaussian random variable with mean δ and variance 2γ 2 .Note that in this case, the value of β is not specified since β tan π=0.However one typically associates the Gaussian distribution with the choice β=0.Then the parameters α, β, γ and δ are unique.
All stable distributions are absolutely continuous, unimodal with an eventually bell-shaped density function.However the density is known in closed form in three cases only: the normal distribution corresponding to α=2 with the PDFf ( , m=δ, σ = √ 2γ ; the Cauchy distribution corresponding to α=1 and β=0 (with the PDF: ) and the reciprocal of a χ 2 1 variable corresponding α=1/2, β=1, γ =1, and δ=0 (with the PDF: for x>0, x 3/2 ).For α<2, these distributions are heavytailed.The heavy tails are asymptotically Pareto-like, which means that for α< 2: lim The non-normal stable distributions have been given less attention than the normal distribution probably because the normal distribution is the only stable distribution which has a finite variance.Among infinite-variance distributions, the non-normal stable distributions play an important role, not only because of their closure properties under convolution, but also because only a stable distribution can be the limiting distribution of sums of independent, identically distributed random variables.In this paper we implicitly consider that the SST is the sum of many small terms for which the only possible limit is a stable distribution.Then the problem is to answer the first question: is the variance finite, leading to a normal distribution, or asymptotically infinite, leading to a stable non Gaussian distribution?Mandelbrot (1963), Granger and Orr (1972) and others such as Nikias and Shao (1995) or Nolan (1999) have proposed a number of graphical procedures in order to choose between a Gaussian distribution and another non-Gaussian stable distribution.Of course, a more general and difficult problem would be to test whether a set of data comes from a stable distribution or a non-stable distribution.For the reasons given above, we assume that the distribution is stable and look at procedures for distinguishing non-Gaussian stable distribution.
First from the inspection of the PDF, one can visually infer the deviation from Gaussianity (hereafter referred to as GT1 for Graphical Test 1, see Fig. 2b).
Another method of analysing the infinite variance feature is to plot the sample variance estimate S 2 n based on the first n observations, against n, i.e. where, where N is the number of points in the time series.If all the X j come from the same distribution, then S 2 n should converge to a finite value, if the population distribution F (x) has a finite variance.Otherwise, S 2 n will diverge.We call it: convergence variance test GT2.Note that the non-convergence of S 2 n does not imply infinite variance, if the hypothetical range of possibilities is expanded to include non-stationary series, with the population variance increasing over the time for instance.
Another test, initiated by Mandelbrot (1963) and called the log-tail test or GT3 hereafter in this paper, is to plot the estimate of log P[X>T] against log T where X is the random variable being estimated.This test examines the shape of the tails of the empirical cumulative distribution function and provides information on the behaviour of the distribution for high temperature T .If the true distribution is stable, with the characteristic exponent α, Eq. ( 4) suggests that the plot should be a straight line with a slope −α.Basically, we plot: where 1 |xj |>T =1 if |x j |>T and 0 otherwise, against log T for T >0.If the plot is linear, it is a strong indication that a stable distribution will provide an excellent fit to the available data.All these tests are only judgmental visual inspections of a graph.They are not precise enough to infer real values of stable parameters.The above tests make it possible to highlight a small but significant deviation from a Gaussian distribution of observed ENSO indices as already noticed from the inspection of Fig. 1, namely the heavy tails of the NINO3 SSTA distribution.The empirical variance test (GT2) does not stabilize (see Fig. 2c, GT2 panel).GT3 (see Fig. 2d, GT3 panel) provides an estimate for α of around 1.80 for the NINO3 SST index of K98.Nevertheless this method is quite imprecise for estimating the index of stability.Several methods have been proposed in order to estimate the parameters of a stable law: graphical methods, quantile methods, maximum likelihood ratio methods for example.We used a regressiontype method (TL1A for Telecom Lille I Algorithm) initiated by Koutrouvelis (1980) and classically used by practitioners mainly to infer the value of α and mainly because the amount of computation involved is minimal.We briefly describe the method below.First it is easy to see that Eq. ( 2) implies: We denote by φ N the sample characteristic function which is obtained from the random sample x 1 , x 2 , . . ., x N by Equation ( 5) depends only on α and γ and suggests estimating α and γ by regressing y=log(-log|φ n (T)| 2 ) on w=log |T| in the model: where (t k ; k=1, 2, . . ., K) is an appropriate set of real numbers, µ=log(2γ α ), and ε k denotes an error term.
Once estimates of α and γ have been obtained and α and γ have been fixed to these values so that they are no longer unknown, estimates of β and δ may be obtained, for α =1 using the following equation: Then, we can estimate the parameters by regressing u on sgn(u)|u| α in the model: where η l denotes an error term and (u l ; l=1, 2, . . ., L) is an appropriate set of real numbers.Then, the set of four parameters governing the stable distribution obtained by the previously described two-step procedure is refined in the next step by introducing certain "standardization" to the data and by appropriately choosing the points t k and u l .This regression method also allows the confidence interval around the estimated values of the stable distribution parameters to be derived, as we do have asymptotic variances of the estimators and can thus put confidence limits on the parameters.They are given in Tables 3 and 5.

Shift Detection Test (SDT)
In order to detect shifts in the time series, the method described in Potter (1981) was used.It is based on a bivariate test developed by Maronna and Yohai (1978).The main improvement over other well known tests is the introduction of another correlated series, assumed to be unchanged.Unlike earlier procedures for detecting a shift in mean out of an independent time series, this method is statistically rigorous and provides estimates of the time and amount of change in the mean.For more details, readers are invited to refer to Potter (1981) and to Appendix A. Performing this test on the empirical variance series in order to detect ruptures in SST variability also provides significant information on statistical characteristics.On the one hand, we can distinguish different regimes, in terms of variability features, within a time set.On the other hand, we can clearly identify EEs (responsible for isolated bursts in the empirical variance time series).

"Naïve" statistical models
In order to assess model performance and the sensitivity of ENSO statistics to changes in mean state and EE occurrence, simple statistical models were proposed, based on naïve strategies.Theoretical ENSO series (hereafter referred to as TGS for Theoretical Generated Series) were therefore generated that built upon the aforementioned properties of ENSO time series distribution.
-First of all, α-stable sets were considered.Chambers et al. (1976) developed an algorithm allowing the simulation of a symmetrical α-stable set.It was widened to general α-stable sets by d' Estampes (2003).The generation method is described in Appendix B. The objective was to simulate a set which follows coherent stable statistics with parameters related to in situ data.We chose α=1.80 and β=1 to match the estimates from the K98 data (see above).
-Secondly, a classical Gaussian process was considered and an associated time series was generated.The aim was to simulate a typical Gaussian process which could also be seen as a "pseudo stable" process with α=2.The Gaussian TGS parameters were chosen to fit with those of the K98 NINO3 index, i.e. 0.0040 • C for the mean and 0.8082 • C for the standard deviation.
-Thirdly, a statistical model that highlights the influence of climate shifts on ENSO statistics was proposed.It considered a non-stationary Gaussian process characterized by a threshold in its mean and standard deviation.These ruptures in the process parameters were concomitant with climate shifts evidenced in K98 (1903 and 1976, see Appendix A and Fig. A1) with the characteristics of the observed shifts (mean difference between one period and the next and standard deviation of each period being imposed) being prescribed in this TGS.
-Finally, we combined the above to generate a nonstationary stable set.The aim of such a generation was to simulate a process able to rectify its high-frequency variability (i.e.EE probability occurrence) according to its low-frequency modulation (changes in mean state).
The generation was performed in the same way as that of the previous TGS.However, warm periods were characterized by α-stable statistics consistent with the observed inter-shift periods whereas cold period statistics remained Gaussian (α=2).This TGS can also be viewed as a purely stable process whose intrinsic main parameters (i.e.α and β) experienced a low frequency modulation.
To compare statistical distributions of TGS against observations and model outputs, the quantile-quantile (percentile in that case) plots were used according to Hannachi et al. (2004).The q-q plot approach and their relevance in evaluating model performance are discussed in Hannachi et al. (2004) and Hannachi (2006).In brief, we plot the percentiles of the data or models outputs versus the percentiles of the TGS, to assess if the 2 series come from the same statistical distribution (the q-q plot is the bisector of the plot in

α-stable character of ENSO in data and models
This section presents the results of the estimation of the α and β parameters and documents the relationships between changes in mean state (as revealed by the shift detection method) and ENSO statistics.Regions of strong stability are found around the eastern edge of the warm pool (180 • W) and along the equator in the far eastern Pacific.The statistical tests described in Sect. 2 and applied to the SSTA time series in these regions corroborate the significant deviation from Gaussianity (not shown).Another interesting feature is evidenced in the map of β (Fig. 3 lower panel) which exhibits a zonal see-saw pattern with positive values in the NINO3 region and negative values in the western Pacific and off the equatorial wave guide.The map of β agrees strongly with the pattern for SSTA skewness (see the contours of SST skewness1 overplotted in Fig. 3b).

Observations: K98 data
These statistics are sensitive to the period under investigation.In order to illustrate such sensitivity, the most robust shifts in the time series at each grid point are estimated according to the SDT method (see Sect. 2.2.2. and Appendix A).A minimum inter-shift period of ten years was arbitrarily chosen.Note that the dates of shifts that were detected are quite consistent while performing the SDT method on SSTA or on empirical variance SSTA series.Actually, mean shifts are often followed by shifts in variance which agrees strongly with Sardeshmukh et al. (2000).The statistical properties of the detected inter-shift periods, assumed to be stationary, were then investigated.We only retained the most relevant shifts, i.e. those with a statistical test significance level greater than 7.90, which corresponds to a 90% significance level for a 100-occurrences set (see Potter (1981) and Appendix A).The results are presented in Fig. 4 and summarized in Table 3.The values of (α, β) in Table 3 were inferred from the TL1A algorithm.Consistently with former studies, the detected shifts took place around 1903 (Karspeck et al., 2004), 1976(Guilderson and Schrag, 1998) and 1998 (Overland et al., 2008) (cf.Table 3).Nonetheless, we did not take into account the last detected shift (1998) 50  as the following period was too short relative to our criterion.The SDT also provided a value for the change in mean state (Fig. 4, lower panel) which was also consistent, though slightly overestimated, with earlier works previously mentioned.Over the inter-shift periods, the results of Fig. 4 indicate that the ENSO statistics experienced significant changes, consistent with the study by Karspeck et al. (2004).In particular, warm periods were characterised by stronger asymmetry and a greater deviation from Gaussianity (smaller α and β around 1 Fig. 4a, b, e, f) whereas the cool period exhibited a Gaussian symmetrical pattern on average over the tropical Pacific (α≈2 and β around 0, Fig. 4c and d, see also Table 3).Note that over the period 1998-2007, there was a significant reduction in the stable nature of the NINO3 SST index and its asymmetry (α=2 and β=0.04, estimated using TL1A), which is consistent with the above.
Figure 4 also highlights spatial pattern differences for the mean SST change for the 1903 and 1976 shifts with the 1976-shift mean SST change having a El Niño like structure (Fig. 4h) whereas the 1903-shift mean SST change corresponded to a reduction in the mean zonal SST gradient near the eastern edge of the warm pool (Fig. 4g).Such changes are likely to be associated with distinct impacts on ENSO dynamics and are consistent with the changes in ENSO statistics.In particular, a decrease in mean zonal SST gradients as observed for the 1976-shifts is concurrent with a flattening mean thermocline (Moon et al., 2004), which impacts ENSO towards larger amplitude modulation (Dewitte et al., 2007a).This also favours the amplification of ENSO nonlinearity (Timmerman and Jin, 2002) and thereby is likely to modify the statistical characteristics considered in this study.We will come back to this issue in the last section of the paper.
Additional tests were performed on the K98 data set in order to investigate further the sensitivity of the statistics to the ruptures in tropical Pacific mean background.
For instance, TL1A is applied to the SSTA series over the 1950-2007 period from which the 1976 shift had been removed.This was achieved by keeping the same SST sets prior to the detected date and removing from the data after the shift, the mean SST amplitude change measured by the SDT at each grid point (for example 0.3 • C for the averaged NINO3 index).The results are presented in Fig. 5.They indicate that removing the 1976 shift from the series led to a reduction in stability in the eastern tropical Pacific.In particular, the α parameter increased up to 1.9 instead of 1.7 for the "unchanged" series in the region of the Humboldt Current System.Actually, an α parameter close to the Gaussian value of 2 was observed over most of the rest of the domain except for the eastern side of the warm pool (not shown).In summary, the above results obtained from observed data suggest that deviation from Gaussianity and asymmetry of SSTA distribution are associated with time-mean state changes.This is consistent with recent studies based on physical model and observations which put forward that changes in ENSO properties are linked to changes in mean state through the nonlinearities of the tropical Pacific (An, 2004;Dewitte et al., 2007a).In the light of the above results with the K98 data set, the following section investigates such an issue from various model outputs.

Models
Two model types are considered below: 1) An intermediate complexity model, the ZC model, that has been widely used for ENSO studies and 2) full-physics coupled general circulation models (CGCM) whose outputs are available to the scientific community within the Program for Climate Model Diagnosis and Intercomparison (see Sect. 2.1.2. and Table 1).

The ZC model
The ZC model was run over 1200 years and the last 1000 years were analysed.Although based on linear dynamics for 52    the ocean and atmospheric parts, the nonlinearity of the system is considered through the ocean thermodynamics and a moisture feedback process for the atmosphere.Due to its formulation, NDH, can easily be diagnosed from the model outputs.Three model fields are analysed below: SST, Thermocline Depth Anomaly (TDA), and NDH.Monthly anomalies were calculated relative to the mean seasonal cycle calculated over the 1000-year period.Similar statistical tests to those applied to the K98 data set were performed on the model outputs.The SDT detected 20 main shifts over the 1000-year period (see Fig. 6 and Table 4 which summarize the results for the NINO1.2index for the first 200 years of simulation, for clarity only).The results are robust since the detected shift dates are comparable for both the mean and the empirical variance and for each proxy: SSTA, TDA and NDH.Nonetheless, the differences in the results of the SDT are smaller between NDH and TDA than between NDH and SSTA (especially in the eastern part of the basin) supposedly because SSTA is influenced by both the direct ENSO asymmetric forcing and the NDH, whereas TDA results directly from the linear response of the wind forcing.
In order to infer the statistical properties of the inter-shift periods, (α, β) were estimated for the "cool" and "warm" periods.The latter periods were detected by applying the SDT to the NINO1.2index.Fig. 6 presents the 30-years running mean of NINO1.2 averaged SSTAs; vertical lines indicate the shift dates, as estimated by the test.Note that the SDT was performed on both the raw and the filtered series and led to similar results.Negative shifts were followed by a cooler period (characterized by blue shading in Fig. 6) whereas positive shifts led to a warmer tropical Pacific (see overlapping red shading in Fig. 6).In a second step, TL1A was performed on each inter-shift period.Statistics for each warm (cool) period were then averaged to derive a mean value for α, β and NDH characterizing a warm (cool) background.The results are displayed in Fig. 7.In order to highlight the deviation of α from 2, the value of exp(4-α) instead of α was considered (Fig. 7a, b, e, f).Blue shading is synonymous of Gaussianity whereas red accounts for non-Gaussian stable statistics.The average value of α over the domain (130 • E-80 • W; 10 • N-10 • S) is indicated at the top of each map.The results indicate that the simulated SSTAs and TDAs tend to be more non-Gaussian during warm periods than during cold periods.The contrast is even more striking when looking at the differences between the spatial patterns for the different periods, particularly visible on the β maps (Fig. 7c, d, g and h).Indeed a warm ocean background emphasizes the contrast in asymmetry between the western/eastern Pacific, particularly clear on the SST field (Fig. 7g and h), whereas asymmetry was less pronounced during cool periods (extended symmetrical patterns, where β≈0).In contrast with K98, the spatial variability of α for the ZC model remains difficult to interpret.Nevertheless, periods for which tropical Pacific mean state is warm were characterized by higher values of NDH mean than cold periods which were associated with a lower NDH mean (Fig. 7i and j).
These results can be interpreted in the light of recent model studies which attribute to nonlinear advection a role in rectifying ENSO variability (An and Jin, 2004;Timmermann et al., 2003;Jin et al., 2003).For instance, Timmermann et al. (2002) suggested that El Niño bursting was associated with an increased NDH whereas La Niña events had a lower value of NDH.Consistent with our results, warm periods regarded as "nonlinear active" manifested more stable statistics in terms of deviation from Gaussianity and asymmetry, whereas cool periods experienced more Gaussian statistics.It is hypothesized that during the periods when NDH contributes to enhanced ENSO amplitude, the growth of warm EEs is favoured.Conversely, during cool periods, NDH (asymmetry) is reduced along with the occurrence of EEs.Still, isolated extreme warm events can take place during cool "pseudo linear" periods but with fewer occurrences than over warm periods.We will examine this hypothesis further in the discussion section.

IPCC models
The IPCC data base offers the opportunity to investigate how more complex models behave regarding the relationship between ENSO statistics and changes in mean state.Note that, despite being full-physics, the IPCC models exhibit numerous biases, especially in the mean states which drastically impacts ENSO variability.Knowing the characteristics of these biases may help interpret the ENSO statistics.For instance, in some models, the ENSO dynamic is dominated by thermocline feedback processes, which overestimates the control of SST by the thermocline depth anomaly and thus ENSO asymmetry in the eastern tropical Pacific.Some others favour the dominance of the zonal advective feedback generally leading to faster and more regular ENSO variability (Van Oldenborgh et al., 2005;Guilyardi, 2006;Dewitte et al., 2007b).Despite these biases, interestingly, there was a clear relationship in all these models between changes in mean state and nonlinearities (as measured by NDH).As a demonstration of this latter statement, a singular value decomposition (SVD) analysis was carried out between the 11-year running skewness and mean of SST over the tropical Pacific domain (10 • S-10 • N) for all models, following similar methodology to An (2004) (his Fig. 1).The results of the SVD analysis are reported in Table 5 for the first dominant mode, which indicates a significant percentage of covariance between skewness and mean SST for all the models, along with a high correlation between associated time series.Differences between models were mostly 54    1. Nº25 is for SODA.Red overlapped curve is the best fitted power law to all IPCC models (Group 1), blue one is the best fitted power law to only "good" (see text) models (Group 2)."Bad" models are represented by black filled triangles.Each number represents a model, as referenced in Table 1.No. 25 is for SODA.Red overlapped curve is the best fitted power law to all IPCC models (Group 1), blue one is the best fitted power law to only "good" (see text) models (Group 2)."Bad" models are represented by black filled triangles.
for all models (NINO3 SST index) and reported in Table 5 along with the number of detected shifts.The results indicated that there was a great diversity of behaviour regarding α-stability and asymmetry.This is in contrast with the results of the above SVD analysis.In order to visualize the differences between models, Fig. 8 is presented, which displays various scatter plots representing the relationships between the different orders of ENSO statistical moments and a proxy of nonlinearity, namely the root mean square (rms) of NDH (nonlinear advection) as diagnosed from LODCA (see Sect. 2.1.2.) (NDH was normalized by the rms of the NINO3 index).Here, the shift number (brought back to a 100-year period) quantifies the variability of the abrupt changes in ocean mean state, and is thus equivalent to the 1 st order statistical moment.β represents the asymmetry of ENSO and consequently can be assimilated to a 3rd order statistical moment; whereas α gives information on EEs or in other words on the abundance of rare values, i.e. higher order statistical moments.For an α-stable distribution, remind that for r≥α, E(|X| r )=+∞ A detailed examination of Fig. 8 indicates that there is a highly nonlinear relationship between statistical moments.Note for instance, the significant number of models being Gaussian ( α=0), while also exhibiting a marked asymme-try (NDH variability) (Fig. 8d and e).In order to quantify the nonlinear dependency between statistical moments, power laws (i.e.f(x)=a+b.xµ Stanley, 1995) were used and fitted to the model ensemble for the various scattered plots in Fig. 8.The power laws were fitted by minimizing the rms residuals for two different groups of models: (1) for all the IPCC models listed in Table 5, hereafter referred to as Group 1 (red curves in Fig. 8) and for a selected smaller group, hereafter referred to as Group 2 (blue curves in Fig. 8).The latter is composed of the most 'realistic' models according to recent works (Van Oldenborgh et al., 2005;Guilyardi, 2006;Capotondi et al., 2006;Belmadani et al., 2009).The models that were excluded from this group (BCCR-BCM2-0, CCCMA-CGCM3.1,GISS-AOM, GISS-E-H) were shown to simulate a biased ENSO variability (see aforementioned studies for more details).Note that these models are not considered either in most analyses of the multi-model studies by van Oldenborg et al. (2005) and Guilyardi (2006).Three other models (IPSL-CM4, NCAR-CCSM3.0 and CSIRO-MK3.0)were arbitrarily removed from this group as they exhibited "exotic" statistical behaviour (namely no EE occurrence despite a realistic/marked positive asymmetry).The "eliminated" models are represented by black triangles in Fig. 8.  2006)'s methodology.The four other rightmost columns provide the results of the SVD between change in mean state and skewness based on SSTA (a 11-year running window is used), namely the percentage of covariance of the first SVD mode, the percentage of the variance for the mean and skewness for the first SVD mode and the correlation between the associated time series.Results of SODA are given as a reference for some parameters (bottom line).Non-Gaussian stable models (α<2) are written in bold.α and β are estimated using TL1A.The bottom panels in Fig. 8 corroborate the previous results from the ZC model, namely the existence in the IPCC models of a power law-type relationship between nonlinear activity and EE occurrence (Fig. 8e) as well as a quasi-linear relationship between NDH and asymmetry (Fig. 8f), consistent with earlier studies (An and Jin, 2004;An et al., 2005).Indeed, these two curves fitted the power laws relatively well, having a coefficient µ equal to 0.38 and 0.05, respectively.These curve fits were consistent as the residual of the rms was quite low, respectively 0.09 and 0.01.The residuals fell respectively to 0.07 and 0.005 when fitted to the models of Group 2. In that case, µ has values equal to 0.25 and 0.03) for the above-mentioned scattered plots (Fig. 8e and f).Nonetheless, despite the low value of residuals in the power law regression method, one can note that for the relationship between EEs and NDH (Fig. 8e), no satisfactory visual fitting was obtained (due to the significant spread of the models), suggesting that NDH cannot fully account for the dynamics of the EE occurrences.Similarly, an absence of any relationship (according to the power law) is evidenced between the number of shift and NDH, suggesting again that nonlinearities associated with nonlinear advection cannot alone explain how climate shifts are triggered (Fig. 8a, µ=0.21,Residual=8.70 for Group 2).In the same way, although displaying similar exponents in the power law fit, no clear relationship was observed between ENSO asymmetry and the number of shifts (Fig. 8c, µ=0.19,Residual=3.32 for Group 2) , and neither between the number of shifts and EE occurrence (Fig. 8b,µ=0.20,Residual=8.75for Group 2).On the other hand, EE occurrence and ENSO asymmetry are clearly related as there was a strong agreement with the fit (Fig. 8d, µ=0.30,Residual=1.29×10 −04 for Group 2).
The fact that the power law can be used to fit the relationship between statistical moments emphasizes the complex scaling relationships associated with the ENSO modulation from low-frequency mean state change to EE occurrence.Interestingly, NDH relates to intermediate statistical moments (i.e.variability and asymmetry) with low values of residuals (which is consistent with the aforementioned studies).However it does not seem to be the main nonlinear term governing interaction between "extreme" statistical moments (1st and high orders, i.e. slowly varying mean state and triggering of EEs).This clearly highlights the variability time scale interactions associated with ENSO that are evidenced here through the various statistical moments.It also suggests that other sources of nonlinearity than NDH are involved in www.nonlin-processes-geophys.net/16/453/2009/ Nonlin.Processes Geophys., 16,2009 the processes leading to climate shifts and EE occurrence.Due to the complexity of the tropical Pacific system, these sources are numerous (see An, 2009 for a review) and are more or less well represented in the current generation of global coupled models.
To summarize, we find that the IPCC models, although all exhibiting a relationship between asymmetry and slowly varying mean state, behave differently with regard to nonlinearity as measured by α and slowly varying mean state.Only a few models have a value for α that is comparable to observations.Interestingly, these models simulate a relatively large number of climate shifts, consistent with the observations (K98, see Sect. 2.1.1.).This suggests that the nonlinear processes involved in the generation of climate shifts are somewhat different to the ones leading to the rectification of ENSO variability by the slowly varying changes in mean state.The physical processes responsible for the occurrence of climate shifts remain unclear.The results reported here suggest that they may not have a signature in the third statistical moment (asymmetry) of the dynamic fields but may involve higher-order statistical moments and related nonlinearities.
The following section is dedicated to the statistical interpretation of the above results proposing "naïve" models of ENSO, and the definition of a measure of "model skill" for simulating extreme ENSO events based on q-q plot estimates.

Statistical parameterisation of ENSO
In this section, four simple or "naïve" statistical models are proposed (see Sect. 2.2.3.) in order to interpret ENSO statistics with regard to their relationships with changes in mean state and EE occurrence.

Alpha-stable TGS
The first type of model is based on stable TGS (see Appendix B and Sect. 2.2.3.).This model's results indicate that such parameterization leads to an under-estimation of cold events (Fig. 9a lower quadrant) and an over-representation of high frequency events.However, the results from GT2 and GT3 applied to the generated series match with those for in situ data (not shown).The q-q plot presented in the left-hand panels of Fig. 9 highlights the lack of statistical representativeness of a purely stable process (over-representation of extreme warm events certainly due to the too high asymmetry that is prescribed, β=1, see Fig. 9e upper quadrant), notably on the NINO3 index.The too high probability of EE occurrence is certainly due to the too high variability of the TGS, especially on short time scales.Actually, the TGS is essentially based on an initial random set generation, unable to account for the observed persistence (inertia) in large-scale ocean circulation.

Gaussian TGS
The results from GT (not shown) indicate that this simple statistical model cannot account for the ENSO asymmetry.However, q-q plots (cf.Fig. 9b) show that it leads to a better statistical distribution of the generated series than alphastable TGS.Still, warm EEs remain underestimated (Fig. 9b upper quadrant).It is noteworthy that cold events are well represented by this generation method.

Non-stationary TGS
This non-stationary process enhances the representativeness of cold events (see Fig. 9c).Actually q-q plot matches quasiperfectly for negative SSTAs whereas warm events remain under-represented (Fig. 9c upper quadrant).

Non-stationary stable TGS
This process clearly enhances the representation of warm events, but without altering the probability of the occurrence of cold events as a purely stationary stable TGS described in the first part of this section (cf.Fig. 9d).This is particularly true when observing q-q plot of the TGS versus K98.Although there is strong agreement between observations and this TGS, the latter does not statistically fit with the ZC model results.It is believed that this is due to the fact that the ZC model cannot account for all the nonlinear dynamics leading to EEs.
The results of the proposed "naïve" models illustrate the complexity of ENSO with regard to its statistical properties.It indicates that ENSO cannot be accounted for by a single statistical law, or at least by a law whose intrinsic parameters are permanent over time.Unlike a purely α-stable distribution, a Gaussian law fails to represent positive SSTAs while cold events have a perfect Gaussian distribution.Finally, the results of this "naïve" parameterisation of ENSO further suggest that ENSO experiences various types of behaviour, which combines Gaussian distribution for cold symmetrical periods (α=2) and α-stable for warm active periods (α<2).Such a parameterisation (a non-stationary stable process with the main parameter values following the slowly varying Tropical Pacific mean state) could account for the rectified effect leading to the occurrence of EEs.
In the following paragraphs, in the light of the above, we investigate the statistics of the full-physics models based on q-q plot analyses.
To quantify the ability of the IPCC models to represent strong warm episodes, we then propose the following criterion C:    found (2.25 • C in K98).Simply put, T model (q T =99%) represents the minimum SSTA for an El Niño to be considered as an EE in each IPCC model.Graphically, C represents the distance between the q-q plot of the model and the bisector at the T data (q T = 99%) (= 2.25 • C) abscissa (see black vertical arrow on the upper quadrant of Fig. 9d).According to Hannachi (2006)'s methodology, we are able to compute error estimates on quantile calculation and can thus provide confidence intervals on C.
The results of the classification are shown in Table 5.The lower the value for C, the better the model sim-ulates ENSO statistics in term of representativeness of EEs.Negative values for C indicate that the model overestimates strong warm events, or in other words gives too much weight to the positive tail of the distribution.We can also point out that the majority of "good" models according to this criterion (BCCR-BCM2, CSIRO-MK3.5,GFDL CM2 0, INM-CM3.0,MRI.CGCM2.3.2A,UKMO-HadCM3) exhibit stable statistics.
Interestingly, these "good" models (including MIUB ECHO G, MPI ECHAM5) also match models simulating realistic inter-decadal variability as shown in Lin (2007) 5).However, unlike K98, these models tend to simulate more numerous cold EEs than warm EEs due to the negative asymmetry of their NINO3 SST index.

Discussion and conclusions
In this paper, ENSO statistics which are accounted for by α-stable distribution are related to some aspects of ENSO's observed characteristics, namely its modulation, its asymmetry and its tendency to produce EEs.As the PDF associated with the NINO3 SST index deviates significantly from the Gaussian distribution, the heavy tailed α-stable distribution is proposed because it better accounts for the occurrence of EEs.Although it is impossible to have access to its PDF in a closed form (except in particular cases), the distribution is characterized by two main parameters, α and β, that provide meaningful information on ENSO statistics, namely EE abundance and asymmetry, respectively.A shift detection method initially developed by Potter (1981) was also used to diagnose the change in mean state of the tropical Pacific and select warm and cool periods in the time series.The observations were first investigated based on the K98 SST data set.Consistent with former studies (Burgers and Stephenson, 1999;Hannachi et al., 2004;An and Jin, 2004), the results indicate that ENSO has α-stable non-Gaussian features and is asymmetrical.Interestingly, cool and warm periods exhibit different statistical behaviour, with cool periods being more Gaussian and having lower asymmetry than warm periods.A comparable tendency was found in the ZC model.
In particular the ZC model had increased (reduced) nonlinearity quantified through NDH (Timmerman and Jin, 2002) during warm (cool) periods.
The full-physics models of the IPCC data base were then investigated.Interestingly, all the models exhibited a clear relationship between changes in mean state and ENSO asymmetry (skewness), in agreement with observations (An, 2004).They did however have contrasting statistics in terms of their propensity to simulate EEs.In particular, in the light of the results of the SDT applied to the IPCC models (Table 5) and a recent study (Lin, 2007), it was shown that only the models simulating a realistic decadal variability also exhibited marked α-stable statistics.
Naïve statistical models were then proposed to interpret these results.The "naïve" model simulations indicated that ENSO could be accounted for by a non-stationary stable process with the typical exponent of the law experiencing variations that mimic the changes in the tropical Pacific Ocean background.As α contains information on EE abundance but also on the decay rate of the ENSO PDF tail, this corroborated the existence of interactions between statistical moments of ENSO time series.The results of q-q plot (see Hannachi et al., 2004) applied to the IPCC model simulations, and of the comparison between the various quantities studied in this paper (NDH, α, β, number of shifts) through power laws suggest that the interaction between statistical moments (variability time scales) does not solely operate through nonlinear advection or the nonlinearities associated with ENSO asymmetry.It is then likely that other nonlinear processes come into play to explain EE occurrences.Investigating the sources of these nonlinearities is beyond the scope of this paper.At this stage it is interesting to note that, although current measures of ENSO nonlinearities (through either skewness or NDH) have provided meaningful information on the rectification of ENSO variability by changes in mean state (An, 2004;An et al., 2005;Dewitte et al., 2007a), they may not fully account for the complexity of the rectified effect.In the light of the results, we can hypothesize that EE occurrences are part of the feedback loop linking changes in mean state and ENSO asymmetry (An, 2004;Dewitte et al., 2007a).The schematic diagram in Fig. 10 summarizes the parallel that has been made in this paper between ENSO statistical moments and the physical processes involved in the rectification of ENSO variability through changes in mean state.It suggests that higher statistical moments contribute to the rectified effect by controlling the triggering of EEs, and supposedly some feedback between EEs and climate shifts.Non-linear regressions between statistical moments from IPCC model time series give similar exponents, which could suggest ENSO chaotic behaviour.This still requires further study.can have unknown parameters, and which are independent of the unobservable noise u j .The u j are iid normal with mean 0 and variance σ 2 u .Maronna and Yohai called it Model II.Under these assumptions x is an ancillary statistic for tests concerning b j and c.The null hypothesis H 0 is that b j =b; j=1; . . .; N for some unknown b; against the alternative H 1 that b j =b for j ≤j 0 b j =b+d for j >j 0 where b; d and j 0 <N are unknown.
The percentile points are extremely difficult to calculate either analytically or numerically and depend on h.That is why; instead of trying to compute p-values we chose a level far larger than all the published quantiles for these tests and decided to reject H 0 if the computed value of the test statistic exceeded this level.
We performed the test on the NINO3 index computed from K98 using various reference series X: SSTAs in the whole tropical Pacific (120 • E-60 • W; 29 • N-29 • S) and NINO1.2 regions, and a random-generated Gaussian set whose mean and standard deviation are the same as the NINO3 index.Consistent with earlier studies (cf.Karspeck et al., 2004, or Guilderson et al., 1998 among others), a shift in April 1976, estimated amplitude: 0.29 • C was detected in each of these experiments.We also evidenced the cold 1903 and 1998 shifts (Overland et al., 2008).
On the other hand, we performed the test on empirical variance series.An example is given for the NINO1.2.index.This way of running the method is not only able to isolate ruptures in variance (making it possible to identify homogeneous periods in terms of variability) but also to highlight isolated extreme bursts in the empirical variance set (allowing EEs to be identified).Following a dichotomic way of performing the test, we were able to evidence the 1903-1905 cold shift, the 1943 neutral shift (Karspeck et al., 2004) and the 1958, 1983and 1998 EEs (see Fig. 4 and Fig. A1).Nonetheless, the SDT applied to empirical variance sets did not permit the 1976 shift to be located as this rupture took place in the mean rather than in the variance.In this study, we assumed a minimum inter-shift period of 10 years in order to parallel the occurrence of ruptures with decadal to interdecadal variability.

Figure 2 .
Figure 2. Graphical statistical tests used to diagnose the deviation from Gaussianity of time series.a) Full Kaplan Nino3 SST time series.b) Smoothed histogram of the data with a

Fig. 2 .
Fig. 2. Graphical statistical tests used to diagnose the deviation from Gaussianity of time series.(a) Full Kaplan Nino3 SST time series.(b) Smoothed histogram of the data with a Gaussian PDF overlapped (K98 data fitted).(c) Empirical variance of K98 (solid line) and empirical variance of the Gaussian fitted distribution (dash line).(d) Asymptotical test (g(T ) as in Sect.2.2.1) for K98 data (solid line) and for the Gaussian fitted distribution (dash line) (see their description in Sect.2.2.1.).

Figure 3
Figure 3 presents the map of α and β parameters in the tropical Pacific for the K98 data set using the TL1A method.It clearly highlights the non-Gaussian α-stable nature of the SSTAs in most regions of the tropical Pacific since α is lower than 2 over most of the basin (except in the south eastern tropical Pacific, around ∼15 • N in the central-eastern Pacific and the far western equatorial Pacific, see Fig. 3 upper panel).Regions of strong stability are found around the eastern edge of the warm pool (180 • W) and along the equator in the far eastern Pacific.The statistical tests described in Sect. 2 and applied to the SSTA time series in these regions corroborate the significant deviation from Gaussianity (not shown).Another interesting feature is evidenced in the map of β (Fig. 3 lower panel) which exhibits a zonal see-saw pattern with positive values in the NINO3 region and negative values in the western Pacific and off the equatorial wave guide.The map

Figure 3 .
Figure 3. α (upper panel) and β (lower panel) maps of SSTA in the tropical Pacific from K98.The TL1A method was used over the 1870-2007 period.The isotherm 28°C, for the mean SST, is overplotted on the α map to locate the Warm Pool region.On the β map, the contours for skewness (-0.2 and 0.2 iso-contours) are overlapped to highlight the consistency between β and asymmetry.

Fig. 3 .
Fig. 3. α (upper panel) and β (lower panel) maps of SSTA in the tropical Pacific from K98.The TL1A method was used over the 1870-2007 period.The isotherm 28 • C, for the mean SST, is overplotted on the α map to locate the Warm Pool region.On the β map, the contours for skewness (−0.2 and 0.2 iso-contours) are overlapped to highlight the consistency between β and asymmetry as measured by the 3rd statistical moment.

Figure 4 .
Figure 4. α (left) and β (right) maps in the tropical Pacific computed from K98 with TL1A over the periods 1870-1903 (a and b) 1903-1976 (c and d) and 1976-2007 (e and f).Inter-shift periods and the amplitude in °C of the shifts for the 1903 shift (g) and the 1976 shift (h).Bottom panel: NINO3 index (thin line) with its means over the inter-shift periods overlapped (thick horizontal lines).

Fig. 4 .
Fig. 4. α (left) and β (right) maps in the tropical Pacific computed from K98 with TL1A over the periods 1870-1903 (a and b) 1903-1976 (c and d) and 1976-2007 (e and f).Inter-shift periods and the amplitude in • C of the shifts for the 1903 shift (g) and the 1976 shift (h).Bottom panel: NINO3 index (thin line) with its means over the inter-shift periods overlapped (thick horizontal lines).

Figure 5 .
Figure 5. Differences in statistics for the α (upper panel) and β (lower panel) parameters between the 'non-filtered' K98 SST and the 'filtered' (i.e. with the 1976 shift removed, see text) K98 SST.

Figure 6 .
Figure 6.Running mean (30 years window) of the NINO1.2.index.The inter-shift periods are indicated (the shifts were detected from the non-filter series) with a shading (blue shading indicates cool periods, red shading indicates warm periods).

Fig. 6 .
Fig. 6.Running mean (30 years window) of the NINO1.2.index as simulated by the ZC model.The inter-shift periods are indicated (the shifts were detected from the non-filter series) with a shading (blue shading indicates cool periods, red shading indicates warm periods).

Figure 7 .Fig. 7 .
Figure 7. Mean value for warm (left) and cool (right) periods for TDA and SSTA and NDH mean.

Figure 8 .
Figure 8. Scatter plots of the statistical moments and nonlinear terms for the IPCC models.Upper left panel (a.):Numbers of shifts (per 100 years of simulations) vs. nonlinearities; upper right panel (b.): numbers of shift vs. deviation from Gaussian models.Middle left panel (c.): numbers of shift vs. asymmetry; middle right panel (d.): asymmetry vs. deviation from Gaussianity.Bottom left panel (e.): nonlinear terms vs. deviation from Gaussianity; bottom right panel (f.): nonlinear terms vs. deviation from asymmetry.Each number represents a model, as referenced inTable 1. Nº25 is for SODA.Red overlapped curve is the best fitted

Fig. 8 .
Fig. 8. Scatter plots of the statistical moments and nonlinear terms for the IPCC models.Upper left panel (a): Numbers of shifts (per 100 years of simulations) vs. nonlinearities; upper right panel (b): numbers of shift vs. deviation from Gaussian models.Middle left panel (c): numbers of shift vs. asymmetry; middle right panel (d): asymmetry vs. deviation from Gaussianity.Bottom left panel (e): nonlinear terms vs. deviation from Gaussianity; bottom right panel (f): nonlinear terms vs. deviation from asymmetry.Each number represents a model, as referenced in Table1.No. 25 is for SODA.Red overlapped curve is the best fitted power law to all IPCC models (Group 1), blue one is the best fitted power law to only "good" (see text) models (Group 2)."Bad" models are represented by black filled triangles.

)Figure 9 .
Figure9.Quantiles of an asymmetrical alpha stable TGS (with α =1.8 and β = 1) versus quantiles of data / model outputs, (a.); quantiles of a Gaussian TGS (with m =0.0040 and σ = 0.8082) versus quantiles of data / model outputs, (b.); quantiles of an non-stationary Gaussian TGS (with m = -0.0639and σ = 0.7464 on the first homogeneous period, m = -0.1241and σ = 0.7738 on the 2 nd and m = 0.3696 and σ = 0.8338 on the 3 rd ) versus quantiles of data / model outputs, (c.); quantiles of stable, asymmetrical, non-stationary TGS (with α =1.85 and β = 0.99 on the first homogeneous period, α =2 and β = -0.17 on the 2 nd and α =1.66 and β = 0.99 on the 3 rd ) versus quantiles of data / model outputs, (d.).Upper quadrants are associated with warm anomalies whereas lower quadrants are related to cold anomalies.Red colour is associated with the NINO3 index computed from K98, purple colour is associated with the NINO1.2index computed from K98, green colour is associated with the NINO3 index computed from ZC, blue colour is associated with the NINO1.2index computed from ZC.On each panel, the ideal q-q plot (NINO3 quantiles from K98 vs. NINO3 quantiles from K98) is indicated in black solid line.The black arrow represented on the upper quadrant of (d) indicates the measurement of the C criterion (see Section 4.4.).

Figure 10 .
Figure 10.Schematic of the mechanism of interaction between ENSO time scales variability and change in mean state and its relationship with ENSO statistics.

Fig. 10 .
Fig. 10.Schematic of the mechanism of interaction between ENSO time scales variability and change in mean state and its relationship with ENSO statistics. 58

Figure A1 .
Figure A1.15 years running mean of the NINO3 index (upper panel) and empirical variance of the NINO1.2.Index (lower panel) from K98. Mean shifts and EE detected by the SDT are reported on the plot.

Fig
Fig. A1.15-years running mean of the NINO3 index (upper panel) and empirical variance of the NINO1.2.Index (lower panel) from K98. Mean shifts and EE detected by the SDT are reported on the plot.

Table 1 .
Description of the CGCM simulations considered in this study.

Table 2 .
Description of the mathematical tests used in this study.

Table 3 .
Main detected mean shifts from K98 (tests performed with K98 NINO3 SSTA index and Gaussian reference set) and estimation of stable parameters on each inter-shift periods with two different methods (Nolan on line estimation and TL1A in bold).
that case) or at least to evaluate if the statistical properties of the sets are close.Table2summarizes and briefly describes the statistical tests used in this study.

Table 4 .
Detected mean and variance shifts for SSTA, TDA and NDH over 200 years of NINO1.2 index from a 1000-year ZC model simulation (mean are indicated in black whereas variance is written in italic).

TDA [m] SSTA [ • C] NDHA [ • C/month]
Variance Mean, Variance Mean, Variance Mean, Variance Mean, Variance Mean, Variance and/or mean, with some models having a low explained variance for the skewness (such as MIROC3.2-MEDRES)andothers having a high explained variance for the mean (such as IPSL-CM4).These results convey the fact that the slowly varying mean state is related to the nonlinearity of the equatorial Pacific system in all the models.Below, we investigate how this translates to the ENSO statistics in terms of stability and asymmetry.Applying the tests described in Sect.2.2.1.(TL1A),the α and β parameters were calculated www.nonlin-processes-geophys.net/16/453/2009/ Nonlin.Processes Geophys.,16,2009

Table 5 .
Statistical properties of the ENSO variability as simulated by the IPCC models.The last column presents the values of the criterion C and the corresponding standard errors computed followingHannachi ( Quantiles of an asymmetrical alpha stable TGS (with α=1.8 and β=1) versus quantiles of data/model outputs, (a); quantiles of a Gaussian TGS (with m=0.0040 and σ =0.8082) versus quantiles of data/model outputs, (b); quantiles of an non-stationary Gaussian TGS (with m=−0.0639 and σ =0.7464 on the first homogeneous period, m=−0.1241 and σ =0.7738 on the 2nd and m=0.3696 and σ =0.8338 on the 3rd) versus quantiles of data/model outputs, (c); quantiles of stable, asymmetrical, non-stationary TGS (with α=1.85 and β=0.99 on the first homogeneous period, α=2 and β=−0.17 on the 2nd and α=1.66 and β=0.99 on the 3rd) versus quantiles of data/model outputs, (d).Upper quadrants are associated with warm anomalies whereas lower quadrants are related to cold anomalies.Red colour is associated with the NINO3 index computed from K98, purple colour with the NINO1.2index computed from K98, green colour with the NINO3 index computed from ZC and blue colour with the NINO1.2index computed from ZC.On each panel, the ideal q-q plot (NINO3 quantiles from K98 vs. NINO3 quantiles from K98) is indicated in black solid line.The black arrow represented on the upper quadrant of (d) indicates the measurement of the C criterion (see Sect. 4.4.).