Decadal ENSO amplitude modulations: a nonlinear paradigm

Global and Planetary Change 37 (2003) 135 – 156 www.elsevier.com/locate/gloplacha

Decadal ENSO amplitude modulations: a nonlinear paradigm Axel Timmermann Institut fu¨r Meereskunde, Universita¨t Kiel, Du¨sternbrooker Weg 20, D-24105 Kiel, Germany Received 20 June 2001; accepted 6 July 2002

Abstract This paper investigates the mechanism which generates decadal modulations in the amplitude of the El Nin˜o-Southern Oscillation phenomenon (ENSO). Our analysis is based on a multicentury present-day climate simulation performed with an ENSO-resolving Coupled General Circulation Model (CGCM). In consistency with observations, it is found that ENSO variance undergoes changes with a time scale of about 10 – 20 years. This decadal beat is closely linked to the second dominant pattern of tropical (sub)surface temperature variability. The dipole-like characteristic of this mode is generated mainly by the interplay of horizontal, vertical advection and mixing. We suggest a nonlinear mechanism, which is capable of generating decadal tropical climate anomalies as well as decadal ENSO amplitude modulations (DEAMs) without invoking extratropical dynamics. This mechanism is based on the idea of homoclinic orbits. This new paradigm is validated using a low-dimensional ENSO model that is derived empirically from the CGCM simulation. D 2003 Elsevier Science B.V. All rights reserved. Keywords: El Nin˜o-Southern Oscillation phenomenon; amplitude modulation; nonlinear paradigm

1. Introduction The El Nin˜o-Southern Oscillation phenomenon (ENSO) is the strongest interannual climate mode (Neelin et al., 1998). It can be characterized by an interannual cooling (La Nin˜a) and warming (El Nin˜o) of the eastern equatorial Pacific. Though it originates in the tropical Pacific, it has an impact on weather and climate globally. ENSO has to be regarded as an inherently coupled atmosphere – ocean mode. Eastern equatorial Pacific temperature anomalies are accompanied by tropical wind anomalies, which in turn reinforce the oceanic temperatures. Besides this positive

E-mail address: [email protected] (A. Timmermann).

feedback (Bjerknes, 1969), a delayed negative feedback (see discussion in Neelin et al., 1998) exists which can be provided by wave propagation, advection and other mechanisms in the tropical ocean. It has been observed (Gu and Philander, 1995; Torrence and Webster, 1998) that the amplitude of ENSO undergoes changes on time scales of 10 –20 years. The physical mechanism, which generates these changes, is still unknown. It has been speculated (Gu and Philander, 1997) that temperature anomalies which are subducted in the extratropics have an impact on decadal and probably also on interannual variability in the tropics after entering the equatorial thermocline. However, the subduction hypothesis has also been criticized strongly (Schneider et al., 1999). Several authors (Schneider et al., 1999; Pierce et al., 1999) have argued that subducted temperature anoma-

0921-8181/03/$ - see front matter D 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0921-8181(02)00194-7

136

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

lies attain only values of about 0.1 K, once they enter the equatorial thermocline and that a more likely cause for decadal-scale climate variability in the equatorial Pacific are remote effects of the Pacific Decadal Oscillation (Pierce et al., 1999) or extratropical spiciness anomalies (Schneider, 2000). Thus, it appears as if the interaction between decadal-scale extratropical and decadal-scale tropical variability has been studied quite extensively (only some very few articles are mentioned here). However, the interaction between decadal-scale climate anomalies and interannual ENSO-like variability has been studied only very briefly (Knutson et al., 1997; Timmermann et al., 1999a) using multicentury-long coarse resolution coupled global circulation model (CGCM) simulations. These early simulations suffer from an underestimation of ENSO-like variability. In their studies, Knutson et al. (1997) and Timmermann et al. (1999a) discuss interdecadal changes in the variance of ENSO, which have been simulated with typical time scales of 60 and 35 years, respectively.1 In a most recent study by Timmermann et al. (1999b), a present-day and a greenhouse warming simulation performed with a state-of-the art CGCM (ECHAM4/OPYC) was analyzed in order to understand how ENSO might change in response to greenhouse warming. The authors found that the ENSO amplitude increases with increasing atmospheric CO2 concentrations. While comparing the greenhouse warming simulation with a present-day control run, decadal amplitude modulations with a similar amplitude and time scale were found as in the observations (Fig. 3 in Timmermann et al., 1999b). Our study here is a follow-up inspired by understanding the nature of decadal ENSO amplitude modulations simulated by the ECHAM4/OPYC3 simulation. Three hypotheses will be explored in this paper:


Decadal ENSO amplitude modulations (DEAMs) are generated by the linear superposition of two ‘‘mistuned ENSO’’ oscillators (Section 3).

1 These time scales coincide also with the major time scales of thermohaline circulation (THC) variations in these CGCMs as discussed by Timmermann et al. (1999a). A possible connection between ENSO and the dynamics of the THC has also been studied by Oberhuber (private communication, 1997) published in Latif et al. (2000) using the CGCM ECHAM4/OPYC3.




DEAMs are generated by the effect of a decadal background cycle on ENSO (Section 4).
Nonlinear interactions in the ENSO system give rise both to DEAMs and a decadal tropical mode (Section 5). This paper concludes with a summary and discussion of our results (Section 6).

2. Data The CGCM simulation data to be analyzed come from a 240-year long present-day simulation conducted with the coupled general circulation model ECHAM4/OPYC3. The atmosphere model ECHAM4 is coupled to the isopycnal ocean model OPYC3 by the exchange of freshwater-, heat- and momentum fluxes and by adopting a soft annual mean flux correction (adjustment) method.2 In two recent investigations (Roeckner et al., 1996; Bacher et al., 1998), the ENSO performance of the control integration is validated. It is shown in these studies that the sea surface temperature anomalies (SSTA) related to the El Nin˜o-Southern Oscillation phenomenon have a realistic amplitude and exhibit the typical irregularity of the ENSO signal. The corresponding correlation (teleconnection) patterns of the winds, pressure and precipitation are very similar to those found in observations. However, the simulated ENSO cycle is too short (period of about 2 years) as compared to observations, which is mainly due to an overestimated atmosphere– ocean sensitivity in the tropics (Gro¨tzner, 1998, private communication). It should be noted here that the decadal ENSO amplitude modulations to be discussed below are not an artefact of our model, nor of the flux corrections because the CGCM ECHO-G simulates (Rodgers, private communication) similar 2 The ECHAM4 and the OPYC3 were coupled and integrated for 100 years, while restoring the SST and SSS toward climatological values. Over the 100-year-coupled simulation, the impact of the associated Newtonian feedbacks on the SST and SSS was gradually decreased to zero, while these feedback terms were integrated at full weight in order to derive annual-mean flux adjustments for heat and freshwater. After 100 years of coupled integration, the annual-mean flux adjustments were frozen, and then were applied throughout the subsequent 240-year control run of the coupled model.

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

amplitude modulations in terms of amplitude, period and regularity. The ECHO-G is based on a different ocean model (HOPE) and does not employ flux corrections of any kind.

3. Hypothesis 1: linear ENSO beat ENSO dynamics is very often characterized by the spatially averaged temperature anomalies in the eastern equatorial Pacific. A very popular index region is the Nin˜o 3 region extending from 5jS – 5jN, 150– 90jW. Fig. 1 (upper panel) shows the Nin˜o 3 temperature anomaly time series as simulated by a multicentury present-day climate integration performed with the ECHAM4/OPYC3 model. Apart from the

137

dominant biannual ENSO cycle, an amplitude modulation of ENSO becomes apparent: eras of high and low ENSO variance alternate quite regularly on a decadal time scale. This is further investigated by computing the interannual (1.5 – 7 years) wavelet variance of the Nin˜o 3 SSTA. The result is displayed in Fig. 1. This procedure allows for an extraction of the envelope curve of an amplitude-modulated signal and is discussed in great detail by Gu and Philander (1995) and Torrence and Webster (1998). The envelope curve of the simulated Nin˜o3 SSTA (Fig. 1) exhibits a rather regular oscillation with a time scale of about 16– 17 years. This regular behaviour indicates already that a robust physical mechanism must operate. It is quite unlikely that decadal-scale oceanic processes in the extratropical Pacific such as

Fig. 1. Upper panel: Nin˜o 3 SSTA simulated by the ECHAM4/OPYC3 model (red), interannual wavelet variance (thick black line). Lower panel: same as upper panel but for the observed Nin˜o 3 SSTA.

138

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

subduction (Gu and Philander, 1995) or Rossby wave propagation (Liu, 1999; Latif and Barnett, 1994; Timmermann et al., 1999a) are responsible for such a regular modulation of ENSO variance because they are strongly dependent on extratropical atmospheric conditions and hence, on stochastic excitation. Extratropical SSTA spectra (not shown) as simulated by the ECHAM4/OPYC3 model are basically consistent with the red noise null hypothesis (Hasselmann, 1976) and do not show any preferred decadal spectral peak. As shown in the lower panel of (Fig. 1), the observational Nin˜o 3 SSTA time series exhibits a similar ENSO amplitude modulation as the simulation. The increase of ENSO variance toward the end of the time series might be contaminated by a spatial

undersampling of early SST observations around 1860 – 1950 and hence a possible underestimation of time variance. A possible simple mechanism for the generation of DEAMs is the linear superposition of two oscillatory ENSO modes with slightly mistuned frequencies. In order to test this hypothesis, a Fourier analysis of the Nin˜o 3 SSTA time series (Fig. 1, upper panel) was performed as well as a Singular Spectrum Analysis (Vautard et al., 1992). The results are shown in Fig. 2. The power spectrum (Fig. 2c) is obtained using a Bartlett window with a length of 36 years. It can be compared to a red noise null hypothesis spectrum shown with its 95% confidence limit in Fig. 2c. One observes enhanced variability on time scales of 2– 3

Fig. 2. (a) Nin˜o 3 SSTA and reconstruction based on two SSA modes, (b) running standard deviation of the time series shown in (a), (c) power spectrum of Nin˜o 3 SSTA and the corresponding red noise spectrum, (d) power spectrum of the principal components of the two leading SSA modes of Nin˜o 3 SSTA.

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

years and a modal splitting corresponding to periods of f1 1 = 1.8 and f2 1 = 2.2 years. These peaks are well separated in terms of band width criteria. However, the simulated emphasis of the biannual variability band is unrealistic. The SSA analysis extract also two significant modes with the same biannual time scales as shown by the spectrum of the associated principal components (Fig. 2d), conforming the modal splitting of close-by frequencies. The Nin˜o 3 SSTA time series based on these two SSA modes is reconstructed. The result is depicted in Fig. 2a. As can be easily seen, the match with the original time series is poor. In particular, the decadal beat in the original time series is not reproduced correctly. This is illustrated in Fig. 2b showing the running standard deviation (which is very similar to the wavelet variance in Fig. 1) of the original time series and its SSA mode reconstruction, retaining only these two biannual SSA modes. The linearly reconstructed beat has a time scale, which is too short (12 instead of 16 – 17 years). The envelope period of 12 years is generated by the term ( f1  f2). Hence, the decadal beat cannot be explained simply by a linear superposition of two slightly mistuned ENSO oscillators. Another interpretation of Fig. 1 (lower panel) is that ENSO variability becomes large when a 16-year mode is in a particular phase. This idea resembles very much the phase locking feature of ENSO to the annual cycle, which is a nonlinear phenomenon.3 Standard spectral analysis does not capture nonlinear phasedependencies among oscillators. The signals x1 ðtÞ ¼ cosðx1 t þ /1 Þ þ cosðx2 t þ /2 Þ þ cos½ðx1 þ x2 Þt þ /3 ;

ð1Þ

x2 ðtÞ ¼ cosðx1 t þ /1 Þ þ cosðx2 t þ /2 Þ þ cos½ðx1 þ x2 Þt þ /1 þ /2 ;

ð2Þ

have the same power spectrum with peaks at x1, x2, x1 + x2. However, if /i(t) are randomly sampled 3

It is important to note the difference between nonlinear phase locking and linear instability of a seasonally varying background state. Whether the real ENSO system is phase-locked or whether the preferred occurrence of El Nin˜o and La Nin˜a events during boreal winter is based on Floquet dynamics has not been resolved yet.

139

phases, the two processes, Eqs. (1) and (2), differ fundamentally. Eq. (2) characterizes a coupled oscillatory system, where an oscillation with frequency x3 = x1 + x2 is generated due to a nonlinear interaction of oscillators with frequencies x1 and x2. Summarizing higher order harmonics at frequencies 2x1, 2x2 and x1  x2 in the term R(t), x2(t) can be generated by interactions of the type x2 ðtÞ ¼ xðtÞ þ ax2 ðtÞ  RðtÞ;

ð3Þ

with x(t) = cos(x1t + /1) + cos(x2t + /2). Because of this relation, the type of phase coupling in Eq. (2) of x2(t) is called quadratic phase coupling. A very convenient way to characterize nonlinear phase interactions among oscillators of different frequencies is bispectral analysis, which is described in the Appendix A. For our purposes, a peak in the bispectrum indicates where phases of individual oscillatory components are statistically linked, and hence, indicates quadratic phase coupling. Eqs. (6) – (10) from Appendix A are applied to the simulated Nin˜o 3 SSTA time series. The result is shown in Fig. 3 One observes strong phase coupling for the on-diagonal frequency pairs (0.0225, 0.0225), (0.028, 0.028), (0.038, 0.038) [month 1], and off-diagonal interactions, e.g. at (0.016, 0.023) [month 1]. Furthermore, one observes enhanced bispectral power at the frequency pairs (0.0045, 0.032), (0.0045, 0.045) [month 1]. This is a manifestation of nonlinear interactions among decadal (18 years) and interannual (2 years) variability. This conforms our hypothesis that ENSO is phase coupled to a decadal oscillation and vice versa. However, it should be noted that there exists also the possibility that phase coupling among oscillators does not manifest itself in amplitude modulations. Here, bispectral analysis was used in order to show that a nontrivial interaction among ENSO and decadal variability occurs. Physical and mathematical origin of this interaction will be discussed in subsequent sections. It should be noted here that statistical significance tests for bispectra can be performed either by using Hinich’s v2 test criterion (Hinnich, 1982) or by running a set of Monte Carlo experiments based on phase-shuffled surrogates (Theiler et al., 1992). It is important to note that the requirements to data statistics are much higher for bispectra than for first order spectra, a problem which was already pointed out by

140

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

Fig. 3. Amplitude of bispectrum computed from the Nin˜o 3 SSTA time series in Fig. 1 (upper panel). x-Axis and y-axis represent frequency f1 and f2. Frequency unit is month 1.

Barnett (1991). Five thousand surrogate data series were simulated with the same spectrum as the Nin˜o3 SSTA. Bispectra are computed for this Monte Carlo ensemble and Fig. 3 is compared with the mean and variance of the ensemble. The features discussed here exceed the 80% confidence limit, obtained from the phase-shuffled surrogates.

4. Hypothesis 2: decadal ‘‘background’’ mode generates DEAMs Before analyzing the physical mechanism which generates DEAMs, the simulated mean state in the tropical Pacific has to be assessed. It might have an impact on the variability and the decadal phase locking of ENSO.

Fig. 4 displays the mean simulated surface temperatures and currents in the tropical Pacific. The position and strength of the North Equatorial Countercurrent agrees well with the observations (Lagerloef et al., 1999; Ji et al., 1994), whereas the South Equatorial Current is somewhat too strong. The meridional velocity components are mainly due to Ekman drift, whereas the zonal velocities have both a significant Ekman and a geostrophic contribution. The coupled model is flux corrected in the annual mean. Due to this correction, the mean SST structure (Fig. 4) is very realistic. The transition between warm pool and cold tongue, characterized by large temperature gradients, is located between 140jW and 120jW. Later on, it will be shown that this transition region is important for the generation of tropical decadal anomalies.

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

141

Fig. 4. Simulated mean tropical surface currents (left, unit m/s) and temperatures (right, unit K).

The interaction between long-term (decadal) climate variability and ENSO is investigated. The aim is to identify a decadal oscillation, which is associated with DEAMs as discussed in the previous section. The discussion of the causality issue will be postponed to a later section. An EOF analysis of the tropical Pacific subsurface temperature anomalies in a 100-m depth is performed (note the isopycnal ocean data have been interpolated on z-coordinates). Similar results can be obtained for the surface and somewhat deeper levels. The depth of 100 m yields an optimal signal to noise

ratio. The leading EOF (not shown) is associated with the typical ENSO mode (similar to Fig. 1 in Timmermann, 1999). The second mode is depicted in Fig. 5. It characterizes an east –west displacement of ENSO’s center of action. Associated with the second EOF of subsurface temperature anomalies is a deepening of the thermocline in the far eastern equatorial Pacific as well as a reduction of the thermocline slope. The corresponding 2nd principal component (normalized and filtered with a 10-year running mean) is shown in Fig. 6. Comparison with the running stand-

Fig. 5. Simulated 2nd EOF pattern of subsurface (100 m) temperature anomalies. Unit is K.

142

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

Fig. 6. Normalized and filtered (10-year running mean) principal component of the 2nd EOF mode of subsurface temperature anomalies (solid) and normalized time series of the running standard deviation of Nin˜o 3 SSTA anomalies (window length 10 years).

ard deviation of Nin˜o 3 SST anomalies yields a very high correlation of 0.96.4 It can furthermore be shown that temperature anomalies averaged over the far eastern equatorial Pacific reveal a close connection between their low pass-filtered decadal first, second and third order running moments (not shown). It turns out that ENSO variability is strong during phases when El Nin˜o events are shifted more to the east, and La Nin˜a events being shifted more to the west. The amplitude of ENSO is small in the reverse case. This result is consistent with the results presented in An and Wang (2000). A spectrum analysis (not shown) of the second principal component of tropical subsurface temperature anomalies reveals a peak at 16– 18 years, which is significant even above the 99% confidence limit. Similar peaks could not be identified in the extratropics. It should be noted here that the decadal dipole 4

This high correlation can have several reasons. (a) A decadal mode causes decadal variations in ENSO variance. (b) Decadal ENSO variance and background state changes are generated by the same mechanism. (c) ENSO is characterized by a skewed probability distribution. Due to the skewness, decadal variations in ENSO variance will automatically manifest themselves in decadal changes of the mean state.

mode has a very close similarity also to the decadal mode found by Kleeman et al. (1999). A region which is dynamically very important for the generation of the ‘‘Kleeman et al. mode’’ is the subtropical east Pacific close to Baja California. In contrast to their findings, no significant decadal mode could be identified in our simulation using only data north of 20jN and different statistical techniques (POPs, EOFs and area averaged indices) thereby favoring a tropical mode hypothesis. Nevertheless, it is interesting to study the dynamics of the extratropical Pacific and how it could be linked to decadal changes in ENSO variability. Therefore, I perform a linear lag regression analysis between the running standard deviation ENSO time series shown in Fig. 1 and oceanic fields (time series of oceanic quantities at each grid point). In Fig. 7, the regression between ENSO variability and the anomalous depth of the density level r = 1026 kg m 3 is shown. This density level captures thermocline dynamics. The subsurface dynamics in the North Pacific (shown in Fig. 7) can be characterized by two main features: a westward propagating anomaly at latitudes of 8j0 –15jN and a southwestward propagating anomaly which is generated at about 160jW and 35jN.

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

143

Fig. 7. Lag regression patterns computed between the running standard deviation of Nin˜o 3 SSTA and anomalies in the depth of the 1026 isopycnal.

These features are reminiscent of the 1st and 2nd mode described in the wave theory of the ventilated thermocline (Liu and Zhang, 1999; Liu, 1999). According to this theory, the 1st mode is most strongly pronounced in the lower thermocline. It is mainly generated by wind forcing and has a characteristic propagation time of 5 years. The second mode is advected by the mean thermocline circulation. It follows the southwestward subduction path in the North Pacific. Typical time scales are on the order of 10 years. This mode is easily excited by surface buoyancy forcing. This scenario is consistent with Fig. 7. A comparison of Fig. 7 with the observed subsurface temperature evolution as described in Zhang et al., 1998 (Fig. 1d, therein) reveals further similarities: during phases of high ENSO variability, the eastern equatorial Pacific is warm with a warm tongue extending from Baja California to the dateline. The western equatorial Pacific as well as the

central North Pacific are relatively cool. Furthermore, the dynamical evolution on decadal time scales is similar. Fig. 1 in Zhang et al. (1998) can be directly compared with the last four lag regression patterns shown in Fig. 7. A crucial question to be answered here is whether these decadal subsurface temperature anomalies are causing DEAMs or are caused by DEAMs. In particular in the light of recent studies by Pierce et al. (1999), the possibility should be taken into account that decadal anomalies in the tropics are mainly generated by tropical wind anomalies (originating from extratropical large-scale atmospheric patterns) rather than by the oceanic subduction bridge. Our study reveals no conclusive evidence that the simulated decadal anomalies in the tropics are generated by extratropical subduction processes. In particular, the absence of decadal spectral peaks in the extratropics weakens the possi-

144

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

bility of extratropical influence on DEAMs. However, the reverse might very well be the case. In particular, the strong surface and subsurface temperature anomalies associated with Fig. 5, which are on the order of 2– 3 K, have global teleconnections and encompass also the extratropics. Changes in the strength of the Aleutian low, which accompany decadal changes of ENSO variability (not shown), modify oceanic forcing by extratropical wind stress curl and heat flux anomalies. Thus, oceanic subduction and Rossby wave propagation can be excited by decadal temperature changes in the tropical Pacific. In the South Pacific, one observes a westward propagation of temperature anomalies into the tropics. A negative anomaly is generated at lag  8 years at about 25jS. This negative temperature (or thermocline depth) anomaly propagates to the west coast of Australia. It appears as if this anomaly ‘‘flushes’’ the equatorial Pacific at lag

 2 and 0. However, an analysis of the squared coherence computed between thermocline depth anomalies off the Australian coast and the running standard deviation of Nin˜o 3 SST anomalies results in values of about 0.3 – 0.6, thus explaining only a moderate percentage of decadal variance in ENSO variability. In order to study further how ENSO variability changes in conjunction with the 2nd EOF mode shown in Fig. 5, more oceanic fields are analyzed. Fig. 8 shows the lag regression patterns between ENSO variance and the equatorial fields of temperature anomalies. One observes that during phases of high ENSO variance (lag 0), the thermocline slope is reduced. Furthermore, one observes a deepening of the thermocline in the eastern equatorial Pacific associated also with anomalous downwelling (not shown). These changes are expected to go along with less ENSO variability (see Zebiak and Cane, 1986; Meehl

Fig. 8. Lag 0 regression pattern computed between the running standard deviation of Nin˜o 3 SSTA and anomalies of subsurface temperatures along the equator. y-Axes represent depth [m], x-axes represent longitude.

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

et al., 2001), in contrast to our findings. Typical temperature anomalies at the surface attain values of about 2 – 3 K. However, a deepening and slope reduction of the equatorial thermocline is expected to go along with reduced ENSO activity (Zebiak and Cane, 1986) in contrast to our findings. This counterintuitive behaviour has also been identified in the nonfluxcorrected Coupled General Circulation simulation performed with ECHO-G (ECHAM4/global HOPE) model (Rodgers, private communication). The time scale of the simulated ENSO amplitude modulation as well as its regularity are very similar to the ECHAM4/ OPYC3 run. Our discussion revealed that an additional mechanism has to operate which overcompensates the effect of thermocline slope reduction and thermocline deepening. One of the most dominant signatures of the decadal tropical mode shown in Fig. 5 is the temperature dipole along the equator with temperature anomalies, which attain temperatures of a mature El Nin˜o or La

145

Nin˜a event. An important question to be answered is which processes determine these decadal temperature anomalies. The temperature tendency equation for the mixed layer is divided into its components: heat flux, mixing, advection (horizontal and vertical), convection and diffusion. For each of the components, a regression is computed onto the time series of the running standard deviation of ENSO. The result is depicted in Fig. 9. It becomes evident that the strong positive temperature anomaly in the eastern equatorial Pacific seen in Fig. 5 is generated by the interplay of two major components: advection and mixing. Heat flux and diffusion generate negative SST tendencies, which are overcompensated by the two other components. This result stresses the role of advection (horizontal and vertical). Summarizing our results, our lag regression analysis has revealed the following: DEAMs are accompanied by a decadal oscillation which is characterized by a tropical temperature dipole extending from the surface to depths of about 150 m. Large ENSO

Fig. 9. Regression at lag 0 computed between the running standard deviation of Nin˜o 3 SSTA and equatorial anomalies of heat flux, mixing, advection, convection and diffusion in the mixed layer.

146

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

variability is correlated with a reduction of the thermocline slope and a deepening in the eastern equatorial Pacific. The eastern equatorial temperature anomaly is generated by advection and mixing. Furthermore, it has been shown that the ocean dynamics in the North Pacific and South Pacific fluctuate in conjunction with the tropical dipole. In the North Pacific, propagation of thermocline depth anomalies is reminiscent of both, subduction of temperature anomalies and propagation of Rossby waves. However, there is hardly any evidence that extratropical signals cause DEAMs as speculated by Gu and Philander (1997), Weaver (1999) and Schneider (2000) and others. Correlations among low passfiltered SSTAs in the tropics and extratropics attain (significant) values of about 0.3– 0.45—too little to build a theory on. It should further be noted that our statistical analysis based on the mean – variance relation (see Hypothesis 2) accounts only for a very small part of possible nonlinear relations, namely for the linear regression hxVhxV2ii. Furthermore, the possibility has to be taken into account that the decadal ENSO modulation and the decadal ‘‘background’’ cycle are two sides of one coin rather than generating one another.5 This idea will be further discussed in the next section.

5. Hypothesis 3: DEAMs are generated by nonlinear tropical dynamics Our goal is to derive an empirical low-order ENSO model from the CGCM data under consideration. It is required that the reduced model captures the most important dynamical features of ENSO as simulated by the CGCM. This task will be achieved by applying multiple nonparametric regression techniques to EOFtruncated temperature and sea level anomalies obtained from the multicentury-long CGCM simulation. This methodology has been applied successfully in Timmermann (2001) and Timmermann et al. (2001). A brief summary of this method can be found in Appendix B.

5

This does not mean that the dependence of ENSO on its mean state is denied. The only thing which is stated here is that this dependence is not the complete story.

An EOF analysis of the simulated sea surface temperature anomalies and sea level anomalies (which correlate very well with thermocline depth anomalies) is computed based upon the multicentury-long ECHAM4/OPYC3 simulation. The two leading temperature and sea level EOFs for the equatorial Pacific are depicted in Fig. 10. The first SSTA mode corresponds to the mature ENSO phase (El Nin˜o or La Nin˜a), whereas the second EOF mode is characterized by a dipole-like structure. It is associated with a zonal displacement of the ENSO-related SSTA. When the second SSTA EOF mode has a positive principal component, the El Nin˜o center will be shifted eastward, whereas a La Nin˜a situation will be focussed on the equator more westward to its normal position. The opposite holds for negative values of the second principal component. The first EOF mode of thermocline variability describes an east – west thermocline see-saw mode, consisting of off-equatorial signals in the warm pool region and equatorial anomalies in the central and eastern equatorial Pacific. The second sea level EOF can be interpreted in terms of a Kelvin wave travelling along the equator and an off equatorial Rossby wave. The corresponding principal components (T1, T2, H1, H2) for 40 simulation years are displayed in Fig. 11. The leading temperature and sea level EOFs are characterized by a noisy oscillation with a period of about 2– 3 years. The second principal component of temperature anomalies is much noisier than the first one. It appears also that when T2 attains large values, the associated interannual variability of T1, H1, H2 is also high. In order to apply the ACE algorithm I set (x1, x2, x3, x4)=(T1, T2, H1, H2). By means of the optimal transformation technique described in Appendix B, a set of dynamical equations is derived from the four principal components of the ECHAM4/OPYC3 model shown in Fig. 11. Lookuptables for the general functions Uji are obtained which are fitted using polynomials of 4th order. The associated ODEs are integrated forward in time using a Runge– Kutta scheme. The resulting ‘‘reconstructed’’ principal components (Fig. 12) reveal a very interesting feature: The amplitude of ENSO, as represented by T1, is amplitude modulated on a time scale of about 16 years. The second simulated principal component T2 follows this amplitude modulation in such a way that high values of T2 correspond to high variance of T1. This is in close agreement with the original CGCM

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

147

Fig. 10. Dominant EOF patterns of sea surface temperature anomalies (upper panel) and sea level anomalies (lower panel). Data are obtained from the first 40 years of a greenhouse warming simulation performed with the CGCM ECHAM4/OPYC. Patterns are dimensionless.

data as discussed in the previous sections (see Fig. 6), although the time series simulated with our empirical model look much more regular. Hence, the empirical dynamical model is able to capture the relation between ENSO variance and the east – west temperature dipole shown in Fig. 5. Furthermore, the time scale of the amplitude modulation is reproduced, despite of the fact that the empirical model is derived only from 40 years of data. This astonishing capability of the model obtained from the ACE algorithm is discussed more extensively in Timmermann et al. (2001). Motivated by this success, sensitivity experiments are performed with the empirically derived ENSO model, rather than with the computationally expensive CGCM. The general caveats of this approach are also discussed in Timmermann et al. (2001).

If our working Hypothesis 2 were true and decadal changes in ENSO variability were triggered by a decadal background cycle6, one should be able to simulate decadal amplitude modulations of T1 by prescribing the T2 component using a decadal sinusoidal forcing with an amplitude of 8– 12. I performed a series of such forced experiments. In none of these experiments, one observes decadal ENSO amplitude modulations in T1. Assuming that our four-dimensional empirical ENSO model captures the main features of the ENSO dynamics simulated by the CGCM, one can conclude that Hypotheses 2 is inappropriate to explain decadal ENSO amplitude 6 Which is associated with the east – west dipole pattern in Fig. 5.

148

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

Fig. 11. Leading principal components for temperature anomalies and sea level anomalies simulated by the ECHAM4/OPYC3 CGCM. The thick line in the upper right panel corresponds to a 50-month running mean of T2.

modulations. Decadal changes in the background state do not trigger decadal ENSO amplitude vacillations. This means that the stated ‘‘counter-intuitiveness’’ of our linear analysis might be a misleading interpretation of the underlying nonlinear dynamics. In our case here, ENSO amplitude modulations can only be simulated when the full nonlinear four-dimensional dynamical equations are used. This illustrates the crucial difference between forced decadal ENSO amplitude modulations and coupled ones. In order to study the sensitivity of the empirical ENSO model to parameter changes, the interaction between the first and second principal components for temperature is modified. This is achieved by introducing a bifurcation parameter l into x˙ 1 = . . . + lU12(x2) + . . . Two additional experiments are performed for l = 1.04, 1.2. The simulated principal components are depicted in Figs. 13 and 14, respectively. As compared to Fig. 12, one observes that an increase in the interaction coefficient l goes along with an increase in the length of the amplitude modulations, because the large decadal changes in T2 are

interrupted by interannual events of small amplitude. This increases the return time of large amplitude events. The overall features, however, agree with Fig. 12. In case of further enhancement of l to values of about 1.2 (and further), one sees that the amplitudemodulated character of ENSO disappears and the T1 component is characterized by an interannual oscillation with a period of 2– 3 years, in accordance with the original ECHAM4/OPYC3 data. From these experiments, it becomes evident that the dynamical regime of ENSO depends strongly on the nonlinear interaction coefficient l. I will study this feature more in detail by computing phase-space sections of the attractors for l = 0.98, 1.04, 1.08, 1.2 through the T1  H1 plane. The results are shown in Fig. 15. 5.1. Homoclinic ENSO dynamics In Fig. 15a, one observes the following phenomenon: starting with large negative T1, H1 values, one observes that the dynamics spirals into the vicinity of a saddle node which had been identified from linear

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

149

˜k Fig. 12. P Simulated leading principal components for temperature and sea level anomalies, using the optimal transformations Ui of the model x˙ k ¼ Ai¼1 U˜ ki ðxi Þ , with A = 4, k = 1,. . .4. The numerical integration of the corresponding equations of motion is performed by 4th order polynomial fitting the optimal transformations. The interaction coefficient is l = 0.98.

Fig. 13. Same as in Fig. 12. The interaction coefficient is l = 1.04.

150

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

Fig. 14. Same as in Fig. 12. The interaction coefficient is l = 1.2.

Fig. 15. Phase space plots in T1  H1 plane for different values of the nonlinear interaction coefficient l = 0.98 (a), 1.04 (b), 1.08 (c), 1.2 (d).

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

151

Fig. 16. (a) Schematic illustration of a homoclinic bifurcation. The bifurcation parameter l is increased from left to right.

stability analysis (Timmermann, 2001). Not hitting this stationary point exactly activates the repelling properties of the real positive eigenvalue and the trajectory is expelled from the saddle node neighborhood, attaining large negative temperature as well as thermocline depth values. Nonlinearities slow down the linear instability and the system spirals back into the origin due to the increasing importance of the complex conjugate eigenvalue pair with negative real part. The physical importance of the saddle node is that it represents an oscillatory stable climate mean state. The oscillatory stability leads to damped ENSO oscillations. In another direction, one observes a nonoscillatory growing instability, which is initially characterized by almost neutral temperature conditions and a deepening of the equatorial thermocline. The dynamical evolution shown in Fig. 15a exhibits very close similarities to a homoclinic orbit.7 It should be noted here that for nonlinear interaction coefficients smaller than 0.95, the dynamics gets unstable and the simulated principal components grow to infinity. For values l greater than 0.95, the simulated principal components stabilize and oscillatory behaviour of different complexity emerges. This is illustrated in Fig. 15b,c,d. For l = 1.04, the phase-space trajectory looks more complicated than for l = 0.98. However, still the dynamics is periodic. When l attains values of about 1.08, the dynamics becomes chaotic. The overall dynamical evolution, however, remains very similar: damped ENSO oscillations, which grow in amplitude, once a critical vicinity of the saddle point has been 7

A homoclinic orbit occurs when a stable system manifold of a saddle focus (or point) merges with one of the ends of its unstable manifold to form a closed loop.

reached. For larger values of l, the dynamics simplifies and the typical amplitude modulation character vanishes. These features are reminiscent of homoclinic bifurcations. Such a homoclinic bifurcation is schematized in Fig. 16 for a three-dimensional case. Our empirical model is four-dimensional, thereby allowing for even more complexity. Nevertheless, the basic structure of a homoclinic bifurcation in three dimensions can be recovered qualitatively from Fig. 15. An important property of homoclinic bifurcations is that the system can exhibit arbitrarily long periods in the parameter vicinity of homoclinic orbits. Very often, however, the associated long-term oscillations are unstable and overlap such as to generate chaos. The main message to be taken from this section is that homoclinic bifurcations provide a concept that might explain the generation of amplitude modulations with long periods.8 In higher dimensions, the existence of a homoclinic orbit in parameter space is very likely to generate chaos in the parameter vicinity of the bifurcation point (Glenndinning and Sparrow, 1984). This scenario is also called Shil’nikov (1965, 1970) scenario.

6. Summary and discussion It is widely accepted that the amplitude of ENSO varies from decade to decade. It has not been understood why. The decadal variations in ENSO amplitude 8 The average return time of large events does not depend solely on the linear eigenvalues but on the strength of the nonlinearity, because it determines how often the system spirals around the saddle node until the amplitude cycle repeats itself.

152

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

might just be an expression of ENSO being partly stochastically excited or of dynamical processes interacting with ENSO. In particular in the light of recent observational studies (Torrence and Webster, 1998; Gu and Philander, 1995), which indicate that the amplitude of ENSO is modulated on a rather distinct time scale of 15 – 20 years, the latter possibility is quite likely and has to be investigated more in detail. Our analysis was based on a multicentury presentday climate simulation performed with ECHAM4/ OPYC3 model. This ENSO-resolving model simulates an amplitude-modulated ENSO signal with a typical envelope time scale of 15 – 20 years, similar to observations. This finding inspired us, despite of some model caveats, to study the simulated amplitude modulation statistically. Our analysis explored firstly on the possibility that the DEAMs arise from two linearly connected ENSO oscillators with slightly different internal frequencies. Evidence from spectral analysis suggested, however, that this hypothesis cannot explain the simulated DEAMs easily. Moreover, a bispectral analysis of Nin˜o 3 SST anomalies revealed that ENSO is interacting with decadal ‘‘background’’ oscillation in a nonlinear manner. Our study continues with the next working hypothesis: A decadal background cycle in the Pacific influences ENSO variability. The crucial idea, which was pursued, is that slow changes in the mean state of the tropical Pacific might have an impact on the characteristics of ENSO variability (structure, frequency, amplitude, growth rate). Our statistical analysis revealed, in fact, the existence of a decadal mode in the tropical Pacific, associated with an east – west displacement of ENSO’s center of action. This decadal mode is highly correlated with the simulated changes in ENSO variance. Furthermore, due to the strong SST anomalies generated by this mode, it encompasses also the extratropics in both Hemispheres and the North Atlantic. Evidence for extratropical dynamics causing this mode is weak. The most striking finding of our statistical analysis was that this decadal mode accompanies changes in ENSO variability in a counterintuitive manner: phases of high ENSO variability go along with reduced zonal and vertical temperature gradients on the equator. This finding turned out to be inconsistent with the second working hypothesis. The third working hypothesis which was explored is the following: decadal changes of ENSO variance as well as the existence of a decadal tropical mode are a

manifestation of one dynamical structure, which is based on coupled nonlinear dynamics and do not interact with each other in a simple linear cause-andeffect relation. Using an empirical low-order ENSO model, which was derived from ECHAM4/OPYC3 data, similar modulations of ENSO amplitude are obtained which are accompanied by decadal variations in the zonal gradient of the thermocline anomalies as simulated by the CGCM, although the time series of the low-order model are naturally less noisy. Our analysis suggested that the long time scale cannot be explained simply in terms of linear dynamics but that nonlinear interactions have to be taken into account. It turned out that due to the interplay of the first and second temperature anomaly, EOF long time scale dynamics can be generated which can be described in terms of the Shil’nikov scenario. Furthermore, the low-order model simulates a similar counterintuitive behaviour as seen in the CGCM. Due to the role of nonlinearities, which can support also chaotic oscillations, the construction of a simple ‘‘intuitive’’ physical concept of DEAMs is highly nontrivial. Similarly, though studied for over a decade now, the physical mechanism generating ENSO’s phase-locking to the annual is still not fully disentangled. It is even not clarified as to whether ENSO’s preferred occurrence during boreal winter is due to nonlinear phase-locking or due to cyclic instabilities (Floquet dynamics). I illustrated that a mathematical scenario, known as the Shil’nikov scenario can generate (arbitrarily) long time scales, due to internal nonlinear dynamics. Translated to the climate context, this means that long-term climate variability can be generated far beyond the time scales suggested by linear analysis.9 Furthermore, the paradigm of homoclinic ENSO orbits predicts the emergence of amplitude modulations and irregular oscillations. This new paradigm for ENSO had been verified in recent papers (Timmermann and Jin, 2002a,b) using the nonlinear ENSO recharge oscillator (Jin, 1997) and the Zebiak and Cane (1986) model. This analysis revealed also that preferred decadal (amplitude modulation) time scales can be generated in the tropics by nonlinear interactions. This finding is consistent with the modeling results of Kirtman and Schopf (1998).

9

Very often used in ENSO theory.

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

In Fig. 1, it can be seen that phases of high ENSO variability are characterized also by more regular ENSO variability, whereas relatively calm phases seem to be noisier. This finding might suggest also that ENSO predictability and irregularity are modulated by DEAMs. The question arises also as to whether the observed changes of ENSO predictability (Balmaseda et al., 1995) are related to observed DEAMs (discussed, e.g. by Torrence and Webster, 1998). These questions will be investigated more in detail in forthcoming studies. Another question, which is important to be addressed, is whether the envelope curve of ENSO is predictable on longer time scales. A prediction of high and low ENS0 variance regimes would be of high societal value (Timmermann and Jin, 2002a). Such a prediction exercise might profit from the theory of heteroscedastic processes (Horn et al., 1975), which are characterized by time-varying variances and which are used in econometrics. Our future research will be devoted to establishing our nonlinear paradigm for ENSO amplitude modulations using models of intermediate complexity such as the Zebiak and Cane (1987) model and others.

where h. . .i is an ensemble averaging operator. The bispectra of x1 and x2 are: Bx31 ðX1 ; X2 Þ ¼ 0;

ð5Þ

Z Z 1 l l ¼ ds1 ds2 ½cosðx2 s1 þ x1 s2 Þ 4 l l þ cosðx3 s1  x1 s2 Þ þ cosðx1 s1 þ x2 s2 Þ

Bx32 ðX1 ; X2 Þ

þ cosðx3 s1  x2 s2 Þ þ cosðx1 s1  x3 s2 Þ þ cosðx2 s1  x3 s2 ÞeiðX1 s1 þX2 s2 Þ :

ð6Þ

It becomes evident that the bispectrum of x2 is nontrivial and exhibits a peak at x1, x2 which is the frequency pair for which quadratic phase coupling occurs. Eq. (4) will be applied to the Nin˜o 3 SSTA time series shown in Fig. 1 (upper panel). Similar to the power spectrum analysis, a two-dimensional Parzen window is used. The resulting discretized equations read (see Nikias and Petropulu, 1993): B3 ðX1 ; X2 Þ ¼

Acknowledgements

153

L L X X

ˆ 1 ; s2 Þwðs1 D; s2 DÞeiðX1 s1 þX2 s2 Þ ; ð7Þ mðs

s1 ¼L s¼L

The author is very grateful to Dr. Henning Voss for sharing his deep insight into the field of nonlinear data analysis and for providing the ACE algorithm. The author wishes to thank Dr. Fei-Fei Jin for the numerous discussions on Homoclinic ENSO dynamics and J.M. Oberhuber for his support in handling the CGCM data. This work is supported by a DFG fellowship. The author is also grateful to Annegret Schurbohm for her graphical assistance.

with ˆ 1 ; s2 Þ ¼ mðs

s2 K 1X 1 X xi ðkÞxi ðk þ s2 Þxi ðk þ s2 Þ; K i¼1 M k¼s 1

ð8Þ wðs1 D; s2 DÞ ¼ dðs1 DÞdðs2 DÞ;

ð9Þ

The Parzen window is defined as: Appendix A . Bispectral analysis The bispectrum for a given signal x(t) is defined as:

B3 ðX1 ; X2 Þ ¼

Z

l

l

Z

8 1  6u2 þ 6AuA3 ; AuAV0:5 > > > > < dðuÞ ¼ 2ð1  AuAÞ3 ; 0:5 < AuAV1 : > > > > : 0; otherwise

ð10Þ

l

ds1 ds2 hxðtÞx l

 ðt þ s1 Þxðt þ s2 ÞieiðX1 s1 X2 s2 Þ ;

ð4Þ

The indices have the following meaning. For averaging the third order moment, the time series of length N is split into K records of M samples (N = KM). xi(k) is

154

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

the data per segment (i = 1,2,. . .K). s1 = max(0, s1, s2), s2 = min(M  1, M  1  H 1, M  1  s2). L denotes the region of support of the estimated 3rd order moment function. The following relation holds: js1j V L, js2j V L. The band width D is determined by D = 1/L.

Appendix B . Deriving an empirical ENSO model In the context of data-based dynamical system modeling, one very often encounters the problem to determine the relation F in the model: x˙ tj ¼ F j ðxt Þ;

ð11Þ

xt=(xt1,. . .,xti,. . .,xtN)T is the N-dimensional data at time t and x˙ tj represents the time derivative of

where vector the j-th component at time t. In the following, I omit the time index t. In order to find the optimal model function F˜ which solves Eq. (1) in a least-square sense, the expectation value of x˙ j conditioned on x has to determined. In general, the optimal function can be estimated as the conditional expectation value: Z ˜F j ðxÞ ¼ pð˙x j j xÞ˙x j d˙x j ; ð12Þ where p(.j.) denotes the conditional probability density function. Due to the relation p(x˙ jjx) = p(x˙ j, x)/ p(x), this problem becomes equivalent to the estimation of a multidimensional probability density function. Even for small system dimensionality N, this is often practically impossible. In climate dynamics, observational datasets are too short and too noisy in order to apply this strategy successfully. The following letout is proposed: For many problems, it is not necessary to solve the general reconstruction problem. On the basis of some a priori knowledge of the probability density function, the reconstruction problem simplifies considerably. It is assumed that the model Eq. (11) can be approximated by a linear combination of general transformations Uij(xi), i.e.: x˙ j ¼

N X

Uij ðxi Þr j nj ;

ð13Þ

i¼1

where the mismatch between Eqs. (11) and (13) is parameterized in terms of white noise hnj(t)nk(tV)i = djkd(t  tV) of variance r j2.

Eq. (13) defines a nonlinear multiple regression problem. It is solved with the Alternating Conditional Expectation Value (ACE) algorithm of Breiman and Friedman (1985). The ACE algorithm converges to consistent estimates of so-called optimal transformations, which are the least-squares solution of Eq. (13). Here, a modified algorithm is used, where the modification consists of fixing the l.h.s. of Eq. (13), whereas in Breiman and Friedman (1985), it is also varied. For details on the numerical implementation of the ACE algorithm, see Voss and Kurths (1997), Voss (2001) and Timmermann et al. (2001). The modified ACE algorithm works as follows: Globally optimal solutions (in a least-square sense, i.e. hð˙xj  SNi¼1U˜ ji ðxi Þ2 i ¼ min) can be obtained iteratively from the sequence: j ðx i Þ ¼ h˙x j j x i i; U˜ i;0

* U˜ i;kj ðxi Þ

¼

ð14Þ  +   ðx Þx i : p;k * 

N X x˙  ; U˜ j j

ppi

p

ð15Þ

The index j corresponds to the component of the differential equation, k to the iteration step (k>0) and i to the sum over the predictor components. The index k* equals k for p < i and k  1 for p>i, and h. . .j. . .i denotes the conditional expectation value. In expression (15), only scalar quantities are involved, and in contrast to Eq. (12), only one-dimensional conditional probabilities (or, equivalently, two-dimensional joint probabilities) have to be estimated. These can be interpreted in terms of the time transition probabilities or as the dynamical contribution to the componentwise Frobenius – Perron operator of the underlying ODE system. Minimizing hð˙xj  SNi¼1 U˜ ij ðxi ÞÞ2 i is equivalent to the maximization of the correlations: * + N X j i U˜ i ðx Þ x˙ j Wj ð˙xj ; x1 ; . . . ; xN Þ ¼ 0 @h˙x j 2 i

i

*"

N X

#2 +11=2 ; A U˜ ij ðxi Þ

i

ð16Þ where it is assumed that all variables have zero mean. Hence, this technique to solve the nonlinear regres-

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

sion problem is also called maximal correlation approach.

References An, S.-I., Wang, B., 2000. Interdecadal changes in the structure of ENSO mode and its impacts on the ENSO frequency. J. Climate 13, 2044 – 2055. Bacher, A., Oberhuber, J.M., Roeckner, E., 1998. ENSO dynamics and seasonal cycle in the tropical Pacific as simulated by the ECHAM4/OPYC3 coupled general circulation model. Clim. Dyn. 14, 431 – 450. Balmaseda, M.A., Davey, M.K., Anderson, D.L.T., 1995. Decadal and seasonal dependence of ENSO prediction skill. J. Climate 8, 2705 – 2715. Barnett, T.P., 1991. The interaction of multiple time scales in the tropical climate system. J. Climate 4, 269 – 285. Bjerknes, J., 1969. Atmospheric teleconnections from the equatorial Pacific. Mon. Weather Rev. 97, 163 – 172. Breiman, L., Friedman, J.H., 1985. Estimating optimal transformations for multiple regression and correlation. J. Am. Stat. Assoc. 80, 580 – 598. Glenndinning, O., Sparrow, C., 1984. Local and global behaviour near homoclinic orbits. J. Stat. Phys. 35, 645 – 697. Gu, D.-F., Philander, S.G.H., 1995. Secular changes of annual and interannual variability in the tropics during the past century. J. Climate 8, 864 – 876. Gu, D.-F., Philander, S.G.H., 1997. Interdecadal climate fluctuations that depend on exchanges between the tropical and extratropics. Science 275, 805 – 807. Hasselmann, K., 1976. Stochastic climate models: Part 1. Theory. Tellus 28, 473 – 485. Hinnich, M.J., 1982. Testing for Gaussianity and linearity of a stationary time series. J. Time Ser. Anal. 3, 169 – 176. Horn, S.D., Horn, R.A., Duncan, D.B., 1975. Estimating heteroscedastic variances in linear models. J. Am. Stat. Assoc. 70, 380 – 385. Ji, M., Leetmaa, A., Derber, J., 1994. An ocean analysis system for seasonal to interannual climate studies. Mon. Weather Rev. 123, 460 – 481. Jin, F.F., 1997. An equatorial ocean recharge paradigm for ENSO: Part I. Conceptual model. J. Atmos. Sci. 54, 811 – 829. Kirtman, B.P., Schopf, P.S., 1998. Decadal variability in ENSO predictability and prediction. J. Climate 11, 2804 – 2822. Kleeman, R., McCreary Jr., J.P., Klinger, B.A., 1999. A mechanism for generating ENSO decadal variability. Geophys. Res. Lett. 26, 1743 – 1746. Knutson, T.R., Manabe, S., Gu, D., 1997. Simulated ENSO in a global coupled ocean – atmosphere model. Multidecadal amplitude modulation and CO2 sensitivity. J. Climate 10, 138 – 161. Lagerloef, G.S.E., Mitchum, G., Lukas, R., Niiler, P., 1999. Tropical Pacific near-surface currents estimated from altimeter, wind and drifter data. J. Geophys. Res. 104, 23313 – 23326. Latif, M., Barnett, T.P., 1994. Causes of decadal climate variabil-

155

ity over the North Pacific and North America. Science 266, 634 – 637. Latif, M., Roeckner, E., Mikolajewicz, U., Voss, R., 2000. Tropical stabilisation of the thermohaline circulation in a greenhouse warming simulation. J. Climate 13, 1809 – 1813. Liu, Z., 1999. Forced planetary wave response in a thermocline circulation. J. Phys. Oceanogr. 29, 1036 – 1055. Liu, X., Zhang, R.H., 1999. Propagation and mechanism of decadal upper ocean variability in the North Pacific. Geophys. Res. Lett. 26, 739 – 742. Meehl, G.A., Gent, P., Arblaster, J.M., Otto-Bliesner, B., Brady, E., Craig, A., 2001. Factors that affect amplitude of El Nino in global coupled climate models. Clim. Dyn. 17, 515 – 526. Neelin, J.D., Battisti, D.S., Hirst, A.C., Jin, F.F., Wakata, Y., Yamagata, T., Zebiak, S.E., 1998. ENSO theory. J. Geophys. Res. 103 (C7), 14261 – 14290. Nikias, C.L., Petropulu, A.P., 1993. Higher-Order Spectra Analysis: A Nonlinear Signal Processing Framework PTR Prentice Hall, Englewood Cliffs, NJ. Pierce, D.W., Barnett, T.P., Latif, M., 1999. Connections between the Pacific Ocean tropics and midlatitudes on decadal timescales. J. Clim. 13, 1173 – 1194. Roeckner, E., Oberhuber, J.M., Bacher, A., Christoph, M., Kirchner, I., 1996. ENSO variability and atmospheric response in a global atmosphere – ocean GCM. Clim. Dyn. 12, 737 – 754. Schneider, N., 2000. A decadal spiciness mode in the tropics. Geophys. Res. Lett. 27, 257 – 260. Schneider, N., Venzke, S., Miller, A.J., Pierce, D.W., Barnett, T.P., Deser, C., Latif, M., 1999. Decadal coupling of northern midlatitude and equatorial Pacific Ocean via the thermocline? Geophys. Res. Lett. 26, 1329 – 1332. Shil’nikov, L.P., 1965. A case of the existence of a denumerable set of periodic motions. Sov. Math. Dokl. 6, 163 – 166. Shil’nikov, L.P., 1970. A contribution to the problem of the structure of an extended neighborhood of a rough equilibrium state of saddle-focus type. Math. USSR Sbornik 10, 91 – 102. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., Farmer, J.D., 1992. Testing for nonlinearity in time series: the method of surrogate data. Physica D 58, 77 – 94. Timmermann, A., 1999. Detecting the Nonstationary Response of ENSO to Greenhouse Warming. J. Atmos. Sci. 56, 2313 – 2325. Timmermann, A., 2001. Changes of ENSO stability due to Greenhouse Warming. Geophys. Res. Lett. 28, 2061 – 2064. Timmermann, A., Jin, F.-F., 2002a. Homoclinic chaos: a paradigm for decadal ENSO amplitude modulations. J. Atmos. Sci. (in press). Timmermann, A., Jin, F.-F., 2002b. A nonlinear mechanism for decadal El Nin˜o amplitude changes. Geophys. Res. Lett. 29 (10.1029/2001GL013369). Timmermann, A., Latif, M., Gro¨tzner, A., Voss, R., 1999a. Modes of climate variability assimulated by the coupled atmosphere – ocean model ECHAM3/LSG, Part I: ENSO-like climate variability and its low-frequency modulation. Clim. Dyn. 15, 605 – 618. Timmermann, A., Latif, M., Bacher, A., Oberhuber, J., Roeckner, E., 1999b. Increased El Nin˜o frequency in a climate model forced by future greenhouse warming. Nature 398, 694 – 696.

156

A. Timmermann / Global and Planetary Change 37 (2003) 135–156

Timmermann, A., Voss, H.U., Pasmanter, R., 2001. Empirical dynamical system modeling of ENSO using nonlinear inverse techniques. J. Phys. Oceanogr. 31, 1579 – 1598. Torrence, C., Webster, P.J., 1998. Interdecadal changes in the ENSO-monsoon system. J. Climate 12, 2679 – 2690. Vautard, R., Yiou, P., Ghil, M., 1992. Singular-spectrum analysis: a toolkit for short, noisy chaotic signals. Physica D 58, 95 – 126. Voss, H.U., 2001. Analyzing nonlinear dynamical systems with nonparametric regression. In: Mees, A.I. (Ed.), Nonlinear Dynamics and Statistics. Birkh. Auser, Boston, pp. 413 – 434.

Voss, H., Kurths, J., 1997. Reconstruction of nonlinear time delay models from data by the use of optimal transformations. Phys. Lett., A 234, 336 – 344. Weaver, A.J., 1999. Extratropical subduction and decadal modulation of El Nino. Geophys. Res. Lett. 26, 743 – 746. Zebiak, S.E., Cane, M.A., 1986. A model El Nino-Southern Oscillation. Mon. Weather Rev. 115, 2262 – 2278. Zhang, R.H., Rothstein, L.M., Busalacchi, A.J., 1998. Origin of upper-ocean warming and El Nin˜o change on decadal scales in the tropical Pacific Ocean. Nature 391, 879 – 883.