(Un)predictability of strong El Niño events
John Guckenheimer, Axel Timmermann, Henk Dijkstra, Andrew Roberts
aa r X i v : . [ m a t h . D S ] O c t (Un)predictability of strong El Ni˜no events John Guckenheimer , Axel Timmermann , Henk Dijkstra , and Andrew Roberts Mathematics Department, Cornell University, Ithaca, NY 14853 IBS Center for Climate Physics, Pusan National University, Busan, South Korea Institute for Marine and Atmospheric Research Utrecht, Department of Physics andAstronomy, Utrecht University, Utrecht, The Netherlands Cerner Corporation, Kansas City, MissouriMay 9, 2019
Abstract
The El Ni˜no-Southern Oscillation (ENSO) is a mode of interannual variability in the coupled equa-torial Pacific coupled atmosphere/ocean system. El Ni˜no describes a state in which sea surface temper-atures in the eastern Pacific increase and upwelling of colder, deep waters diminishes. El Ni˜no eventstypically peak in boreal winter, but their strength varies irregularly on decadal time scales. There wereexceptionally strong El Ni˜no events in 1982-83, 1997-98 and 2015-16 that affected weather on a globalscale. Widely publicized forecasts in 2014 predicted that the 2015-16 event would occur a year earlier.Predicting the strength of El Ni˜no is a matter of practical concern due to its effects on hydroclimate andagriculture around the world. This paper discusses the frequency and regularity of strong El Ni˜no eventsin the context of chaotic dynamical systems. We discover a mechanism that limits their predictabilityin a conceptual “recharge oscillator” model of ENSO. Weak seasonal forcing or noise in this model caninduce irregular switching between an oscillatory state that has strong El Ni˜no events and a chaotic statethat lacks strong events, In this regime, the timing of strong El Ni˜no events on decadal time scales isunpredictable.
The problem in predicting El Ni˜no events is that they occur quite irregularly, and their development seems tobe different each time [20]. In addition to observations, the amplitudes of El Ni˜no events in a long simulationof the GFDL CM2.1 coupled ocean-atmosphere General Circulation model with constant forcing are highlyvariable on decadal time scales [31]. The predictability horizon of individual El Ni˜no events is thought tobe about 6-9 months, depending on the season and the phase of the ENSO cycle. During boreal springthe coupled ocean-atmosphere system is thought to be at its frailest state [30]. Then the system is mostsusceptible to perturbations [29] which leads to a ‘spring’ predictability barrier in April/May [18]. The roleof the initial error pattern has been emphasized and in particular its interaction with the seasonal cycle andthe internal ENSO cycle [21, 8, 34]. However, since detailed surface and ocean subsurface observations ofENSO began in the 1950s, very strong El Ni˜no events have occurred only once every 15-20 years (1982, 1997and 2015) and the processes controlling their predictability horizon are largely unknown.Similar to investigations of tipping points [24], we would like to identify precursors that predict thestrength of impending events and their frequency on both shorter and longer time scales. Processes deter-mining ENSO variability and its limits of predictability have been studied in conceptual models [25, 28, 16].Such models seek to identify key dynamical processes in the complex, coupled ocean-atmosphere system inthe tropical Pacific.This paper discovers a new dynamical mechanism that limits the predictability of strong El Ni˜no eventsin the low-dimensional conceptual Jin-Timmermann (JT) model of ENSO originally proposed by Jin [16]and then extended by Timmermann et al. [27]. The state space variables of this deterministic model are sea1urface temperatures of the equatorial western Pacific and eastern Pacific, and the thermocline depth (e.g.the depth of the 20 ◦ C isotherm) of the western Pacific. Nonlinear terms in the model are associated withthe wind-induced anomalous advection of sea surface temperature anomalies and thermocline dynamics thataffect these state variables. Comparisons of this reduced order model with observational data confirm thatthe model embodies the basic features of ENSO in the tropical Pacific [1].Recently, a multiple time scale analysis of the JT model was performed to gain new insight into itsdynamics [22]. This work transformed the model equations into an equivalent dimensionless system whoseequations are x ′ = ρδ ( x − ax ) + x ( x + y + c − c tanh( x + z )) y ′ = − ρδ ( ay + x ) z ′ = δ ( k − z − x ) , (1)The variables x, y and z of this model represent the sea surface temperature difference between the easternand western equatorial Pacific, the sea surface temperature of the western equatorial Pacific relative to anominal reference temperature, and the thermocline depth of the western Pacific, respectively. Each of thefive parameters δ, ρ, c, k and a represents a combination of physical characteristics of the tropical Pacific.The tuning of these parameters and the sensitivity of the El Ni˜no variability to these parameters is discussedin [22]. Before we present additional analysis of this model, we recall relevant background material fromdynamical systems theory. Strong El Ni˜no events are typically characterized by a threshold of average sea surface temperature anomalyin the eastern equatorial Pacific, e.g., the Nino 3.4 region (5N-5S, 170W-120W). [2] Within the setting ofthe JT model, we set the threshold x > − . x represents the difference insea surface temperature of the eastern and western tropical Pacific. The events are then recurrences to thephase space region x > − .
5. Predictability of strong El Ni˜no events in the JT model is then a matter ofthe frequency and regularity of recurrences to this region. We use a simple discrete time dynamical systemto illustrate the type of unpredictability that we subsequently find in the JT model.Sensitive Dependence to Initial Conditions [23] has been a key concept in the study of chaotic dynamicalsystems: pairs of nearby initial conditions yield trajectories that separate from each other. Statistical theoriesof chaotic attractors [33] emphasize asymptotic properties like invariant measures, Lyapunov exponents andentropies. There are dynamical systems in which multiple time scales are an emergent phenomenon. Sometrajectories spend long times near metastable invariant sets, occasionally making transitions between them.We use the term epoch to denote a maximal time interval when the system is close to one of the metastableinvariant sets.Our viewpoint is strongly influenced by the Wentzell-Freidlin theory of stochastically perturbed dynamicalsystems [10]. Imagine a deterministic dynamical system with several attractors and a small stochastic orrandom perturbation of this system. Wentzell and Freidlin introduce a Markov chain whose states arethe attractors of the deterministic system with transition probabilities given by the observed frequencies oftransitions between their basins of attractions in trajectories of the stochastic system. This same Markovmodel can be used for deterministic systems which have an attractor comprised of several almost invariant sets with infrequent transitions from one to another [11].We illustrate these concepts with iterations of the one dimensional map g ( x ) = α ∗ x (1 − x ) , x ∈ [ −∞ , ∞ ].When 2 < α < √ / ≈ . − ,
0] and [0 , α = 1 / √ α , g will have stable periodic orbit(s) that are its attractors,but for a positive measure set of α , the map has chaotic attractor(s). Figure 1 illustrates the dynamicswhen α = 2 .
6, seemingly a value that gives rise to a single chaotic attractor. The intervals [ g ( − / √ , , g (1 / √ iterates. The2 x -202 g ( x ) -2 -1 0 1 2 x -202 g ( g x )) t -202 x ( t ) run length l og ( f r equen cy ) a bc d Figure 1: Mode switching in a discrete one dimensional iteration: (a) Graph of the function g ( x ) = 2 . ∗ x (1 − x ). The intervals [ g ( − / √ ,
0] and [0 , g (1 / √ ± / √ g together with an expanded plot near 1 / √ x = 0 .
2. Epochsin which the sign of the iterates remain constant are clearly visible. (d) log(number of epochs of length n)against n in a trajectory segment of 10 iterates.longest epoch has length 610; the shortest has length 8. The approximately linear slope of the graph suggeststhat the epoch lengths can be modeled by a Poisson distribution. In a longer trajectory of 10 iterates, wefind that all pairs of integers ( m, n ) with 8 ≤ m, n ≤
60 occur as successive epoch lengths, indicating thatthe length of an epoch does not enable one to predict the length of the next epoch. Superficially, successiveepoch lengths seem to be statistically independent of each other.The unpredictability of epoch lengths in this model is related to its sensitive dependence to initialconditions. Nearby initial conditions yield trajectories that separate from one another so that informationabout their approximate location is lost after a moderate number of iterates. In this example, recurrencesto the intervals [ −∞ , −
1] and [1 , ∞ ] mark the transitions from one basin to the opposite one. The set ofinitial conditions where a transition occurs at the n th iterate is a finite union of intervals, but the number ofintervals increases geometrically with n and, for a positive measure set of α , their distribution approachesan absolutely continuous invariant measure [5].We also computed the epoch lengths for sample trajectories of the stochastic map g r ( x ) = 2 . ∗ x (1 − x )+0 . ∗ ξ (where ξ is a normally distributed random variable) to reside in one of the two attractor basins.Here, we are certain that the dynamics have an absolutely continuous invariant measure and exponentialdecay of correlations [3]. The resulting distribution of epoch lengths has no apparent differences from thoseshown in Figure 1(d) for the deterministic system with two almost invariant sets. Both the deterministicand random versions of this model have epoch lengths that are unpredictable from past history. Several different types of attractors are found in the JT model. We are especially interested in those thathave strong El Ni˜no events. Since the variable x represents the difference in sea surface temperatures betweenthe eastern and western tropical Pacific, these events are marked by small values of | x | . We used bifurcationanalysis in our search for different types of attractors. The program MATCONT [7] was used to locate pointsof Hopf bifurcation [14] and then to identify points of “generalized Hopf” bifurcations where the bifurcationis neither subcritical nor supercritical. Starting at parameters with a supercritical Hopf bifurcation andvarying parameter δ , periodic orbits were continued until they lost stability at a period doubling bifurcation3 a -0.4 8-0.6 -1 y a -27-1 x -3 -4b Figure 2: (a) (Blue) Hopf bifurcation and (red) period doubling curves in the ( δ, a ) plane. The point near(7 . , . δ, a ) parameter values (7 . , . ρ, c, k ] = [0 . , . , . a showing a curve of equilibria (black) that passes through aHopf bifurcation point (black dot) near δ = 7 . δ = 7 . δ, ρ, c, k ] = [0 . , . , . , . δ led to a regime with an MMO attractor. Decreasing δ showed that the ranges of δ with MMO and chaotic attractors overlapped. At the end of this section, wediscuss the bifurcations that occur on the boundaries of the overlap region.The parameters [ δ, ρ, c, k, a ] = [0 . , . , . , . , . x remains far from 0. As illustrated in Figure 4, this attractor can be analyzedby introducing a cross-section and studying its return map. Figure 4a shows the intersection of a trajectorywith the cross-section x = x eq , where x eq ≈ − . x at the equilibrium point. Black filledcircles in the figure are intersections of a periodic orbit with the cross-section. This periodic orbit appearsto be on the boundary of the basin of attraction of the attractor, and it plays a central role in bifurcationsof the attractor. The return map contracts in one direction, and stretches and folds in a second direction.Figure 4b shows the folding by plotting the value of z at each return to the cross-section against the value4 x z 1 y2 0-0.5-1-1.5 a 0 10 20 30 40 50t-4-3-2-10 x b Figure 3: (a) Phase portraits of a chaotic attractor (red) and an MMO periodic orbit (blue) that coexistat parameters [ δ, ρ, c, k, a ] = [0 . , . , . , . , . x coordinate of the two trajectories. y a z(t) z ( t + ) b Figure 4: (a) Intersections ( y n , z n ) of a trajectory in the chaotic attractor with the plane x = x eq ≈ − . z n +1 against z n for this trajectory.of z at the previous return. There is sensitivity to initial conditions within the chaotic attractor, but it isrelatively mild: nearby initial conditions separate along the attractor without leaving it.The basins of attraction of the two attractors are intertwined with each other in an intricate way, asillustrated in Figure 5. Trajectories starting in a 500 ×
500 initial condition grid in the cross-section x = x eq ≈ − . x = x eq ≈ − . x = x eq . These initial conditions were chosen alonga curve of length approximately 0 . y coordinate for the third return of points in aregion W similar to the region displayed in Figure 5. In the middle of this strip, there is a “ridge” of pointsfor which the y coordinate of the return is close to − .
1. These are points whose third return is associatedwith a strong El Ni˜no event. The sides of the ridge are points where the third iterate of the return maphas large stretching. This can be quantified by computing the finite time Lyapunov exponent or FTLE [15]of trajectories in W . The FTLE is the largest singular value of the variational equations of the flow mapalong a trajectory segment of specified length. It measures the maximal stretching of infinitesimally closetrajectories along the trajectory segment. Figure 7 plots this quantity on a log scale in a strip of initialconditions in the plane x = − . . ≈
5. Values are on abase 10 log scale. The center of the strip, indicated by its blue color, consists of points in the MMO basin ofattraction and their FTLEs are approximately 5. Trajectory segments along two flanking strips have FTLEtwo orders of magnitude larger, demonstrating the high sensitivity to initial conditions there.These properties contrast with the one dimensional map g discussed in the previous section. As theparameter a of that example increases to its critical value 3 √ /
2, the basins of two attractors touch ata single point and the two attractors merge into one. In the JT model, bifurcations of the MMO andchaotic attractors occur at different parameter values and the two attractors do not merge to become one.Apparently, when one of the invariant sets ceases to be an attractor, most points that were in its basin ofattraction have trajectories that now flow to the other attractor. We investigated how this happens as theparameter a is varied.In the case of the chaotic attractor of the JT model, a boundary crisis bifurcation [12] occurs where theattractor becomes a non-attracting invariant set. As a decreases to a value close to 7 . a , a chaoticinvariant set persists, but it is no longer attracting. Trajectories can “escape” from the former attractorbasin along the unstable manifold of the periodic orbit.The bifurcations associated with the MMO attractors also involve a boundary crisis. The MMO attractorfor a = 7 . − .
46 and 0 . ± .
47. Because the ratio of the negative eigenvalue to the real part of the complex eigenvalueshas large magnitude, a large volume of phase space is drawn toward the stable manifold of the equilibriumbefore their trajectories spiral away slowly along the two dimensional unstable manifold of the equilibrium.Moreover, trajectories starting close to the equilibrium in the unstable manifold increase their distance fromit by a factor of only about 1 . z1-3-2.5 y-2 2 x -1.5-1-0.5-1 -0.5 0 a Figure 6: (a): Ensemble of 100 trajectory segments in the cross section x = x eq ≈ − . . .
4. (b): The image of the y -coordinate under the third iterate of the return map of a thin strip of initial conditions (dark blue) ofwidth 0.035 in the cross section x = x eq . Points in the middle of the strip make a large amplitude excursionrepresenting a strong El Ni˜no. The map stretches the y -coordinate by a factor of order 10 where the returnshave intermediate amplitude between those on the ridge and the strip boundaries. z -0.8-0.75-0.7-0.65-0.6 y Figure 7: The largest finite time Lyapunov exponent for trajectory segments of length five. Initial conditionsare a 50 ×
50 grid in the plane x = − . x = − .
5, which we set (arbitrarily) as a thresholdfor strong El Ni˜no events. This value is larger than the maximum of x in the chaotic attractor (see Figure3.) Figure 8a plots intersections of the unstable manifold of the equilibrium with this cross-section as largeblue dots. There is a fold that occurs near z = 0 .
835 dividing the intersection into two branches that areclose together. To investigate the next return of the unstable manifold to the cross-section, we approximatethe intersection by a quadratic curve, plotted as black dots in the figure. The returns of the black pointsare plotted as green x’s. This image lies close to the intersection of the unstable manifold of the equilibriumwith the cross-section, so we regard the return map as approximately one dimensional. Figure 8b shows thegraph of this one dimensional map, parametrized by the z coordinates of domain and range. There is aninterval (roughly z ∈ [0 . , . p at ( y, z ) ≈ ( − . , . − , confirming that the return map is approximately rank 1. Since themagnitude of one eigenvalue is larger than 1, the periodic orbit containing p is not stable, precluding that itis the MMO attractor. Nonetheless, the MMO attractor lies close to this periodic orbit. As the parameter a varies, the one dimensional approximations of the return map to the cross-section z = 1 . a interval [7 . , . .8346 0.8354 0.8363 z -0.076-0.0758-0.0755 y a z in z r e t b in z r e t c a chaotic attractor d MMO attractor
Figure 8: (a) Intersections of the unstable manifold of the equilibrium with the cross-section x = − . z -coordinateplotted against their initial values on the black curve of subplot (a). The black line is z ret = z in . Theparameters ( δ, a ) = (0 . , . . , . δ, a ) = (0 . , . x = 1 . δ, a ) = (0 . , . δ, a ) = (0 . , . x = 1 . a increases [9]. At a = 7 . δ, a ) parameter space where there isbistability. Values of δ were selected and then a divide and conquer strategy was utilized with varying a . Inthe region with chaotic attractors, the distance between the attractor chaotic attractors and the periodic orbitshown in Figure 3 indicates the distance to the bifurcation boundary. The scale of Figure 3 does not allowone to see this distance, but it is clearly visible for values of a farther from the boundary. Newton’s methodwas applied to the return map in order to compute the periodic orbits with high precision. The extent of theattractor was determined by computing intersections of along trajectory segment with the cross-section. Oncrossing the bifurcation curve, we observed a discontinuous increase in the set of intersections of trajectoryand cross-section. For the MMO attractors, we computed approximate return maps like those displayed inFigure 8b,c. The boundary crisis occurs where the image of the minimum value of the return map is its fixed8 .48 1.5 1.52 1.54 1.56 1.58 1.6 1.62 1.64 1.66 z -0.95-0.9-0.85-0.8-0.75-0.7-0.65-0.6-0.55-0.5 y a in z r e t b Figure 9: (a): Two trajectories for parameter value a = 7 . z -coordinate under return maps for a curve of initial conditions that lie closeto the image of the return map. Blue x’s are images of the return map; magenta circles are images of thethird iterate. The third iterate has subintervals (shown by black boxes) that are mapped into themselves inaddition to a stable fixed point near 1.534.point of positive slope. This is similar to the bifurcation displayed for the one dimensional map g in Section2. In Figure 8c, the image of the minimum value is clearly larger than the fixed point, so that the map has ahorseshoe. Figure 8d displays three points on each of the two bifurcation boundaries.This subplot also showsthe parameter values ( δ, a ) = (0 . , . δ, a ) = (0 . , . a > . Mode Switching
The model of Jin and Timmermann [27] includes terms representing annual forcing from the seasonal cy-cle of solar insolation and stochastic forcing (e.g., due to rapidly varying wind stress on the sea surface).We found that the addition of either of these components to our simulations is capable of inducing modeswitching between the two attractors discussed in the previous section. Figure 10 illustrates a simulationof the periodically forced system. The parameter a was made a sinusoidally varying function of time withperiod one year and amplitude 0 .
500 1000 1500 2000 t - . x + . y + . z Figure 10: Time series of the function g ( x, y, z ) ≈ − . x + 0 . y + 0 . z on a portion of a trajectorycomputed with a = 7 . .
002 sin(1 . t ) and [ δ, ρ, c, k ] = [0 . , . , . , . g = 1 .
306 and g = 1 .
22 as described in the text. The oscillations of g lie between these thresholds during chaotic epochs.between two almost invariant sets that approximate the attractors of the unforced model. The figure showsa representative time series of 2000 years selected from a simulation of 100,000 years. To visualize the almostinvariant sets in a time series, we defined an “observable” g to be the projection onto the line orthogonalto the eigenspace of the unstable complex eigenvalues at the equilibrium point. The oscillations of g nearthe chaotic attractor have much smaller magnitude than those of the MMO attractor. The epochs duringwhich the trajectory is close to one of the two attractors are readily apparent in the time series of g . Theswitching between modes is quite abrupt, and the duration of the epochs is irregular. We sought to investi-gate the distribution and independence of these durations, but found it difficult to define precise phase spaceboundaries to use in segmenting the trajectory into epochs automatically.We set two thresholds for the function g and computed the times when g decreased through a threshold.The values of the thresholds were chosen to be intermediate between the range of oscillations of g in thechaotic and MMO attractors. They are plotted as black lines in Figure 10. Most of the chaotic attractorepochs lie between two successive events, with the upper threshold crossed at the beginning of the epochand the lower threshold crossed at the end of the epoch. Moreover, the duration of most of these epochs issubstantially larger than the roughly decadal duration between strong El Ni˜no events during their epochs.This prompted us to select pairs of successive events separated by a time duration longer than 15 years asbracketing chaotic epochs. In most cases, blue crosses on the upper threshold mark the onset of chaoticepochs, while those on the lower threshold mark the termination of the epochs. The complements of thesetime intervals were regarded as strong El Ni˜no epochs. In some instances, such as a time interval near t = 20200, there are shorter oscillations in the chaotic range that have not been identified as chaotic epochs.Figure 11(a) displays a histogram of the number of epochs of duration less than 400 years in bins of width10. (The longest epoch has a duration of approximately 540 years.) These durations are quite long comparedto those observed in simulations of large climate models [32]. However, we chose a very small amplitude forthe periodic component of a to make it easier to segment a trajectory into epochs. With larger magnitudeof the periodic variation of a , the epochs are shorter and our criteria for segmenting a trajectory into epochsbreaks down. The distribution of epoch durations is roughly exponential. Figure 11(b) is a scatter plotwhich shows pairs of successive epoch durations. Apart from higher density of shorter durations, there is noapparent pattern. We conclude that epoch durations and ENSO regime predictability in this JT model withsmall annual forcing are unpredictable [26]. Discussion
Reconstructions of past climate from historical data indicate that ENSO is highly variable with long periodsof larger and smaller variance on decadal time scales [4, 19]. Simulations of global climate models reachsimilar conclusions [32]. There are many potential explanations for why ENSO is so variable. We havedemonstrated here that interactions of only a few degrees of freedom in the highly reduced Jin-Timmermannmodel can produce unpredictability of strong El Ni˜no events and a new type of ENSO complexity on a10
100 200 300 400 epoch n epo c h n + b Figure 11: (a): Histogram of epoch durations smaller than 400, sorted into bins of width 10 years. Thenumber of occurrences in each bin are displayed on a log scale to highlight the resemblance to an exponentialfrequency distribution. (b): Scatterplot of pairs of successive epoch durations smaller than 400.decadal time scale. This is a stronger statement than saying that that ENSO is chaotic. It places a timescale on the sensitive dependence on initial conditions. More specifically, our results show that epochs ofrelatively regular cycles of strong El Ni˜no events like those observed during the past century might notcontinue. If they do cease, they may resume at any time.We compare this unpredictability of ENSO regimes with ensemble weather forecasting. Ensemble fore-casting is applied in situations where one expects sensitive dependence to initial conditions; i.e., nearbyinitial conditions yield trajectories that diverge from one another. This divergence happens gradually, soforecasting gives predictions with uncertainty that grows with time. Here, switching is abrupt and the timescale for predicting when a mode switch will take place seems to be comparable to the decadal time scaleof strong El Ni˜no events. Rapid divergence of nearby trajectories from one another is concentrated in smallregions of the state space of the model. If we understand where this divergence takes place, then we canexploit that knowledge in our forecasts. We can also use what we have learned about the complex dynamicsof the low dimensional JT model as a guide to gain deeper insight into the coupled oceanic and atmosphericdynamics of ENSO in large climate models. Whether the large scale dynamics of these large models have lowdimensional attractors and almost invariant sets displaying different ENSO dynamics remain unansweredquestions.
Acknowledgments:
Henk Dijkstra’s work was partially supported by a Mary Shepard B. Upson VisitingProfessor position at the College of Engineering of Cornell University, Ithaca, NY. He thanks for Prof. PaulSteen (Cornell University) for being his host and the many interesting discussions. Axel Timmermann wassupported by the Institute for Basic Science (project code IBS-R028-D1).
References [1] Soon-Il An and Fei-Fei Jin. Nonlinearity and asymmetry of enso.
Journal of Climate , 17(12):2399–2412,2004.[2] Anthony G. Bamston, Muthuvel Chelliah, and Stanley B. Goldenberg. Documentation of a highlyENSO related SST region in the equatorial pacific: Research note.
Atmosphere-Ocean , 35(3):367–383,1997.[3] Michael Benedicks and Lai-Sang Young. Absolutely continuous invariant measures and random pertur-bations for certain one-dimensional maps.
Ergodic Theory Dynam. Systems , 12(1):13–37, 1992.[4] K.M. Cobb, C. D. Charles, H. Cheng, and R. L. Edwards. El Ni˜no-Southern Oscillation and tropicalPacific climate during the last millenium.
Nature , 424:271–276, 2003.[5] Welington de Melo and Sebastian van Strien.
One-dimensional dynamics , volume 25 of
Ergebnisse derMathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)] . Springer-Verlag,Berlin, 1993. 116] Mathieu Desroches, John Guckenheimer, Bernd Krauskopf, Christian Kuehn, Hinke M. Osinga, andMartin Wechselberger. Mixed-mode oscillations with multiple time scales.
SIAM Rev. , 54(2):211–288,2012.[7] A. Dhooge, W. Govaerts, and Yu. A. Kuznetsov. MATCONT: a MATLAB package for numericalbifurcation analysis of ODEs.
ACM Trans. Math. Software , 29(2):141–164, 2003.[8] W. Duan, X. Liu, K. Zhu, and M. Mu. Exploring the initial errors that casue a significant “springpredictability barrier” for El Ni˜no events.
J. Geophys. Res , 114:C04022, 2009.[9] J.-P. Eckmann. Roads to turbulence in dissipative dynamical systems.
Rev. Modern Phys. , 53(4, part1):643–654, 1981.[10] Mark I. Freidlin and Alexander D. Wentzell.
Random perturbations of dynamical systems , volume 260 of
Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences] .Springer, Heidelberg, third edition, 2012. Translated from the 1979 Russian original by Joseph Sz¨ucs.[11] Gary Froyland and Michael Dellnitz. Detecting and locating near-optimal almost-invariant sets andcycles.
SIAM J. Sci. Comput. , 24(6):1839–1863, 2003.[12] Celso Grebogi, Edward Ott, and James A. Yorke. Crises, sudden changes in chaotic attractors, andtransient chaos.
Phys. D , 7(1-3):181–200, 1983. Order in chaos (Los Alamos, N.M., 1982).[13] Celso Grebogi, Edward Ott, and James A. Yorke. Chaos, strange attractors, and fractal basin boundariesin nonlinear dynamics.
Science , 238(4827):632–638, 1987.[14] John Guckenheimer and Philip Holmes.
Nonlinear oscillations, dynamical systems, and bifurcations ofvector fields , volume 42 of
Applied Mathematical Sciences . Springer-Verlag, New York, 1990. Revisedand corrected reprint of the 1983 original.[15] G. Haller. Distinguished material surfaces and coherent structures in three-dimensional fluid flows.
Phys. D , 149(4):248–277, 2001.[16] FF Jin. An equatorial ocean recharge paradigm for ENSO .1. Conceptual model.
J. Atmos. Sci. ,54(7):811–829, APR 1 1997.[17] M. T. M. Koper and P. Gaspard. The modeling of mixed-mode and chaotic oscillations in electrochemicalsystems.
J. Chem. Phys. , 96(10):7797–7813, 1992.[18] M. Latif, T. P. Barnett, M. A. Cane, M. Fl¨ugel, N. E. Graham, H. von Storch, J.-S. Xu, and S. E.Zebiak. A review of ENSO prediction studies.
Clim. Dynam. , 9:167–179, 1994.[19] S. McGregor, A. Timmermann, M. H. England, O. Elison Timm, and A. T. Wittenberg. Inferredchanges in El Nino-Southern Oscillation variance over the past six centuries.
Climate of the Past ,9(5):2269–2284, 2013.[20] Michael J McPhaden, Axel Timmermann, Matthew J Widlansky, Magdalena A Balmaseda, and Tim-othy N Stockdale. The Curious Case of the EL Ni˜no That Never Happened: A Perspective from 40Years of Progress in Climate Research and Forecasting.
Bull. Amer. Meteor. Soc. , 96(10):1647–1665,October 2015.[21] M. Mu, W. Duan, and B. Wang. Season-dependent dynamics of nonlinear optimal error growth andENSO predictability in a theoretical model.
J. Geophys. Res. , 112:D10113, 2007.[22] Andrew Roberts, John Guckenheimer, Esther Widiasih, Axel Timmermann, and Christopher K. R. T.Jones. Mixed-Mode Oscillations of El Ni˜no–Southern Oscillation.
J. Atmos. Sci. , 73(4):1755–1766,2016.[23] David Ruelle. Sensitive dependence on initial condition and turbulent behavior of dynamical systems.In
Bifurcation theory and applications in scientific disciplines (Papers, Conf., New York, 1977) , volume316 of
Ann. New York Acad. Sci. , pages 408–416. New York Acad. Sci., New York, 1979.1224] Marten Scheffer.
Critical Transitions in Nature and Society: (Princeton Studies in Complexity) . Prince-ton University Press, 2009.[25] MJ Suarez and PS Schopf. A delayed action model for ENSO.
J. Atmos. Sci. , 45(21):3283–3287, NOV1 1988.[26] Axel Timmermann and Fei-Fei Jin. Predictability of coupled processes.
Predictability of weather andclimate , pages 251–274, 2006.[27] Axel Timmermann, Fei-Fei Jin, and Jan Abshagen. A nonlinear theory for El Ni˜no bursting.
J. Atmos.Sci. , 60(1):152–165, 2003.[28] E Tziperman, MA Cane, SE Zebiak, Y Xue, and B Blumenthal. Locking of El Nino’s peak time to theend of the calendar year in the delayed oscillator picture of ENSO.
J. Climate , 11(9):2191–2199, SEP1998.[29] P. J. Webster. The annual cycle and the predictibability of the tropical coupled ocean-atmospheresystem.
Meteor. Atmos. Phys. , 56:33–55, 1995.[30] PJ Webster and S Yang. Monsoon and ENSO: Selectively interactive systems.
Q. J. R. Meteor. Soc. ,118:877–926, 1992.[31] Andrew T. Wittenberg. Are historical records sufficient to constrain ENSO simulations?
GeophysicalResearch Letters , 36(12):n/a–n/a, 2009. L12702.[32] Andrew T. Wittenberg, Anthony Rosati, Thomas L. Delworth, Gabriel A. Vecchi, and Fanrong Zeng.ENSO Modulation: Is It Decadally Predictable?
Journal of Climate , 27(7):2667–2681, 2014.[33] Lai-Sang Young. Ergodic theory of attractors. In
Proceedings of the International Congress of Mathe-maticians, Vol. 1, 2 (Z¨urich, 1994) , pages 1230–1237. Birkh¨auser, Basel, 1995.[34] Y. Yu, M. Mu, and W. Duan. Does model parameter error cause a significant spring predictabilitybarrier for El Ni˜no events in the Zebiak-Cane model.