Understanding the origin of extreme events in El Niño-Southern Oscillation
Arnob Ray, Sarbendu Rakshit, Gopal K. Basak, Syamal K. Dana, Dibakar Ghosh
UUnderstanding the origin of extreme events in El Ni˜no-Southern Oscillation
Arnob Ray , Sarbendu Rakshit , Gopal K. Basak , Syamal K. Dana , , and Dibakar Ghosh ∗ Physics and Applied Mathematics Unit, Indian Statistical Institute, Kolkata 700108, India Stat-Math Unit, Indian Statistical Institute, Kolkata 700108, India Department of Mathematics, Jadavpur University, Kolkata 700032, India Division of Dynamics, Technical University of Lodz, 90-924 Lodz, Poland (Dated: June 23, 2020)We investigate a low-dimensional slow-fast model to understand the dynamical origin of El Ni˜no-Southern Oscillation. A close inspection of the system dynamics using several bifurcation plotsreveals that a sudden large expansion of the attractor occurs at a critical system parameter via atype of interior crisis. This interior crisis evolves through merging of a cascade of period-doublingand period-adding bifurcations that leads to the origin of occasional amplitude-modulated extremelylarge events. More categorically, a situation similar to homoclinic chaos arises near the criticalpoint, however, atypical global instability evolves as a channel-like structure in phase space ofthe system that modulates variability of amplitude and return time of the occasional large eventsand makes a difference from the homoclinic chaos. The slow-fast timescale of the low-dimensionalmodel plays an important role on the onset of occasional extremely large events. Such extremeevents are characterized by their heights when they exceed a threshold level measured by a mean-excess function. The probability density of events’ height displays multimodal distribution with anupper-bounded tail. We identify the dependence structure of interevent intervals to understand thepredictability of return time of such extreme events using autoregressive integrated moving averagemodel and box plot analysis.
PACS numbers: 05.45.-a, 92.70.Gt
I. INTRODUCTION
El Ni˜no-Southern Oscillation (ENSO) is one of thepowerful climatic and oceanic events and it is a result ofinterannual climate variability due to an interaction be-tween atmospheric and ocean circulations [1]. It has twophases: El Ni˜no, a warming phase of sea surface temper-ature (SST) in the eastern and central equatorial Pacificocean and La Ni˜na, an occasional cooling of ocean sur-face waters in that area. El Ni˜no usually occurs once inthe range of 2-7 years [2]. The key features of ENSO areits irregularity of occurrence and amplitude modulationdue to ocean-atmospheric instability [3, 4], which drawattention of the scientific community and planners overthe last few decades for its catastrophic effect and socioe-conomic impact. SST of eastern Pacific ocean is cold inthe normal period along the west coast of South America,but sudden changes of wind patterns, ocean circulationpatterns, and rise of the oceanic surface temperature ofeastern Pacific are observed during the El Ni˜no phase ofENSO episodes [5]. It becomes necessary to understandthe extreme nature of ENSO from a purely dynamicalsystem perspective.Strong El Ni˜no events have devastating effects on na-ture [6]. During El Ni˜no periods, normal patterns oftropical precipitation are heavily disturbed and ENSOalso incurs global climate change [7]. Drought may occurin Australia, Indonesia, India, Kenya, Morocco, Alge-ria, Canada, Mexico, etc., and floods in the central and ∗ Electronic address: [email protected] eastern Pacific regions, parts of South America close toArgentina, Chile, Peru, Ecuador, etc., during this phase.El Ni˜no significantly disrupts the local economy relatedto ecosystems, fisheries, and agriculture.Enormous research efforts have been exerted on de-veloping an effective dynamical model dealing with theatmosphere-ocean interactions, exploring the decadaltropical climate variability [4, 8, 9]. A mechanism ofpositive feedback and delayed negative feedback of theocean-atmosphere interaction in the equatorial Pacific,illustrates the dynamic features of the ENSO [10, 11].A coupled atmosphere-ocean model has been studied toobtain self-sustained oscillations without anomalous ex-ternal forces for understanding the characteristic featureof the ENSO phenomenon. A low- dimensional modelbased on recharge mechanism was tried to explain theexistence of decadal ENSO amplitude modulation anddecadal changes in the tropical mean state [9]. Thelow-order nonlinear ENSO model characterizes a strongasymmetry which is an intrinsic characteristic propertyof the ENSO phenomenon influenced by the nonlineareffect of the subsurface temperature [3].A kind of irregular spiking behavior has been reportedin the low-order ENSO dynamical model and more cate-gorically, in the dynamical sense, a homoclinic connectionto a saddle focus was found, which was responsible forthis kind of complexity in dynamics [9]. This explainedthe decadal recurrence of El Ni˜no as chaotic mixed-mode oscillation (MMO) [12] switching irregularly be-tween large- and small-amplitude oscillations. Althoughit was expected that predictability of MMO kind of ElNi˜no and La Ni˜na situations could be possible, it wasshown, using a modified low-dimensional slow-fast model a r X i v : . [ n li n . C D ] J un that strong El Ni˜no events on decadal timescales werecompletely unpredictable [13]. Still researchers are insearch of appropriate methods to address the question ofpredictability of El Ni˜no events to mitigate its harmfulimpact on the society [14–18].In this paper, we investigate the low-dimensional slow-fast ENSO model proposed by Timmermann et al. [9]to understand the origin of self-terminating occasionallarge spiking in the time evolution of the SST leadingto a type of extreme events mainly from the dynamicalsystem point of view. We have adopted numerical aswell as analytical techniques to extend the earlier results[9, 12] to develop further insights on the dynamical fea-tures of the system and explore possibility of extremelylarge events in the system with a parameter variation.We acknowledge the earlier observation [9] that a ho-moclinic connection plays a crucial role in the origin ofalmost decadal recurrence of El Ni˜no events of almostequal height, which we classify here, in our study, aschaotic MMOs. Additionally, we observe intermittentlarge spiking events in the temporal evolution of SSTwith variations in both return time and amplitude neara critical value of the bifurcation parameter. Such a kindof recurrent events, occurring rarely, is considered hereas extreme events which has been focused here from thedynamical system perspective [19–22]. Underlying mech-anism of the extreme events, in our model, is a type ofinterior crisis [23, 24] that occurs via a collision of theperiod-doubling (PD) cascades and period-adding (PA)bifurcations at a critical system parameter and it is dif-ferent from the typical interior crisis elaborated, in theliterature [25–28], as a collision of a chaotic orbit anda saddle orbit or saddle point. Experimental evidenceof extreme rogue waves has also been studied in nonlin-ear optics [20, 29]. A closer inspection, in the region ofthis merging of PD cascade and PA bifurcation, reveals aglobal instability that exists in the form of a channel-likestructure [24, 30, 31] in state space of the system. Thevolume of this channel introduces a variability of ampli-tude and return time of intermittent large spiking events,which we classify here as extreme events and are clearlydifferent from the El Ni˜no events, which information ismissing in the earlier report [9, 12]. We present numericalevidence that the slow-fast timescale of the model systemplays a crucial role in the onset of such extreme eventsin the model. Extreme events are identified here using amean excess function [32] and thereby, we also estimatethe probability density function (PDF) of events’ heightas usually done for extreme events, in general, in sin-gle or coupled systems [20, 28, 33–35]. We confirm rareoccurrence of this new kind of extreme event (large vari-ability of height and return time) using their interval ofreturn time. We have collected interevent interval datafrom numerical experiments and plotted histograms ofreturn time and then try to fit with the known distribu-tions, Weibull, Gamma, and Log-normal distributions fora comparison and thereby to find an appropriate distribu-tion of events. Besides, we try to address the question of predictability of extreme events using autoregressive in-tegrated moving average (ARIMA) model for time-seriesforecasting [36] and box plot analysis [37] of intereventinterval of extreme events.This paper is organized as follows. In Sec. II, we brieflydescribe the slow-fast ENSO model. In Sec. III, we ex-plain the system dynamics and try to understand themechanism of extreme events in Sec. IV. In Sec. V, wecharacterize the extreme events elaborately using statis-tical tools and address the question of predictability ofthe observed extreme events in Sec. VI. Finally, discussan important point in Sec. VII that the slow-fast param-eter plays a crucial role on the onset of extreme events.Results are summarized in Sec. VIII. II. SLOW-FAST ENSO MODEL
A general formulation of a slow-fast dynamical system[38] can be written as η dxdt = f ( x, y, λ ); dydt = g ( x, y, λ ) , (1)where x ∈ R m and y ∈ R n are fast and slow variables, re-spectively, λ ∈ R p is a model parameter, and 0 < η << η =0, the tra-jectory of Eq. (1) converges to the solution of differentialalgebraic equation f ( x, y, λ ) = 0 and dydt = g ( x, y, λ ),where S = { ( x, y ) ∈ R m × R n | f ( x, y, λ ) = 0 } is a criticalmanifold. Now we introduce the slow-fast ENSO modelproposed by Timmermann et al. [9]. If h is the ther-mocline depth of the western Pacific, then its evolutionequation is, dh dt = r (cid:18) − h − bLw s (cid:19) , (2)where r represents the dynamical adjustment timescale, w s denotes the zonal wind stress, b is the efficiency of w s in driving thermocline slope, and L is the basin width.The thermocline depth h of the eastern Pacific is relatedto h by the relation h = h + bLw s .The dynamical equations of equatorial sea surface tem-peratures T and T of the western and eastern Pacific,respectively, are represented by dT dt = − α ( T − T r ) − u ( T − T ) L ,dT dt = − α ( T − T r ) − w (cid:0) T − T sub ( h , T , T ) (cid:1) H m , (3)where α measures a typical thermal damping timescale, T r is the thermal relaxation towards a radiative-convective equilibrium temperature, H m denotes depthof the mixed layer, u and w are the zonal advection veloc-ity and equatorial upwelling velocity, respectively. T sub is the subsurface temperature being upwelled into themixed layer.The zonal wind stress is related to SST as w s = − µ ( T − T ) β , where β and µ are the coupling coefficientbetween SST and wind stress. Hence the expression of h becomes h = h + bLµ ( T − T ) β . (4)The zonal advection velocity u and equatorial upwelling velocity w are assumed to be proportional to the zonalwind stress anomalies w s as u L = (cid:15)βw s ; wH m = − ζβw s . (5)Here (cid:15) and ζ , respectively, quantify the strength of zonaland vertical advection. The expression of T sub can bewritten as T sub ( h , T , T ) = T r −
12 ( T r − T r ) (cid:34) − tanh (cid:16) H + h + bLµ ( T − T ) β − z (cid:17) h ∗ (cid:35) . (6)Here T r and H , respectively, are the mean easternequatorial temperature and eastern thermocline referencedepth, z measures the depth at which upwilling veloc- ity w takes its characteristic value, and h ∗ is sharpnessof thermocline. Hence we obtain the three-dimensional(3D) slow-fast model, dh dt = r (cid:18) − h − bLµ ( T − T )2 β (cid:19) , (7a) dT dt = − α ( T − T r ) − (cid:15)µ ( T − T ) , (7b) dT dt = − α ( T − T r ) + ζµ ( T − T ) (cid:32) T − T r + 12 ( T r − T r ) (cid:20) − tanh (cid:16) H + h + bLµ ( T − T ) β − z (cid:17) h ∗ (cid:21)(cid:33) , (7c)when h represents a slow variable for r < T and T are fast variables. Throughout this paper, we makea choice of fixed parameters [9], ζ =1.3, r = day − , α = 1 /
180 day − , T r = 29 . ◦ C, L = 15 × m, H m = 50m, bLµβ = 22 mK − , µ = 0 . K − day − , H = 100 m, h ∗ = 62 m, T r = 16 ◦ C, and z = 75 m. We analyze thelinear stability of the equilibrium points of system (7)and calculate the invariant manifold as detailed in theAppendix A. The system has one saddle (0 , T r , T r ) and aninterior equilibrium point ( h ∗ , T ∗ , T ∗ ), which is a saddlefocus. We integrate the system (7) using fourth-orderRunge-Kutta algorithm and fixed time step 0 .
01. Ourmain emphasis is to investigate the complex dynamicalbehavior of the system (7) by varying the strength ofzonal advection (cid:15) . III. COMPLEX DYNAMICS OF ENSO MODEL
Emergence of El Ni˜no and extreme events is shownin Fig. 1 with a series of temporal dynamics and their phase portraits in a range of (cid:15) ∈ [0 . , . (cid:15) = 0 . (cid:15) = 0 . T remains bounded ( T max <
24) here, however, for a smallincrease beyond (cid:15) ≈ . T max >
24) start appearing as shown in Fig. 1(e)for (cid:15) = 0 . Time T
26 28 30 T T ( a ) ( b ) Time T
26 28 30 T T ( h )( g ) Time T
26 28 30 T T ( j )( i ) Time T
26 28 30 T T ( l )( k ) FIG. 1: Evolution of different dynamics in ENSO model with a variation of (cid:15) . Temporal evolution of T and correspondingphase space diagram in ( T , T ) plane: [(a) and (b)] small-amplitude period-1 limit cycle at (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . ) at (cid:15) = 0 .
15, and [(k) and (l)] large-amplitude period-1 limit cycle at (cid:15) = 0 . highly irregular interval of time. The time spent duringspiraling out varies and it depends on how close the rein-jected trajectory reaches a vicinity of the saddle focusand thereby it makes a large variation in the number ofsmall oscillations that makes irregular return time of thelarge events. A tendency to develop a homoclinic chaoswith a local instability of the saddle focus is seen here,but the typical global stability of homoclinic chaos [39] isnever achieved due to the presence of a channel like struc-ture. The trajectory revolves around the saddle focus forquite some time due to local instability (attracted alongthe stable eigendirection, pushed away along the unstableeigenplane of the saddle focus) while a globally instabil-ity due to the channel-like structure induces a variationof the return path (global excursion) of the trajectorymaking a wide variation in both the amplitude and re-turn time of large spikes. Our main focus is the originof this exceptional case of large variation in return timeand amplitude of extreme events illustrated in Figs. 1(e)and 1(f), which has not been reported so far, to the bestof our knowledge. For higher values of (cid:15) = 0 . (cid:15) we notice atransition from chaotic MMOs to periodic MMOs (1 )shown in Figs. 1(i) and 1(j) for a larger (cid:15) = 0 .
15. Theperiod-1 limit cycle returns for larger (cid:15) , but with a largeramplitude of oscillation as shown in Figs. 1(k) and 1(l).
IV. ORIGIN OF EXTREME EVENTS
For a closer inspection of the dynamical evolution ofextreme events with both local and global instabilitiesin the system, a bifurcation diagram of T max against (cid:15) is drawn in Fig. 2(a). The stable interior equilibrium(FP) evolves into a limit cycle (LC) via Hopf bifurca-tion (HB) and it continues until (cid:15) ≈ . (cid:15) (a broader scenario in param- FIG. 2: (a) Bifurcation diagram of maximum peak values of T against (cid:15) . It shows different kind of oscillatory behavior ofthe system. (b) Period-parameter bifurcation: Period of periodic mixed mode oscillations are plotted with intervals. Zoomedversion of the bifurcation diagram (c) for (cid:15) ∈ [0 . , . (cid:15) ∈ [0 . , . eter space presented in Fig. 12 in Appendix A). At acritical value (cid:15) ≈ . T max suddenly increases to alarge value where the cascading sequence of PD mergeswith a PA sequence. The bounded chaotic mode explodeshere when we notice intermittent large-amplitude spikingoscillations with a variation in amplitude intercepted bya varying number of small oscillations as shown in thetemporal evolution of T in Fig. 1(e). We classify suchintermittent large events as the extreme events, whichare followed by the typical chaotic MMOs identified asEl Ni˜no events [12] and MMOs as shown in Fig. 2(c) [azoomed version of Fig. 2(a)] for increasing (cid:15) ; exemplarytemporal evolution of chaotic MMO and periodic MMOare shown in Figs. 1(g) and 1(i), respectively. For larger (cid:15) ,the dynamics return to LC, but with a large-amplitude asshown in Fig. 1(k). The periodic MMO sequence is bettervisualized from the right side of the bifurcation diagramin Fig. 2(a) when we decrease (cid:15) and observe a Farey se- quence L s [39] consisting of large-amplitude oscillations(denoted here by L = 1) followed by small-amplitudeoscillations (denoted by s = 1 , , , ..., ∞ ). The Fareysequence of periodic MMOs emerges in successive pa-rameter windows intermediate to chaotic windows of thecontrol parameter (cid:15) . Figure 2(b) shows the sequence ofMMOs (1 − − ... − ) in a period-parameter bifur-cation plot. It first emerges with a period-1 limit cycle(1 ) from (cid:15) ≈ . MMO that shows a devil’s staircase [39] with an asymp-totic increase in the time period of oscillations. Fromour numerical simulations, we are able to record a MMOwith a maximum number of s = 16 (with our best ef-fort) near (cid:15) = 0 .
1, although MMOs with larger s valuesreally exist in the system. The parameter window of (cid:15) for each MMO becomes narrower with increasing s andthe time period of MMOs increases asymptotically. Theinterval between any two successive MMO windows in (cid:15) [Fig. 2(b)] is a chaotic window and hence the MMOsemerge in a sequence of alternate PA bifurcations. As wedecrease (cid:15) , we find an increasing number of small oscil-lations s = 1 , , , ...,
16 in the periodic MMO windowsuntil the sequence collides with bounded chaos.To make a clear picture of this scenario, a small rangeof (cid:15) values is zoomed in two separate bifurcation dia-grams in Figs. 2(c) and 2(d). Figure 2(c) shows an ex-ample of PA sequence in the range of (cid:15) ∈ [0 . , . (cid:15) until it reaches a chaotic MMO(1 ∞ ) and arrive at a critical point of transition wherethe amplitude drops down to small-amplitude boundedchaos at (cid:15) ≈ . (cid:15) ∈ [0 . , . (cid:15) . This sudden change in amplitude of chaosoccurs near (cid:15) ≈ . (cid:15) value as shown in Fig. 1(g) and as re-ported earlier [12]. For a macroscopic view, we draw thecritical manifold of the slow-fast system and search forthe regions of instabilities in 3D phase space that pro-vides a minefield for the generation of extreme events.The critical manifold is obtained by solving a truncatedsystem of the fast variables T and T , given by the Eqs.(7b) and (7c), S = (cid:40) ( h , T , T ) ∈ R (cid:12)(cid:12)(cid:12)(cid:12) α − µ (cid:32) (cid:15) ( T − T ) + ζ (cid:16) T − T r + 12 ( T r − T r ) (cid:104) − tanh H + h + bLµ ( T − T β − z h ∗ (cid:105)(cid:17)(cid:33) = 0 (cid:41) . (8)The Jacobian J S of the truncated system is now de-rived at each point on the critical manifold S defined by using Eq. (8). The Jacobian J S = (cid:2) J ij (cid:3) × , where J = − α + 2 (cid:15)µ ( T − T ); J = − (cid:15)µ ( T − T ); J = − ζµ (cid:32) T − T r + 12 ( T r − T r ) (cid:20) − tanh (cid:16) H + h + bLµ ( T − T β − z (cid:17) h ∗ (cid:21)(cid:33) + ζµ ( T − T ) (cid:32) bLµ h ∗ β ( T r − T r ) (cid:20) sech (cid:16) H + h + bLµ ( T − T β − z (cid:17) h ∗ (cid:21)(cid:33) ; J = − α + ζµ (cid:32) T − T r + 12 ( T r − T r ) (cid:20) − tanh (cid:16) H + h + bLµ ( T − T β − z (cid:17) h ∗ (cid:21)(cid:33) + ζµ ( T − T ) (cid:32) − bLµ h ∗ β ( T r − T r ) (cid:20) sech (cid:16) H + h + bLµ ( T − T β − z (cid:17) h ∗ (cid:21)(cid:33) . (9)We plot the critical manifold [38] in a 3D surface inFig. 3 (green) and mark separate sections (dashed lines)of attracting ( S a ), repelling ( S r ) and saddle ( S s ) regionsof the critical manifold using the sign of the eigenval-ues of J S at each point of the surface. The saddle focus(76 . , . , .
33) (red circle) is lying on the border ofthe attracting and the repelling regions of the criticalmanifold. Starting from a precrisis point at (cid:15) = 0 . S r ) of the critical manifold. In contrast,Fig. 3(b) shows a postcrisis chaotic attractor on the criti-cal manifold S for (cid:15) = 0 . (cid:15) parameter near the critical point when the trajec-tory is destined to travel a close vicinity of the saddlefocus (red circle) in the repelling region. It takes longertime to spiral out since the trajectory gets slower nearthe saddle focus and makes many small-amplitude os- FIG. 3: Critical manifold S (green surface). S a , attracting part; S r , repelling part; and S s , saddle part of the critical manifold. T = T plane in brown, and a red circle indicates the position of saddle focus (76 . , . , . (cid:15) =0.0984 (precrisis) and (b) expanded attractor (blue line) for extreme events (cid:15) =0.0985 (postcrisis). cillations before moving out for a global excursion, butit is reinjected to reach another close proximity of thesaddle focus along the stable eigendirection. This is ex-actly the way a homoclinic chaos originates. However, incontrast to homoclinic chaos, the trajectory is no moreglobally stable rather traces different trajectories duringthe re-injection period due to the presence of a global in-stability in the form of a channel-like structure. Therebyit originates large spikes with a variability of amplitudeand return time intervals that depend upon the numberof small oscillations.Now we explain here the role of the channel-like struc-ture that creates a global instability. The geometric sin-gular perturbation theory [40] explains the existence of alocally invariant slow-manifold S ( r ) close to the criticalmanifold S within a domain O ( r ), where attracting orrepelling properties of the critical manifold remain un-changed. This invariant manifold consists of stable (at-tracting) and unstable (repelling) slow manifolds [41]. Asusual a slow increase in the size of a bounded chaotic at-tractor occurs with the increase of (cid:15) and it enters therepelling part of the critical manifold at a critical valueof (cid:15) . After entering the repelling section, the trajectorytends to move toward the saddle focus along the sta-ble manifold as mentioned above. When the trajectoryreaches a close vicinity of the saddle focus, it spirals outalong the 2D unstable manifold and attempts a large ex-cursion to originate a large spike, but passes through achannel-like structure. The global stability of homoclinicchaos with a fixed global trajectory, as expected in sucha situation, is broken by the channel-like structure. Thisstructure appears due to changes in the alignment of themanifolds in phase space as the control parameter (cid:15) isvaried [24, 42, 43]. To illustrate the influence of thischannel-like structure, we follow the trajectory micro-scopically. We project collected data on a segment ofthe trajectory as shown in Fig. 4(a) from a long timeseries of large events shown in Fig. 1(e). The trajectory while spiraling out, passes through a very narrow chan-nel in phase space. The system trajectory escapes fromits small-amplitude oscillation while spiraling near thesaddle focus and passes through this leaky channel.This instability region explicitly prescribes the locationof the leak in phase space and most importantly, its vol-ume is inversely proportional to the rarity of large spikingevents. The trajectory spends longer time in the narrowchannel when the number of small oscillations are largerand, thereby it creates longer interval of the large eventsand also a wide spreading in the reinjection path of thetrajectory that makes an amplitude variation of the largespikes or events. The role of the channel is further illus-trated by increasing (cid:15) = 0 . (cid:15) )occurs at a crisis point [23] that makes a sudden changein the size of attractor. An additional global instabilityof trajectories is created by a channel-like structure inphase space. The channel volume is modulated by thevariation of (cid:15) and as a result, two kinds of events, theextreme events and the typical El Ni˜no events, evolve asshown in Figs. 1(e) and 1(g), respectively, for two dif-ferent values of (cid:15) , one with large variability in amplitudeand return time making its rare occurrence and anotherone with almost identical heights and more frequent oc-currence of events. FIG. 4: Variation of the channel structure. Segments of a trajectory (blue lines) plotted for (a) (cid:15) =0.0985 (a narrow channel)and (b) (cid:15) =0.1076 (a wider channel) on the critical manifold S (green surface). Red circles are the positions of saddle focus(76 . , . , .
33) for (a) and (73 . , . , .
52) for (b).
V. EXTREME EVENTS: STATISTICALPROPERTIES
We explore here the statistical properties of the ex-treme events, the probability distributions of eventheights and interevent intervals. At first, a mean ex-cess function is defined to identify a threshold level forqualifying extreme events.
A. Mean excess function: Extreme event qualifier
We use a mean excess function [32, 44] to define athreshold of an event. Any event larger than this thresh-old is declared as extreme. For a random variable X , itis defined as e ( u ) = E ( X − u | X > u ) , (10)where E ( X ) denotes expectation of the random variable X and u is the threshold value which varies betweenthe minimum and maximum values of the observed data.This function gives an expected value of excess of a ran-dom variable over a certain threshold. Then we plot thefunction e ( u ) with varying u and there we estimate athreshold value u ∗ until this functional relation is a bet-ter fit to a straight line. If any event crosses u ∗ , thenthe event is considered as an extreme event. The ex-treme event qualifier threshold value is also obtained byconsidering few times more standard deviation than themean value [30, 45], deduced by using Peak Over Thresh-old approach [46].As an example, for (cid:15) = 0 . e ( u ) increases monotonically (blueline) with decreasing u and reaches a peak, and finallydecreases as seen in Fig. 5. However, a best straightline fit (black line) is possible until u = u ∗ = 25 .
5, be-yond which the plot deviates slowly in the beginning andthen largely (open blue circles) for lower u values fromthe straight line. We mark u ∗ = 25 .
23 24 25 26 27 28 u -101234 e ( u ) =0.0985=0.0986=0.0987=0.0988u * = 25.5 FIG. 5: Mean excess function plot against u for different val-ues of (cid:15) : (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) . Vertical dashed line is the threshold value u ∗ = 25 . (cid:15) = 0 . event qualifier (solid blue circle). Similarly, we obtainextreme event qualifiers as u ∗ = 25 . . . (cid:15) = 0 . , . . u ∗ increases (shifts to theright) with increasing (cid:15) and the rarity of occurrencesof extreme events gradually decreases. The slopes ofstraight line fits are − . − . − . − . (cid:15) = 0 . , . , . .
20 24 28
Event height -4 -2 P D F
20 24 28
Event height -4 -2 P D F
24 26 28 u R T ( u )
24 26 28 u R T ( u ) Time T Time T ( a ) ( b )( c ) ( d ) FIG. 6: Top: Probability density function of the event height in semi-log scale, inset figure shows the corresponding time series.Bottom: Corresponding return times: [(a) and (c)] (cid:15) = 0 . (cid:15) = 0 . B. Probability distribution of event heights
Now, probability density functions of all the peaks ina long run of time evolution of T are drawn in semi-log scale in Figs. 6(a) and 6(b) for two different valuesof (cid:15) , i.e. , (cid:15) = 0 . (cid:15) = 0 . . × iterations) so that the distributions are saturated. Ver-tical lines (red) indicate the threshold level of extremeevents estimated by the mean excess function (Fig. 5).Clearly, both PDFs show multimodal non-Gaussian dis-tributions with extreme events larger than the thresholdlevel. Figure 6(a) shows two dominant modes with anadditional less dominant mode (marked by a rectangularbox) for (cid:15) = 0 . (cid:15) push theadditional mode to the right (rectangular box) and PDFapproaches a bimodal distribution as shown in Fig. 6(b)at (cid:15) = 0 . T . One mode corresponds to small amplitude oscilla-tions and another mode corresponds to large amplitudeoscillations. Both the time signals of T (insets) show analternate large-amplitude spiking oscillations (related toextreme events) in alternate time sequence of spiralingsmall-amplitude oscillations. However, the height andreturn of large events show wide variability.Once the intermittent large events are identified as ex-treme from a simulated time series, which are above athreshold u ∗ . Now we introduce a function which is capa-ble to capture the spatial status as well as temporal fluc- tuation of occurrence of the events. For this, we estimatethe return time by taking an average of time-intervals ofall the events at certain height, say, u , and define it by RT ( u ), RT ( u ) = N N (cid:80) i =1 IEI i ( u ) , (11)where IEI i ( u ) is the i -th interevent interval correspond-ing to all the events of the fixed height u , and N is thenumber of interevent intervals collected from a long timeseries. Figure 6 (lower panel) shows the return time withrespect to the event height u for two exemplary valuesof (cid:15) used in the upper panels. We notice that the PDFsin Figs. 6(a) and 6(b) and corresponding return times inFigs. 6(c) and 6(d) are inversely proportional in naturein the interval 24 ≤ u ≤
28. So a larger return time ofa particular event u indicates that extreme events occurrarely. For both the cases, the return time of events in-creases with decreasing probability of number of events.Extreme events occur due to dynamical instability,which are responsible for the nonuniform behavior of thetail of the distribution. PDFs have different character-istics in different parameter regime and the tail of thedistribution cannot be easily fitted with the well-knownstatistical functions [47]. We attempt here to find thecharacteristic feature of the tail of our observed PDFs inFigs. 6(a) and 6(b) with the help of the shape parame-ter in a generalized Pareto distribution (GPD) [32]. The0 -0.06-0.04-0.020 A C F IEI N o . o f o cc u rr en c e s -5 OriginalWeibullGammaLog-normal
IEI N o . o f o cc u rr en c e s -4 OriginalWeibullGammaLog-normal -0.06-0.04-0.020 A C F -0.0100.01 1 10 20 -0.0100.01 1 10 20 ( a ) ( b ) ( d )( c ) FIG. 7: Auto-correlation function (ACF) with varying lag parameter ( τ ) for (a) (cid:15) = 0 . (cid:15) = 0 . { IEI n } and { IEI n } for n = 1 , , , . . . . PDFs of IEI are plotted (blue bars) in (c) (cid:15) = 0 . (cid:15) = 0 . probability density function of the Pareto distribution is G ξ,γ ( x ) = (cid:40) γ (cid:0) ξxγ (cid:1) − ξ − , ξ (cid:54) = 0 γ exp( − xγ ) , ξ = 0 (12)where x ≥ ξ is nonnegative, and 0 ≤ x ≤ − γξ otherwise. Here γ > ξ respectively, are the scaleand shape parameters of the distribution. The sign ofthe shape parameter ξ delineates the nature of the dis-tribution at the tail. In particular, if ξ ≥
0, then thedistribution has no upper limit, while ξ < e ( u ) against u curve by a straight line, when we can express the meanexcess function [44] for the case of GPD as e ( u ) ≈ γ − ξ + ξ − ξ u , (13)where ξ − ξ is the slope of the best fitted straight line.For parameter values (cid:15) = 0 . . − . − . − . − . ξ , tell us that the distributions haveupper-bounded tails, which are shown in Figs. 6(a) and6(b). C. Probability distribution of interevent interval
We study here the distribution and dependence struc-ture of interevent interval (IEI). Figures 7(a) and 7(b)depict that the autocorrelation functions (ACF) [48] ofIEI are out of small band ( i.e. , significant, taking 95%confidence interval) only for the first lag at (cid:15) = 0 . (cid:15) = 0 . { IEI n : n = 1 , , , . . . } is an uncorrelated process[see the inset figure in Fig. 7(a)] for the first case. Alsothe process { IEI n : n = 1 , , , . . . } is an uncorrelatedprocess [see the inset figure in Fig. 7(b)] for the secondcase. The similar comment is applicable for the process { IEI n − : n = 1 , , , . . . } [similarly to that of Fig. 7(a)]as well as { IEI n − , IEI n − , IEI n − : n = 1 , , , . . . } [similarly to that of Fig. 7(b)]. We then usethe Kolmogorov-Smirnov (KS) test [50] to checkwhether the distributions of corresponding elementsof these two sets { IEI n − , IEI n : n = 1 , , , . . . } and { IEI n , IEI n − , IEI n − , IEI n − : n =1 , , , . . . } are same or not. The p -values cor-responding to (cid:15) = 0 . . (cid:15) = 0 . p -values for six pairs( IEI n , IEI n − ), ( IEI n , IEI n − ), ( IEI n , IEI n − ),( IEI n − , IEI n − ), ( IEI n − , IEI n − ), and( IEI n − , IEI n − ) are 0 . . . . . . p -value denotesthe probability of obtaining test results at least asextreme as that of the observed one(s) under the nullhypothesis. Since the p -values are very high, we fail toreject the null hypothesis of equality of any of thesepairs of distributions. The KS test also confirms thatthe distributions of all the processes are also same asthat of the whole data of IEI. For different values of (cid:15) ,the processes can be thought as a renewal process (ofevents) [51] or more generally as a point process.We now try three known parametric family of distribu-tions (Weibull, Gamma, and Log-normal distributions)to fit the histogram of IEI (data collected from our nu-merical experiment) for two different aforementioned val-ues of (cid:15) in Figs. 7(c) and 7(d) to help us understand thenature of occurrence of next extreme events. We havealso observed that Log-normal distribution is the best fitfor both the cases. In Appendix B, we write the expres-sions of PDFs of IEI and present a chart of the estimatedparameters of these distributions (Table I). VI. PREDICTABILITY OF EXTREME EVENTS
The dependence structures of IEI are observed in Figs.7(a) and 7(b) which indicate about predictability of ex-treme events [52]. For predictability, we use ARIMAmodel [36], histogram of the IEIs and a correspondingbox plot analysis [37]. A preliminary description re-lated to the predictive property of box plot analysis isdescribed in Appendix C.ARIMA time-series forecasting model is used to fit thedata at first. For this, we find that the best fitted modelfor data is ARIMA (1 , , i.e. , a model with only au-toregressive part of lag one and no integrated or mov-ing average part. We use half of the 400 samples tofind a best fit of the model and on the remaining halfwe use the one-step-ahead prediction using the aboveAR(1), or, equivalently, the ARIMA (1 , ,
0) model. Weplot the iteration of our simulated IEIs (blue solid curve)and predicted IEIs (red dash curve) for (cid:15) = 0 . (cid:15) = 0 . , , for Fig. 8(a) and 10 for Fig. 8(b). This is reasonablefor predicting such extreme events, whose distribution oftime intervals of occurrence, IEI, have a very wide range.Here we observe that if the value of (cid:15) is shifted away fromcritical crisis point, the predictability of extreme eventsbecomes finer because the chaotic fluctuation of time sig- nal of T decreases varying with parameter and it reachesto MMO, where fluctuation of temporal evolution of sys-tem drops significantly. As a result, it may quite easyto predict such El Ni˜no events using ARIMA method.We investigate further the question of predictability in adifferent way. n I E I n I E I ( b )( a ) FIG. 8: Prediction plots of one-step-ahead forecast usingARIMA (1 , ,
0) time-series model based on the simulatedIEI’s for (a) (cid:15) = 0 . (cid:15) = 0 . T before appearing the eventheight 27 . (cid:15) = 0 . . A careful observation of the simulated time evolutionof T in Fig. 1(e) and the inset in Fig. 6(a) reveals thattwo types of spiking patterns appear, which may be char-acterized as extreme events. A class of large spikingevents called here as extreme events emerge after spirallyoutward with small amplitude oscillations; another classof extreme events emerge after a relatively longer dura-tion of bounded chaotic oscillations around the saddlefocus. To distinguish these two types of extreme events,we plot a set of last 10 local maxima of T before oc-curring an event of height larger than or equal to 27 . IEI N o . o f o cc u rr en c e s IEI N o . o f o cc u rr en c e s IEI N o . o f o cc u rr en c e s IEI N o . o f o cc u rr en c e s I E I I E I I E I I E I ( c ) ( d )( b )( a ) FIG. 10: Histogram of interevent intervals at (cid:15) = 0 . . .
5. The inset figures show the corresponding box plots. as shown in Fig. 9(a). Extreme events that emerge viatwo different processes, small amplitude oscillations andbounded chaotic oscillation, are clearly distinguishable.The peak value of the spiraling oscillation ( T max ) mono-tonically increases (black line) until they reach an eventheight larger than T max ≥ .
5. However, during theemergence of large events from the bounded chaotic os-cillation, there is no such increasing trend before the ex-treme events, rather a serpentine pattern (magenta line)is seen. Inspired from this fact, we separately investigatethe question of predictability of extreme events arisingout of the two cases.Figure 9(b) presents a histogram of IEI for all theevents of height 27 . lower adjacent = 12604 . firstquartile Q = 13291 . median = 22822 . third quartile Q = 46548 .
3, and upper adjacent = 96434 .
04. These measures tell usthat once an event of height 27 . . , . . , . . , .
3] with respective probabilities 0 .
25, 0 . .
75. It is almost sure that these particular type ofevents will occur within the interval [12604 . , . Q − Q = 33257 .
16, and approxi- mately 6.55% of the events are not predictable with thisstandard methodology, which in statistical terms may of-ten be called outliers .Next we investigate all extreme events of height 27 . . .
99 (right box) which are signifi-cantly small compared to bounded chaos. For both thesecases, we do not get any outliers which indicate betterpredictability of such events and the events occur throughspiral outward oscillation is almost predictable. Figure10(b) shows an long-tail unimodal distribution for theinterevent intervals of the extreme events of height 27 . T max exceeds the pre-defined threshold 25 .
5) instead of a specific height. Forsuch cases, similarly we have separated out all the eventswhich are occurring through spiral outward and boundedchaotic oscillations. Figures 10(c) and 10(d) respectively3 T m a x Time T
20 25 30
Event height -5 -4 -3 -2 -1 P D F ( d ) ( e ) ( f ) FIG. 11: Bifurcation diagram of the variable T with respect to (cid:15) for (a) r = and (d) r = . Temporal evolution of T at(b) (cid:15) = 0 . r = , and (e) (cid:15) = 0 . r = and corresponding probability density functions are shown in (c) and (f),respectively. Red lines in (e) and (f) are the threshold values at u ∗ = 25 .
65 (measured using mean excess function). represent the distributions of IEI for these two cases. Forthe spiral outward induce extreme event which exceedsthe length 25 .
5, almost every interevent intervals are con-centrated within a very small abscissa length. It showsthe feature of periodicity for the occurrence of extremeevent followed by small amplitude oscillations. How-ever, bounded chaos induces events that exhibit a uni-modal long-tail distribution. The IQR and outlier forsuch events are respectively 1971 .
79 and 5 . .
18 and 5 .
64% for bounded chaos case. This studyclearly indicates that the spiral outward induces eventsare more predictable than bounded chaos induces events.
VII. EFFECT OF SLOW-FAST TIMESCALE
The parameter r controls the timescale of the slow vari-able of the system (2). It plays a crucial role in the on-set of a sudden large expansion of the attractor from abounded state and in the origin of extreme events. Forelaboration, we fix all the parameters as considered aboveand draw bifurcation diagrams of T for two additionalchoices of r . Figures 11(a) and 11(d) show two bifurca-tion diagrams against (cid:15) for r = and r = , onelarger and one smaller than r = . We have alreadydrawn the scenario of PD and PA cascade collision in Fig.2(d) in Sec. IV for r = . The extreme event occurs at a postcrisis region after transition point (cid:15) ≈ . r = , an expansionof T max occurs at a transition point (cid:15) ≈ . r = value, then the critical (cid:15) for the onset of asudden large expansion shifts to a larger value, which isclearly observed in Fig. 11(d) [cf. Fig. 2(d)]. The corre-sponding time series at (cid:15) = 0 . u ∗ = 25 .
65 (horizontal red line),which is calculated by the mean excess function plot.PDF also changes significantly at (cid:15) = 0 . r = as show in Fig. 11(f) clearly shows a multimodalnon-Gaussian distribution. So it is clear that timescale r plays an important role for generating extreme events byshifting (cid:15) for the onset of extreme events. For a smallervalue of r , i.e. , for a slower rate of change in h , a largervalue of the strength of zonal advection (cid:15) is necessary toreach the crisis point where the enlargement of the at-tractor in phase space occurs leading to extreme events.4 VIII. CONCLUSION
A low-dimensional slow-fast ENSO model has been in-vestigated to understand the dynamical origin of extremeand related El Ni˜no events. Using bifurcation diagramsof SST against a system parameter (cid:15) connected to thestrength of difference in SST in the southern and east-ern Pacific, we identified two different kinds of interest-ing events, the typical chaotic MMO type El Ni˜no eventsand another extreme events. El Ni˜no events emerge morefrequently and almost in identical heights while the sec-ond type of events shows a wide variability in amplitudeand interevent intervals. For a demonstration of boththe events, besides drawing bifurcation diagrams of thesystem by varying the system parameter, we analyticallyderived a critical manifold of the system and located theregions of instabilities in phase space of the system. Wefound a critical value of the control parameter (cid:15) , where asudden expansion of the chaotic attractor was observeddue to an interior crisis. This interior crisis emerged dueto a merging of PD and PA cascade of bifurcations at acritical point when the extreme events emerge by a self-induced switching between small-amplitude oscillationsand large-amplitude events. In fact, a homoclinic chaoswas expected, in such a situation, with a local instabil-ity of the saddle focus of the slow-fast system along witha global stability of the trajectory when large identicalspiking events were usually seen alternately with ran-domly varying number of small oscillations. However, anadditional instability exists in phase space of the systemin the form of a channel-like structure that prevents theglobal stability of the trajectory leading to large varia-tion in amplitude and return time of the large spikingevents. This variation in amplitude depends on the vol-ume of this channel. For a relatively narrow channel,a trajectory while spirally moving out of a saddle fo-cus of the system, spent a longer period of time insidethe channel and finally made a large global excursionto form an extreme event. The reinjected trajectoriesor the repeat global excursions were no more stable andirregular leading to a large variation of amplitude andreturn time of events. For a wider channel, a trajectoryeasily passed through it without spending much time in-side it and hence the reinjected trajectories in the repeatglobal excursion were more stable to make almost iden-tical height and more frequent chaotic MMO-like typicalEl Ni˜no events as reported earlier. We classify the largespiking events as extreme events, which are larger thanan estimated a threshold defined by a mean-excess func-tion.PDF of event heights showed a multimodal distributionwhich tends to be bimodal with an increase in the sys-tem parameter (cid:15) . For characterization of the occurrenceof large events, PDF of interevent intervals has been fit-ted by Weibull, Gamma, and Log-normal distributionsand out of the three, the Log-normal distribution is bestfitted. A dependence nature of the interevent intervals isobtained from ACF plot that illustrates possible impli- cation of the predictability of (resulting) extreme events.We have discussed the predictability of extreme events,first using ARIMA models, and then in more details withbox plot analysis. We show that the ARIMA models canpredict the events with a standard error of the order of10 . We then indulge in in-depth analysis via box plot.For our predictability purpose, we categorize each eventsin two different classes based on their emergence. Oneclass of events takes much shorter time to occur, whichoriginate via spiraling out from the saddle set and from aknowledge of such events and their occurrence, their re-turn time is predicted with a good accuracy; the secondclass of events on an average takes much longer time andthey originate via bounded chaos where they spend moretime before a large event. Our box plot analysis reveals,knowing that such events occur, predictability is foundless accurate than the former one, because of the dura-tion of the time spent in bounded chaotic motion. Theslow-fast ratio parameter of the system plays a signifi-cant role in the onset of extreme events. While real dataare available on chaotic MMO-like decadal emergence ofextreme events, we are yet to find a real data set to verifythe existence of our classified extreme events. However,we explain how and when two events, El Ni˜no events andextreme events may emerge in the low-dimensional slow-fast climate model against a system parameter variation,and explained the phenomena using a common dynamicalmechanism, which has so far been missing. Furthermore,we tried to address the question of the predictability ofextreme events. Acknowledgments:
DG was supported by SERB-DST (Department of Science and Technology), Govern-ment of India (Project No. EMR/2016/001039). SKDacknowledges support by the Division of Mechanics,Lodz University of Technology, Poland. AR thanksArindam Mishra for useful discussions.
APPENDIX A
We make a linear stability analysis of the equilib-rium points of the system (7), which has one axialequilibrium point (0 , T r , T r ) and one interior equilibriumpoint ( h ∗ , T ∗ , T ∗ ). The stability of the axial equilibrium(0 , T r , T r ) is obtained from the Jacobian matrix J of the3D system (7), J = − r bLrµ β − bLrµ β − α − p − α + p , where p = ζµ T r − T r0 (cid:8) − tanh H − z h ∗ (cid:9) . Its eigenvalues are − r, − α, p − α . For our choice of system parameters, thevalue of p becomes 0 . p > α and hence5the axial equilibrium is a saddle point. Now we obtain,[ J + rI ] m = p − α ) m − bLrµ β − ( p − α ) m − bLrµ β − α m − α m − ( p − α ) m ( p − α ) m , when the generalized eigenvector of order m correspond-ing to the eigenvalue − r is [1 0 0] tr (here tr denotes thetranspose of the matrix), for all m ∈ N . For the othereigenvalue − α , one can write[ J + αI ] m = ( α − r ) m ( α − r ) m − p m α − r − p bLrµ β ( α − r ) m − p m p + r − α bLrµ β − p m p m , when its eigenvector is [0 1 1] tr .The stable subspace of the saddle point (0 , T r , T r ) isthe eigenspace generated by the generalized eigenvector { [1 0 0] tr , [0 1 1] tr } , which, in reality, is the T = T plane. The T = T plane is the invariant manifold of thesystem since when we have ddt [ T − T ] = 0 and, for all thechoices of initial conditions, the system remains confinedto this plane. The tangent space of the T = T plane is T = T itself, which is a stable subspace of system (7).In this stable manifold, our original 3D system can bereduced to a 2D system,˙ h = − rh , ˙ T = − α ( T − T r ) , (14)where T = T = T . It is easy to check that the axialequilibrium point (0 , T r , T r ) is a stable node in this plane.Since our system is 3D and the stable manifold is 2D, sothe unstable manifold is 1D, which is spanned by thegeneralized eigenvector [ α − r − p ) βbLrµ p − α of J .For T (cid:54) = T , the interior equilibrium point ( h ∗ , T ∗ , T ∗ )satisfies h ∗ = − bLµ ( T ∗ − T ∗ )2 β , α ( T ∗ − T r ) + (cid:15)µ ( T ∗ − T ∗ ) = 0 , αµ + ζT r ) − ( (cid:15) + ζ ) T ∗ + (cid:15)T ∗ ] ζ ( T r − T r ) = 1 − tanh (cid:16) H + bLµ ( T ∗ − T ∗ )2 β − z (cid:17) h ∗ . (15)Solving Eq. (15), we derive the interior equilibriumpoint with one real and two complex conjugate eigen-values. A bifurcation diagram of this interior equilibriumpoint is plotted by varying (cid:15) in the range [0 ,
1] as shown inFig. 12(a). The interior equilibrium point is stable until (cid:15) ≈ . (cid:15) ≥ . (cid:15) ≈ . (cid:15) is verified by plotting the real parts of threeeigenvalues of the Jacobian matrix of the system at inte-rior equilibrium point is shown in Fig. 12(b). It is noticedthat, for 0 ≤ (cid:15) ≤ . λ i =1 , , remain negative which signify stability ofthe interior equilibrium point in this parameter interval.Above (cid:15) = 0 . (cid:15) = 0 . (cid:15) , the system becomes oscil-latory either in periodic [green dots] or in chaotic [bluedots] state [cf. Fig. 12(a)]. At (cid:15) = 0 . (cid:15) ∈ (0 . , . (cid:15) values, which is our main focus of study anddescribed in the Sec. IV, in detail.The stable manifold, T = T plane (gray color)is drawn in Fig. 13(a) on which the saddle point(0 , . , .
5) (black circle) lies. For all initial conditionsfrom the left side of this plane, any trajectory of the sys-tem goes unbounded along the unstable manifold of thesaddle point (0 , . , . h ∗ , T ∗ , T ∗ ), which lies not faraway from the T = T plane. Some exemplary trajec-tories are shown in different colors for various choice ofinitial conditions. The trajectories move toward the at-tractor along the stable manifold (eigendirection corre-sponding to the negative real eigenvalue) of the saddlefocus, but spirally moves away as repelled by the un-stable manifold (corresponding to the complex conjugateeigenvalues with positive real parts) of the saddle focus.This is an ideal situation for the origin of homoclinicchaos, when the trajectory travels a close vicinity of thesaddle focus along the stable eigendirection, but stronglyrepelled and spiral out along the unstable plane of thesaddle focus and makes large excursions on the samepath repeatedly and return to another close vicinity of6 T R e ( i ) ( a )( b ) FIG. 12: (a) Bifurcation diagram of the interior equilibriumpoint ( h ∗ , T ∗ , T ∗ ) with respect to (cid:15) by plotting the maxima of T . (b) Variation of the real parts of the eigenvalues for theinterior equilibrium point by changing (cid:15) . PD, period doubling;HB, Hopf bifurcation; IPD, inverse period doubling; LPC,limit point bifurcation of cycles. the saddle focus. Instead, the trajectory of ENSO sys-tem is locally unstable in the close vicinity of the saddlefocus as usual, but becomes globally unstable too. Asa result, the global trajectory is not stable, but tracesdifferent paths as shown in Fig. 13 (a) that creates vari-ation in the amplitude of large spiking events and theinterspike intervals, which we explain in the Sec. IV. Theright side of the plane T = T is the basin of attrac-tion of those events. Fig. 13(b) is a 2D projection of thevector field, which clearly enunciates that all the initialconditions from the lower-half of the T = T line havepossibility to originate extreme events. For other initialconditions from the upper-half, i.e., T > T , the systemgoes unbound along the unstable manifold of the saddlepoint. In reality, T > T is also impossible situation.In this way, the equilibrium point (0 , T r , T r ) also playsa significant role to determine the basin of attraction ofthe extreme events. APPENDIX B
Here, we have described three PDF. PDF of theWeibull distribution is P ( r ) = kθ (cid:16) rθ (cid:17) k − e − ( rθ ) k , r ∈ [0 , ∞ ) . (16)PDF of the Gamma distribution is P ( r ) = 1Γ( k ) θ k r k − e − rθ , r ∈ (0 , ∞ ) . (17) T T ( b ) FIG. 13: (a) Trajectories of the chaotic attractor (blue line)in 3D plane. Trajectories for different initial conditions arerepresented by different colors. The diagonal T = T plane(brown) denotes a stable manifold of the axial equilibriumpoint (0 , . , .
5) (black circle). (b) Two-dimensional vec-tor field at h = 0 plane where black line (in Z shape) rep-resents the projection of the critical manifold on this plane.Parameter value: (cid:15) = 0 . PDF of the Log-normal distribution is P ( r ) = 1 √ πrσ exp ( − (log r − µ ) σ ) , r ∈ (0 , ∞ ) . (18)Here, k ( >
0) and θ ( >
0) are shape and scale parametersfor the first two distributions [16 and 17], µ and σ ( > r represents interevent interval. We es-timate the parameters which we have recorded are givenin the Table I:Table I (cid:15) = 0 . (cid:15) = 0 . θ =23463.1 θ =15815.20 k =1.49827 k =3.5059Gamma distributon θ =7322.98 θ =706.46 k =2.8527 k =20.4212Log-normaldistributon µ = 9 . µ = 9 . σ = 0 . σ = 0 . Appendix C: Box plot analysis
FIG. 14: Schematic diagram of a box plot: A blue box isdrawn from the first quartile (Q ) to the third quartile (Q ).Red horizontal straight line is traced at the median of the dataset. Two black dashed vertical lines are drawn respectivelyfrom lower adjacent to Q and Q to upper adjacent . Red (+)markers denote the outlier events, which are least predictable. A box plot [37] is a standard contrivance to displaythe distribution of a given data. It basically producesa five measures summary, namely lower adjacent , firstquartile Q , second quartile ( median ), third quartile Q ,and upper adjacent . A box plot gives more informationthan the measures of central tendency. For informationon variability or dispersion of IEI, it gives an impressionabout the spread of IEI values. Here median is the second quartile, which is the middle-most value of the dataset.The first quartile Q , is the middle point between the median and the smallest number of the whole data set.While the third quartile Q is the middle point betweenthe median and the highest value of the whole data set.IQR is defined as the difference of Q and Q . The loweradjacent and upper adjacent are not always the smallestand largest values in the data set, respectively. Here, lower adjacent is the smallest value of the data set or(Q − . × IQR), whichever is the maximum. While upper adjacent is the minimum of the maximum value ofthe data set and (Q +1 . × IQR). The events from loweradjacent to upper adjacent are predictable events, and therange of upper adjacent to lower adjacent is known asthe range of predictable events. Although, we use largesize of data-set as available from our simulations and asnecessary for such kind of box-analysis.A schematic diagram of the box plot is shown in Fig.14. The five measures approximately divide the entiredata set into four sections, each one approximately con-tains 25% of the data sets. Beyond the upper adjacent and lower adjacent values, events are outliers, which areseemingly unpredictable and much occurred rarely. Theoutliers (red plus) tell us which events are out of our pre-dictable range. The probability of occurrence of an eventwithin the interval [ lower adjacent , Q ] is 0 .
25. Whilethe probability of an event that lies within the intervals[ lower adjacent , median] and [ lower adjacent , Q ] are 0 . .
75, respectively. It is almost sure to lie the pre-dictable events within the interval [ lower adjacent , upperadjacent ]. Five measures of box plot corresponding tothe Fig. 9(b) and Figs. 10(a) and 10(d) are given in theTable II:Table II Box plot Lower Adjacent Q Median Q Upper Adjacent OutliersFig. 9(b) 12604.38 13291.14 22822.3 46548.3 96434.04 6.55%Fig. 10(a) 12604.38 12615.708 12626.15 12641.57 12659.5 0.013267.38 13276.74 13284.7 13299.73 13322.18 0.0Fig. 10(b) 14980.68 23718.93 36832.34 62234.63 120008.18 5.05%Fig. 10(c) 8345.9 11303.59 12609.16 13275.38 13322.18 5.74%Fig. 10(d) 8955.99 16956.65 22821.38 37095.83 67282.32 5.64% [1] H. A. Dijkstra,
Nonlinear Physical Oceanography: A Dy-namical Systems Approach to the Large Scale Ocean Cir-culation and El Ni˜no (Springer Science, New York, 2005).[2] E. S. Sarachik, and M. A. Cane,
The El Ni˜no - SouthernOscillation Phenomenon (Cambridge University Press,Cambridge, U.K., 2010).[3] S. -I. An, and F. -F. Jin, J. Clim. , 2399 (2004).[4] A. Timmermann, and F.-F. Jin, Geophys. Res. Lett. ,3-1 (2002).[5] A. S. Sharma, A. Bunde, V. P. Dimri, and D. N. Baker, Extreme Events and Natural Hazards: The Complexity Perspective (American Geophysical Union, Washington,D.C., 2013).[6] O. Rojas,Y. Li, and R. Cumani,
Understanding theDrought Impact of El Ni˜no on the Global AgriculturalAreas (FAO, Rome, Italy, 2015).[7] A. Timmermann, J. Oberhuber, A. Bacher, M. Esch, M.Latif, and E. Roeckner, Nature , 694 (1999).[8] F. -F. Jin, J. Atmos. Sci. , 811-829 (1997); F. -F. Jin,J. Atmos. Sci. , 830 (1997).[9] A. Timmermann, F. F. Jin, and J. Abshagen, J. Atmos.Sci. , 152 (2003). [10] S. E. Zebiak, and M. A. Cane, Mon. Wea. Rev. , 2262(1987).[11] J. Bjerknes, Mon. Wea. Rev. , 163 (1969).[12] A. Roberts, J. Guckenheimer, E. Widiasih, A. Timmer-mann, and C. K. R. T. Jones, J. Atmos. Sci. , 1755(2016).[13] J. Guckenheimer, A. Timmermann, H. Dijk-stra, and A. Roberts, Dyn. Stat. Clim. Sys., ,https://doi.org/10.1093/climsys/dzx004 (2017).[14] M. Ghil et al. Nonlin. Process. Geophys. , 295 (2011).[15] M. Latif, D. Anderson, T. Barnett, M. Cane, R. Kleeman,A. Leetmaa, J. O’Brien, A. Rosati, and E. Schneider, J.Geophys. Res. , 375 (1998).[16] G. G. Nobre, S. Muis, T. IE Veldkamp, and P. J. Ward,Progress in Disaster Science , 100022 (2019).[17] M. Ghil, and N. Jiang, Geophys. Res. Lett. , 171(1998).[18] M. Ghil, M. R. Allen, M. D. Dettinger, K. Ide, D. Kon-drashov, M. E. Mann, A. W. Robertson, A. Saunders, Y.Tian, F. Varadi, and P. Yiou, Rev. Geophys. , 1003(2002).[19] A. N. Pisarchik, R. Jaimes-Re´ategui, R. Sevilla-Escoboza, G. Huerta-Cuellar, and M. Taki, Phys. Rev.Lett. , 274101 (2011).[20] C. Bonatto, M. Feyereisen, S. Barland, M. Giudici, C.Masoller, Jose R. Rios Leite, and J. R. Tredicce, Phys.Rev. Lett. , 053901 (2011).[21] Hugo L. D. de Souza Cavalcante, M. Ori´a, D. Sornette,E. Ott, and D. J. Gauthier, Phys. Rev. Lett. , 198701(2013).[22] V. Lucarini et al. Extremes and Recurrence in DynamicalSystems (John Wiley & Sons, New York, 2016).[23] Y. S. Fan, and T. R. Chay, Phys. Rev. E , 1012 (1995).[24] R. Karnatak, G. Ansmann, U. Feudel, and K. Lehnertz,Phys. Rev. E , 022917 (2014).[25] C. Grebogi, E. Ott, and J. A. Yorke, Physica D , 181(1983).[26] C. Grebogi, E. Ott, F. Romeiras, and J. A. Yorke, Phys.Rev. A , 5365 (1987).[27] A. Ray, S. Rakshit, D. Ghosh, and S. K. Dana, Chaos , 043131 (2019).[28] S. Leo Kingston, K. Thamilmaran, P. Pal, U. Feudel, andS. K. Dana, Phys. Rev. E , 052204 (2017).[29] N. Akhmediev et al. J. Opts. , 063001 (2016).[30] J. A. Reinoso, J. Zamora-Munt, and C. Masoller, Phys.Rev. E , 062913 (2013).[31] J. Zamora-Munt, B. Garbin, S. Barland, M. Giudici, J.R. R. Leite, C. Masoller, and J. R. Tredicce, Phys. Rev. A , 035802 (2013).[32] S. Coles, An Introduction to Statistical Modeling of Ex-treme Values (Springer, Berlin, 2001).[33] G. F. de Oliveira, Jr, O. Di Lorenzo, T. P. de Silans, M.Chevrollier, M. Ori´a, and Hugo L. D. de Souza Caval-cante, Phys. Rev. E , 062209 (2016)[34] A. Mishra, S. Saha, M. Vigneshwaran, P. Pal, T. Kapi-tanaiak, and S. K. Dana, Phys. Rev. E , 062311 (2018).[35] A. Ray, A. Mishra, D. Ghosh, T. Kapitaniak, S. K. Dana,and C. Hens, Phys.Rev. E , 032209 (2020).[36] C. Chatfield, Time-series Forecasting (CRC Press, BocaRaton, FL, 2000).[37] J. W. Tukey,
Exploratory Data Analysis (Addison-Wesley, London, 1977); F. M. Dekking, C. Kraaikamp,H. P. Lopuha¨a, L. E. Meester,
A Modern Introduction toProbability and Statistics: Understanding Why and How (Springer, London, 2005).[38] J. Guckenheimer, SIAM J. Appl. Dyn. Syst. (4), 1355(2008).[39] S. K. Dana, S. Chakraborty, and G. Ananthakrishna,Pramana, J. Phys. , 443 (2005); S. Chakraborty, andS. K. Dana, Chaos , 023107 (2010).[40] N. Fenichel, J. Diff. Eqns. , 53 (1979).[41] M. Desroches, J. Guckenheimer, B. Krauskopf, C.Kuehn, H. M. Osinga, and M. Wechselberger, SIAM Rev. (2), 211 (2012).[42] G. Ansmann, R. Karnatak, K. Lehnertz, and U. Feudel,Phys. Rev. E , 052911 (2013).[43] G. Ansmann, K. Lehnertz, and U. Feudel, Phys. Rev. X , 011030 (2016).[44] B. Das, and S. Ghosh, Extremes , 325 (2016).[45] S. Nag Chowdhury, S. Majhi, M. Ozer, D. Ghosh, andM. Perc, New J. Phys. , 073048 (2019).[46] S. R. Massel, Ocean Surface Waves: Their Physics andPrediction (World Scientific, Singapore, 1996).[47] T. Sapsis, Phil. Trans. R. Soc. A , 20170133 (2018).[48] J. D. Hamilton,
Time Series Analysis (Princeton Univer-sity Press, Princeton, NJ, 1994).[49] P. J. Brockwell, and R. A. Davis,
Time Series: Theoryand Methods (Springer, Berlin, 1991).[50] J. W. Pratt, and J. D. Gibbons,
Concepts of Nonpara-metric Theory (Springer, New York, 1981).[51] K. V. Mitov, and E. Omey,
Renewal Process (Springer,Berlin, 2014).[52] S. Albeverio, V. Jentsch, and H. Kantz (eds.),