Research articlesThe effect of state dependence in a delay differential equation model for the El Niño Southern Oscillation
Abstract
Delay differential equations (DDEs) have been used successfully in the past to model climate systems at a conceptual level. An important aspect of these models is the existence of feedback loops that feature a delay time, usually associated with the time required to transport energy through the atmosphere and/or oceans across the globe. So far, such delays are generally assumed to be constant. Recent studies have demonstrated that even simple DDEs with non-constant delay times, which change depending on the state of the system, can produce surprisingly rich dynamical behaviour. Here, we present arguments for the state dependence of the delay in a DDE model for the El Niño Southern Oscillation phenomenon in the climate system. We then conduct a bifurcation analysis by means of continuation software to investigate the effect of state dependence in the delay on the observed dynamics of the system. More specifically, we show that the underlying delay-induced structure of resonance regions may change considerably in the presence of state dependence.
This article is part of the theme issue ‘Nonlinear dynamics of delay systems’.
1. Introduction
Feedback loops are an important component of climate systems (see [1] and references therein), whereby changes or inputs are either amplified or dampened by nonlinear processes. Such feedback loops are crucial ingredients in climate models [2–7] and they are associated with delay times due to transportation times of mass or energy over large distances. These time delays are generally large compared with the characteristic forcing time scales of the system under consideration. As such, time delays tend to play a vital role in the behaviour of the system.
This is certainly true for delayed feedback loops that arise in the equatorial coupled ocean-atmosphere system governing the El Niño Southern Oscillation (ENSO) phenomenon. This system consists of an atmospheric component, known as the Southern Oscillation, and an oceanic component, El Niño. The El Niño events are characterized by large increases in sea surface temperature (SST) in the equatorial Pacific Ocean, and they occur sporadically approximately every 4–7 years. Large drops in eastern Pacific SST represent the cool phase and are known as La Niña events. It is well known that El Niño events generally tend to occur near Christmas time. The seasonal cycle (with a period of 1 year) represents the characteristic forcing time scale of the ENSO system. Feedbacks due to ocean-atmosphere coupling processes in the eastern and central equatorial Pacific Ocean have delay times of many months due to the time needed for equatorially trapped waves to propagate across the Pacific Ocean.
El Niño and La Niña events have major consequences all over the globe, yet they are notoriously difficult to predict accurately, even with highly sophisticated global climate models [8,9]. Over many decades, conceptual climate models have proven extremely useful in providing fundamental insights into the inner workings of climate systems; for example, see [10–14]. Conceptual ENSO models have played an important role in building a theory of how El Niño events are created and how they eventually die away, although there are many aspects still open to debate and investigation [12].
A prominent paradigm for how El Niño events form and abate is the so-called delayed action oscillator (DAO), first introduced by Suarez & Schopf [15]. According to this paradigm, feedback loops due to the ocean-atmosphere coupling are crucial in shaping ENSO dynamics. A central ingredient of these feedback loops is the depth of the thermocline, which is a thin layer in the ocean separating deeper cold waters from shallower warm water. Another ingredient is the equatorial upwelling which is caused by the surface winds (trade winds) through a mass divergence away from the equator caused by the Earth's rotation (Ekman upwelling). Let the variable h(t) denote the deviation from a long-term thermocline mean depth in the eastern equatorial Pacific. A positive value of h leads to less water from below the thermocline reaching the sea surface through upwelling, thus an increased SST in the eastern Pacific. This provides a positive feedback because the warm SST anomaly slows down the easterly trade winds since warm air rises and has to be compensated, leading to westerly wind anomalies, further deepening the thermocline. The positive anomalous thermocline depth signal at the equator is carried to the eastern boundary of the ocean via so-called Kelvin waves. Off the equator, a corresponding negative anomalous thermocline depth signal is carried to the western boundary of the ocean via so-called Rossby waves. These waves are reflected as Kelvin waves, which carry the shallow thermocline perturbation back to the eastern boundary of the ocean, completing a negative feedback loop.
There exist a number of models based on the DAO paradigm; for example, see [15–18]. The DDE model introduced by Ghil et al. [19] is one of the simplest in that it focuses on the interaction between the negative delayed feedback and additive seasonal forcing. Neglecting the above-mentioned positive feedback is not necessarily justified in [19], rather it is done in order to demonstrate that the other two terms are sufficient to produce very rich dynamical behaviour relevant to ENSO. The effect of re-introducing the positive delay feedback to this model is studied in detail in [20]. The model from [19] takes the form
Modelling with DDEs is built on the idea that only the delayed effect of a chain of complex processes needs to be modelled, and not the processes themselves—resulting in a convenient model reduction. The resulting DDE model is capable of complicated behaviour, yet is still amenable to sophisticated mathematical analysis. In the case of model (1.1), only the delayed effect of the negative feedback is included, while a full description (for example, by partial differential equations) would involve many more variables and parameters relating to zonal wind anomaly, surface pressure anomaly, thermal damping coefficients, etc. The GZT model is deceivingly simple in that it features only delayed feedback and periodic forcing. Yet, it was shown in [19,23] to be capable of dynamics that could replicate certain ENSO features, such as variability on appropriate time scales and phase locking with the seasons. These are features that have been observed earlier in more complex ENSO models [24–26]. The balance between the complicated dynamics produced by DDE models and their (relative) simplicity, allowing them to still be mathematically tractable, explains their success in areas as diverse as, for example, laser dynamics [27], neural networks [28], traffic flow [29], cell population dynamics [30] and epidemiology [31].
DDEs are infinite-dimensional dynamical systems and their analysis generally poses analytical and numerical challenges. Nonetheless, the theory of DDEs with a finite number of constant delays is well developed [32,33]. In fact, the bifurcation theory of constant-delay DDEs is analogue to that of ordinary differential equations: their solutions depend smoothly on parameters and initial conditions, and their linearizations around equilibria and periodic solutions have at most finitely many unstable eigendirections. Furthermore, advanced numerical tools for their simulation and bifurcation analysis are available and have been applied successfully to analyse models arising in various applications; for example, see [34] and references therein.
In [20,35], we conducted a bifurcation analysis of the GZT model (1.1) in order to understand how it is able to produce such rich dynamics. This was done by means of the MATLAB-based continuation software DDE-Biftool [36,37], which allows the user to track (or continue) steady-state and periodic solutions of DDEs while varying a chosen model parameter. It can calculate the stability of such solutions, allowing one to detect bifurcations. It is then possible to continue detected bifurcations of steady-state and periodic solutions in two-parameter planes. The bifurcation analysis in [35] revealed that in the (c, τn)-plane the dynamics of the GZT model can be divided roughly into three parameter regions where the dynamics is: (i) periodic due to the negative delayed feedback alone, (ii) periodic and dominated by the seasonal forcing and (iii) an interplay of both mechanisms. In the latter case, one finds dynamics on an invariant torus, which may be locked (periodic) or unlocked (quasi-periodic). The periodic solutions on the torus are organized in the parameter plane into p:q resonance tongues, where they wind around the torus to form p:q torus knots. It was shown in [35] that there exist small regions in the (c, τn)-plane where multiple resonance tongues overlap. In these regions, period-doubling cascades and chaotic solutions exist.
Of course, the use of a constant value for any delay in a DDE model, such as the GZT model, is a modelling assumption. A recent review of the use of DDEs in conceptual climate models [38] confirms that the delays in DDE climate models have generally been assumed to be constant. While this assumption is well justified in certain applications, such as machining [39] and laser dynamics [40], delays in many applications are definitely not constant. Often the delay times will depend on the state of the system itself, which leads to DDE models with state-dependent delays. The theory of and numerical methods for state-dependent DDEs is more difficult and still under development [37,41–43].
Recently, Calleja et al. [44] studied the dynamics of a simple scalar DDE where two delay times depend linearly on the state of the system. The associated constant-delay DDE is linear and capable of only trivial dynamics. However, with state dependence of the delays, the system generates surprisingly complicated dynamics, as was demonstrated by a combination of numerical continuation and normal form calculations. In light of this study, we must ask whether, after nearly 40 years of DDE climate models, is it still a safe modelling assumption to take their delay times as constant?
This paper contributes by studying the effect of state dependence on the dynamics of model (1.1). We introduce these effects in a heuristic fashion in the hope that our results will provide motivation for a more thorough and quantitative derivation of state dependence in ENSO. This would involve either the analysis of more complex model equations or data from such models or observations. More generally, this paper can be considered as a first step in addressing two fundamental questions:
| (i) | When do state-dependent delays arise in climate variability phenomena and what mathematical forms do they take? | ||||
| (ii) | Does the state dependence in a delay have a significant effect on the behaviour of the solutions of delay-equation models of these phenomena? | ||||
One reason to choose model (1.1) as a test-case example is that its bifurcation analysis in the (c, τn)-plane, presented in [35], provides a good starting point for observing the influence of state dependence on model behaviour. In spite of its simplicity among DAO models, the constant-delay GZT model already features complicated resonance phenomena. Moreover, model (1.1) can be viewed as prototypical in that it describes the interplay between delayed negative feedback and periodic forcing. Therefore, the results presented here may provide insight into whether or not state dependence of delays is important also in the context of other climate models, as well as applications from other fields. For example, DDE models of similar form can be found in human motion control [45], network dynamics [46] and laser dynamics [47]. Recently, the importance of incorporating state dependence in modelling has been recognized in turning processes [48], physiology [49] and lasers with frequency-dependent feedback [50].
We begin with considering in some detail how non-constant delay times in the negative feedback loop arise from the physics of the ocean-atmosphere coupling. Here we focus on how the surface of the ocean couples with the thermocline below. This leads to two different delay terms with state dependence. The first, denoted τc, is due to the coupling process in the central equatorial Pacific Ocean, and it depends linearly on the current thermocline depth h(t). The second, denoted τe, is due to coupling processes in the East and depends linearly on the state h in the past. Overall, we obtain a state-dependent version of model (1.1) with two extra parameters that allow us to ‘switch on’ these two types of state dependence; as part of our modelling, we also discuss realistic ranges for the respective state dependence. We present bifurcation sets in the -plane that demonstrate the effect of the two types of state dependence, first individually and then in combination; here is the constant part of the now state-dependent delay. As the starting point for our study, we first briefly review the bifurcation set in the -plane for the constant-delay case from [35]. It turns out that, over the range we consider, the state-dependent term τe alone has a negligible effect on the bifurcation set. The term τc, on the other hand, has a significant impact on the resonance regions in the -plane: we find reconnecting torus bifurcation curves, considerable increases in certain resonance regions and much more overlapping between resonance regions and associated period-doubling cascades. Interestingly, in the presence of state dependence in τc, the state-dependent term τe does have a considerable influence on the bifurcation set, increasing the observed complexity even further. This includes the breaking of a reflection symmetry, which in the past was only achieved by introducing an asymmetric coupling function [18,20,21]. Overall, we find for realistic levels of state dependence a bifurcation set that is much more complex than its constant-delay counterpart. In particular, a computation of Lyapunov exponents shows that there are considerably larger parameter regions with chaotic behaviour associated with overlapping resonance regions. Since this type of complicated and irregular behaviour has been identified as showing important characteristics of ENSO, we argue that state dependence in model (1.1) emerges as an important ingredient for a more realistic description of its underlying feedback mechanisms.
The paper is organized as follows. In §2, we discuss how the delay time in model (1.1) can be modified to include two types of state dependence that arise from surface-thermocline interaction. The resulting dynamics is studied in §3, where we review the constant-delay case in §§3a, consider each state-dependence term individually in §§3b, and then in combination in §§3c. Finally, in §4, we make some concluding remarks and point out future research directions.
2. State dependence of the delay in model (1.1)
We now motivate the inclusion of state-dependent terms into the delay time of model (1.1). There are various potential sources of state dependence in the ENSO system. The type of state dependence we present here arises from the interaction times between the thermocline and the sea surface, which are not described by and go beyond the original DAO mechanism.
The idea that both wave dynamics and thermocline-SST interactions are important for accurately describing ENSO behaviour is supported by [51–54]. More specifically, the state-dependent DDE model (1.1) presented in this section, with the state-dependent delay given by (2.6), is related to previous conceptual-level models for ENSO with underlying physics that can be traced back to the Zebiak-Cane model [55]. As summarized by Dijkstra [4, section 7.5.4], applying slightly different assumptions to equations derived from the intermediate-complexity Zebiak–Cane model will lead to the DAO, as well as the coupled wave oscillator [22] and the recharge oscillator [56,57]. Each of these ENSO paradigms describes the negative feedback and the propagation of oceanic waves, albeit in different mathematical forms. The coupled wave oscillator model [22] assumes that the effect of the thermocline depth on the SST in the East is instantaneous. The DAO and recharge oscillator do not make this assumption. However, neither attempt to resolve thermocline-SST interactions to include the effect of an evolving thermocline. As we will see in the following section, this is an effect that can have important consequences for the behaviour of the system.
(a) Components of the negative feedback
Figure 1 is a schematic diagram of ocean–atmosphere coupling processes in the central and eastern equatorial Pacific Ocean that are involved in the negative feedback loop of ENSO. The horizontal direction represents longitude and the vertical direction represents depth. The top curve represents the ocean surface and the lower curve the thermocline layer, the side curves indicate the basin boundaries of the Pacific Ocean, and the ellipse represents a coupling interaction zone. In its mean state the thermocline depth, which is indicated by the dashed line, is shallower in the East (about 50 m) compared with the depth in the central equatorial Pacific Ocean (about 150 m). The arrows represent the four components of the negative feedback: upwelling in the East, easterly winds, ocean adjustment in the central Pacific and wave propagation. Each component is associated with a delay time that appears in brackets in figure 1.
Figure 1. Schematic diagram of ocean-atmosphere interaction in the central and eastern equatorial Pacific Ocean following perturbations to the thermocline depth in the eastern equatorial Pacific. The dashed line represents the mean thermocline depth. (Online version in colour.)
The overall negative feedback loop is represented in figure 2 as a flow diagram. Its first component consists of a positive perturbation in the thermocline depth in the eastern equatorial Pacific that increases the SST in the eastern equatorial Pacific. This happens via an upwelling process with associated delay time τe. Next, the increase in SST slows down the easterly trade winds, creating a westerly wind anomaly over the central equatorial Pacific. This is a fast process and it is assumed here not to be subject to a delay (as is the case in all DAO models). The wind anomaly induces an ocean mass adjustment that decreases the thermocline depth in the central equatorial Pacific. This occurs in a reasonably localized interaction zone, and in mathematical derivations of DAO models [22,57], this zone is simplified to a point. The associated surface-thermocline interaction comes with a delay time τc. Rossby waves then carry the signal to the western basin boundary and are reflected as Kelvin waves, which carry the signal back to the East with an associated delay time τw. The total delay time associated with the negative feedback loop is therefore
Figure 2. Flow diagram of the negative delayed feedback due to ocean–atmosphere coupling; all quantities involved represent deviations from their respective means.
In past studies of DAO models, the delay time τn associated with the negative feedback has been assumed to consist only of τw, which is taken to be constant. We now take into account also the processes of upwelling and ocean adjustment and the associated delay times τe and τc, as sketched in figures 1 and 2. Importantly, while they are on average smaller than τw, we argue now that τe and τc are relevant and, in particular, that they depend on the state. The underlying principle is that the distance the signal must travel to or from the surface depends on the depth of the thermocline as represented by its deviation h from its long-term average.
(b) State dependence due to upwelling
There is a delay, indicated by τe in figure 2, between the thermocline anomaly at depth and the SST anomaly at the surface; this is due to upwelling of colder water. We model this effect by
(c) State dependence due to ocean adjustment
We now address the delay between the SST and thermocline depth anomalies in the central equatorial Pacific Ocean, represented by the third component of the feedback loop in figure 2. This can be attributed to a basin adjustment involving a collection of waves, which can be viewed as an adjustment to Sverdrup balance at off-equatorial latitudes. This adjustment is a crucial element in the so-called recharge-oscillator view of El Niño [56,57]. Due to the westerly wind-stress anomaly, which is arising through a positive SST anomaly, mass is diverging out of the equatorial zone because of this Sverdrup adjustment. Note that, although the associated arrow in figure 1 points downward to show the direction of the feedback loop, there is actually mass adjustment both in and out of the equatorial zone. Only after an adjustment time τc will this have an effect on the thermocline depth. Because the mass transport is dependent on the thermocline depth, so is τc. We model this effect here by
(d) Total delay time
From (2.1), the total state-dependent delay then takes the form
(e) Suitable parameter ranges
Based on the correlation analysis of observational SST and thermocline depth data by Zelle et al. [58], we estimate the constant delay times and to be two weeks and four months, respectively.
The maximum deviations in the thermocline depth in the eastern equatorial Pacific Ocean corresponds to about 50 m [59], which is also the mean thermocline depth, suggesting that weeks. Since in the GZT model, the maximum values of h(t) are of order 1 and time is in years, we obtain the nominal value of ηe ≈ 2/52 ≈ 0.04. Similarly, the maximum deviations in thermocline depth in the central equatorial Pacific Ocean correspond to about one third of its mean depth, suggesting that months. Therefore, we obtain the nominal value ηc ≈ (4/3)/12 ≈ 0.11.
We use these nominal values for ηe and ηc in the bifurcation analysis in the next section. In light of various approximations and uncertainties, and to present a more comprehensive picture of the possible dynamics, we also consider the bifurcation set for twice these values, namely for ηe = 0.08 and ηc = 0.22.
Finally, based on oceanic wave speeds calculated using TOPEX/POSEIDON satellite data in [60,61], realistic values of τw lie between 5.2 and 7.2 months, that is, in the range [0.43, 0.6] when scaled to years. This places the estimated range for at [0.80, 0.97]. Nonetheless, in the following section, we choose to consider the parameter range [0, 2] for two reasons: (i) to allow a comparison with previous literature [19,20,23,35]; and (ii) to understand the global origins of observed complex dynamics.
3. Bifurcation analysis
To demonstrate how the two types of state dependence affect the behaviour of model (1.1), we now present bifurcation sets in the -plane of (1.1) with (2.6) for different values of ηe and ηc. These consist of curves of torus bifurcations, saddle-node of limit cycle bifurcations and period-doubling bifurcations, all computed with the continuation package DDE-Biftool [36,37]. These curves are overlaid onto so-called maximum maps, where the maxima of long time series (after transients have died down) are computed and represented in colour for fixed and increasing c (where the previous solution serves as initial history); we typically compute a trajectory of 1000 years and disregard the first 500 years, where the integration step is 10−4 years. Maximum maps allow for a ready comparison with previous studies of the GZT model [19,23,35], and we fix throughout b = 1 and κ = 11 as in these previous works. Moreover, the bifurcation set in the -plane provides a clear overview of how the dynamics are organized more globally.
In §§3a we give a brief review of the findings in [35] for the constant-delay case ηe = ηc = 0. We then ‘switch on’ and increase the state dependence by presenting the bifurcation set for non-zero ηe and ηc. We first consider their effects individually in §§3b, and then in combination in §§3c. Finally, we briefly discuss in §§3d how the results could be interpreted in the context of the ENSO phenomenon.
(a) The constant-delay case
Figure 3 shows the bifurcation set in the -plane for the constant-delay case ηe = ηc = 0, overlaid onto the maximum map. The maximum map provides a convenient indication of the respective behaviour. Notice that sudden changes in colour, that is, in the maximum value of h on an attractor, align well with shown bifurcations. Periodic solutions born along the line c = 0 are stable and have non-zero amplitude for [62–64]; these are dominated by the effect of negative feedback and have a period of . On the other hand, for sufficiently large values of c there are stable periodic solutions that are dominated by the seasonal forcing. These solutions appear to be close to sinusoidal and have a period of one; for decreasing c, they lose their stability at a torus bifurcation curve. Between the line c = 0 and the torus bifurcation curve solutions are characterized by an interplay between negative feedback and periodic forcing, resulting in dynamics on invariant tori. As was mentioned in §1, the dynamics on the torus may be locked (periodic) or unlocked (quasi-periodic). The locked solutions are organized into p:q resonance tongues, where p/q is the winding number of the periodic solutions, of which there is generally a stable and an unstable one on the torus. Theory predicts an infinite number of resonance tongues, one for every rational winding number [65], and figure 3 shows a set of the largest resonance tongues. As could be expected from the mathematical form of model (1.1), the dependence of the solutions on the parameters appears to be highly nonlinear.
Figure 3. Bifurcation set with maximum map (a) and positive maximal Lyapunov exponent (b) in the -plane of equation (1.1) with (2.6) for the constant-delay case ηe = ηc = 0 in equation (2.6). Shown is a black curve of torus bifurcations, white curves of saddle-node bifurcations of periodic orbits that bound p:q resonance tongues, and grey curves of period-doubling bifurcations. The simulation scans here and throughout the paper are performed for fixed and increasing c. (Online version in colour.)
Inside some of the shown resonance tongues, period-doublings are found. This is a common feature when resonance tongues overlap sufficiently; for a detailed study of accumulating resonance tongues and associated complex behaviour see [66]. The resulting chaotic behaviour is of interest in the context of ENSO, since such behaviour shows a mixture of irregularity with locking to the seasons, which is seen as characteristic for El Niño and La Niña events in the real world. Note that period-doubling cascades and associated chaotic dynamics are found only in small regions of the parameter plane; see figure 3(b), which shows where the respective maximal Lyapunov exponent is positive, as calculated by following the algorithm for DDEs in [67]. Regions of positive maximal Lyapunov exponents are also observed for parameters where we show no curves of period-doublings. This is simply because there exist higher-order resonance tongues in between those we have shown and they are also expected to overlap. Furthermore, as demonstrated for a related DDE model in [20], chaotic regions are often bounded not only by period-doubling cascades, but also feature intermittent transitions that are characterized by the sudden appearance of chaos at a saddle-node bifurcation [68]. A more comprehensive discussion of the bifurcation set in figure 3 and associated other types of dynamics, including multi-stability, changing criticality of the torus bifurcation and folding tori, can be found in [35,69].
(b) Individual effects of the two types of state dependence
Considering only the state dependence of the upwelling delay τe, we find the anticlimactic result that it does not appear to have an influence on the observed dynamics. More specifically, for ηe = 0.04 and ηe = 0.08 and with ηc = 0 in (2.6), bifurcation set and maximum map in the -plane do not show any observable change on the scale presented in figure 3.
State dependence of the ocean mass adjustment alone, on the other hand, does have a dramatic effect. Figure 4 shows the bifurcation sets and maximum maps for fixed ηe = 0 and for ηc = 0.11 in figure 4a and for ηc = 0.22 in figure 4b. Already for ηc = 0.11, the difference with the case of ηc = 0 in figure 3 is striking: resonance tongues have become wider and overlap more, with larger regions where one finds period-doublings. The torus bifurcation curve now extends to c = 10 near the top of figure 4a, and there is a notable change in the size of the observed maxima. Moreover, there is a region with complex dynamics and additional saddle-node and torus bifurcation curves that entered the frame at large values of c above . These are caused by the state dependence — large c values mean large amplitude, which leads to larger fluctuations of the delay time τn. In fact, the scaling relationship between c and ηc as a result of equation (2.6) implies that these curves will diverge off to c → ∞ as ηc → 0; similarly, they move closer to the other bifurcation curves as ηc increases. The existence of non-trivial dynamics for large values of c is in stark contrast to the dynamics in the constant-delay case, where solutions for sufficiently large c are trivially dominated by the periodic forcing.
Figure 4. Bifurcation set and maximum map in the -plane of equation (1.1) with (2.6) for ηc = 0.11 (a) and ηc = 0.22 (b) and for fixed ηe = 0. Shown are black curves of torus bifurcations, white curves of saddle-node bifurcations of periodic orbits and grey curves of period-doubling bifurcations. (Online version in colour.)
Figure 4b, now for ηc = 0.22, shows a considerably more complicated bifurcation set compared with figure 4a. The bifurcation curves that emerge from the right-hand side of figure 4a have shifted towards lower values of c in figure 4b, such that they now interact with the pre-existing structure of resonance tongues. This gives rise to many more and overlapping resonance tongues that extend all the way to large values of c near . In this region, we find many and overlapping curves of subsequent period-doublings, which we refer to here as a period-doubling cluster. In this part of the -plane, there are not only period-doubling cascades inside resonance tongues but also period-doubling cascades of the 0:1 periodic solutions, that is, those dominated by the seasonal forcing. Also near , as well as near , we observe curves of saddle-node bifurcations of 0:1 periodic solutions, which meet at a cusp bifurcation at (which is not so easy to see due to the many surrounding bifurcation curves) and at , respectively.
Although the resonance tongues are not individually labelled, a careful comparison between figures 3 and 4 reveals that, with increasing state dependence, the roots of resonance tongues along the line c = 0 have moved to lower values of . Beyond our numerical approach and the scope of our analysis, the resonance tongues can be studied in the limit of small c. For example, the positions and boundaries of the resonance tongues could be estimated by means of a response or sensitivity function approach [70,71].
We remark that in a parameter region near in figure 4a,b the delay becomes negative during simulations; this is indicated by grey shading. To resume calculations for the next increment in c, we start the integration from a constant initial history of h(t)≡10−3. Similarly, the continuation of the lower period-doubling curve in figure 4b stops when the state-dependent component of the delay becomes so large that the delay becomes negative. The intricate boundary of the grey region of the -plane already indicates that it depends very sensitively on the initial history whether the delay becomes negative along a trajectory.
The overall resonance structure of the bifurcation sets in figures 3 and 4 is in some sense organized by the curves of torus bifurcations. In order to provide insight into how the bifurcation set changes qualitatively with ηc for fixed ηe = 0, figure 5 presents with 12 panels how the curves of torus bifurcations evolve in the -plane when ηc is increased gradually from zero. Here the colour scheme indicates the winding number α along the curve, which is given by the complex conjugate pair of Floquet multipliers e ± 2πiα that lie on the unit circle at the moment of torus bifurcation (also known as a Neimark–Sacker bifurcation). As ηc is increased, new curves of torus bifurcations enter the frame from the right [figure 5a–c]. They then interact with the torus bifurcation curve that exists already for ηc = 0 in several codimension-one singularities (of the surface of torus bifurcations in -space); see, for example, [72] for more information on singularity theory. The first such interaction is a saddle transition at ηc ≈ 0.115, where two curves meet at a single point and then reconnect differently (figure 5d,e). As a result, both curves of torus bifurcations extend to large c near . Two further saddle-transitions just past ηc ≈ 0.130 create two isolas, that is, closed curves of torus bifurcations [figure 5f –h]. Note that all these transitions through singularities of the surface are consistent with the colouring by winding number. Subsequently, these isolas shrink down to points and disappear in what we refer to as minimax transitions. Increasing ηc further, we find that a curve of period-doubling bifurcations (bottom right) enters the frame, to which the lower branch of torus bifurcations connects at a 1:2 resonance point [figure 5i–l].
Figure 5. Curves of torus bifurcations in the -plane in the region [2, 10] × [0.6, 1.8], where colour indicates the winding number. Here ηe = 0 and ηc takes the values (a) 0.080, (b) 0.095, (c) 0.110, (d) 0.115, (e) 0.120, (f ) 0.130, (g) 0.135, (h) 0.140, (i) 0.155, (j) 0.185, (k) 0.190 and (l) 0.200. (Online version in colour.)
Overall, figure 5 gives a good indication of how the bifurcation set becomes more complex with increasing ηc. Indeed, the torus bifurcation curves shown come with infinitely many resonance tongues, leading to the overall complexity that was shown in figure 4. Focusing on only the torus bifurcation curves avoids the considerable computational cost associated with computing many more bifurcation curves and maximum maps. We remark that computing all of the bifurcation curves shown in figure 4 is not fully automatic but requires an interactive process with DDE-Biftool that takes about 2–4 days per panel; computing the associated maximum map, on the other hand, can be fully parallelized, yet still takes about 4–5 h on the New Zealand eScience Infrastructure Pan cluster high-performance computing facility.
As was stated above, one of the effects of state dependence is that some of the resonance tongues become larger. To illustrate this effect in more detail, figure 6 shows how the 1:2 resonance tongue in the -plane changes for fixed ηe = 0 and increasing ηc. It demonstrates that the widening of some resonance tongues coincides with a breaking of symmetry. As discussed in [35], the constant-delay GZT DDE model features a -symmetry: in p:q resonance tongues with p or q even, there exist two distinct, symmetrically related locked periodic solutions. The 1:2 resonance tongue in figure 4a is such an example and there are actually two sets of curves of saddle-node bifurcations that lie on top of each other. As ηc is increased, however, they separate more and more, as figure 4b–d show. For the 1:2 resonance tongue, this symmetry breaking also results in a closed curve of period-doubling bifurcations being created at the right tip where the resonance tongue connects with the torus bifurcation curve; this curve of period-doubling bifurcations has two 1:2 resonance points where the respective branches of torus bifurcations connect. This phenomenon is actually part of the complicated bifurcation structure near 1:2 resonances predicted by theory; for example, see [65, section 9.5.3]. We remark that this bifurcation structure is also relevant for the period-doubling bifurcations of the seasonal forcing dominated solutions that emerge near ; compare with figures 5l and 4b.
Figure 6. Maximum maps and the 1:2 resonance tongue for ηc = 0 (a), ηc = 0.01 (b), ηc = 0.04 (c) and ηc = 0.11 (d) and for fixed ηe = 0. Shown are the white bounding curves of saddle-node bifurcations of periodic orbits, dark grey curves of torus bifurcations and light grey curves of period-doubling bifurcations. For the background colour scheme, see figure 4. (Online version in colour.)
(c) Combined effects of the two types of state dependence
We now study the combined effects of both types of state dependence, that is, we consider ηe > 0 and ηc > 0. Figure 7 demonstrates that state dependence due to upwelling has a large effect on the overall bifurcation set when the ocean-adjustment state dependence is also present. This is surprising since state dependence due to upwelling on its own does not appear to change the bifurcation set.
Figure 7. Bifurcation set and maximum map in the -plane of equation (1.1) with (2.6) for ηc = 0.11 and ηe = 0.08 (a) and ηc = 0.22 and ηe = 0.04 (b). Shown are black curves of torus bifurcations, white curves of saddle-node bifurcations of periodic orbits and grey curves of period-doubling bifurcations. (Online version in colour.)
Figure 7a considers the effect of the maximal nominal value ηe = 0.08 of upwelling on the case of a smaller ocean adjustment of ηc = 0.11 shown in figure 4a. The comparison shows that figure 7a features a dramatic change in the way the different types of solutions are distributed across the -plane. In particular, the torus bifurcation curve appears to have split into two segments around . In fact, as we checked but do not show here, even for c > 40, the two curves remain disconnected. This change of the bifurcation set is the result of a transition similar to the interacting torus bifurcation curves in figure 5. Note in figure 7a how some resonance tongues have become significantly wider and overlap to a greater extent, resulting in larger regions of chaotic behaviour. We observe again the emergence of the period-doubling cluster near , as well as the cusp bifurcations of the forcing dominated solutions near and near .
Figure 7b considers the effect of the nominal value ηe = 0.04 of upwelling on the case of a maximal ocean adjustment of ηc = 0.22 shown in figure 4b. The changes are not quite so pronounced, because the bifurcation diagram for ηc = 0.22 and ηe = 0 is already quite complex with torus bifurcation curves extending to large c. However, note that in figure 7b, the period-doubling cluster near is now found for smaller values of c.
Finally, figure 8 shows the bifurcation set for the maximal nominal ocean adjustment of ηc = 0.22 and the maximal nominal upwelling of ηe = 0.08. Here figure 8a shows the bifurcation set with the maximum map and figure 8b shows it with a map of where the maximal Lyapunov exponent of the respective solution is positive, that is, where the dynamics is chaotic. Comparing figure 8 with figure 3 clearly drives home the point that state dependence has a large effect on the bifurcation set of model (1.1) and, hence, on the observable dynamics. When comparing figure 8a with figure 7b, we observe that increasing the upwelling to ηe = 0.08 results in the two large period-doubling bifurcations of the 0:1 solutions to undergo a saddle transition near τn = 0.5 and connect differently. Moreover, the period-doubling cluster now appears for even lower values of c. As figure 8b shows, chaotic dynamics is associated with this period-doubling cluster. Indeed, comparison with figure 3b clearly demonstrates that state dependence results in considerably more and larger regions where chaotic dynamics can be found. Interestingly, there is also a large region of chaotic dynamics for very low values of , near the boundary where the delay becomes negative.
Figure 8. Bifurcation set and maximum map (a) and positive maximal Lyapunov exponent of solutions (b) in the -plane of equation (1.1) with (2.6) for ηc = 0.22 and ηe = 0.08. The small circle in (a) indicates the parameter choice for the time series in figure 9. (Online version in colour.) Figure 9. (a1) shows a time series of the Nino3 index calculated from the observational dataset NOAA Optimum Interpolation SST V2 (January 1982 – December 2017); the data is detrended linearly. (a2) shows the corresponding power spectrum in arbitrary units [a.u.] calculated with the Welch method. (b1) and (b2) show the same for an example attractor of equation (1.1) with (2.6) for ηc = 0.22 and ηe = 0.08 at the parameter point , as indicated in figure 8, and starting from the constant initial history h≡1. Data for row (a) is provided by the Physical Sciences Division, Earth System Research Laboratory, NOAA, Boulder, Colorado. (Online version in colour.)
(d) Interpretation for the ENSO system
By introducing a physically motivated non-constant delay time, we have made the DAO model more detailed. This model modification has led to dynamics that are expected to represent realistic aspects of the ENSO system. In the constant delay case, irregular (chaotic) behaviour could only be found for small pockets of the -plane and, as such, could not be considered a prominent feature of the model behaviour. In the state-dependent delay case, on the other hand, this type of behaviour is more prominent. Of particular interest here is the period-doubling cluster, which is the relatively large region with chaotic dynamics near . Note that our discussion of suggested nominal values of the parameters , and τw in §2 places the period-doubling cluster within the range of physically relevant values of .
Since the emergence of chaotic behaviour in ENSO models has often been attributed to specific routes to chaos (i.e. the period-doubling [21,22,73], quasi-periodic [17,24,74] and intermittent [75] routes), we mention here that all these routes may be observed within the period-doubling cluster, depending on the chosen path through parameter space.
To show that the associated dynamics is indeed relevant for ENSO, figure 9 displays time series and power spectra for an observational dataset and an example attractor from model (1.1). Figure 9a1 shows a time series for the Nino3 index, which is the averaged SST over 5°N–5°S and 150°W − 90°W. It is calculated from the observational dataset NOAA Optimum Interpolation SST V2 (January 1982 – December 2017) based on the analysis of in situ and satellite data [76]. Typically, for Nino3, the data are detrended by using a five-month running mean. Here, we simply detrend linearly in order to preserve the seasonal cycle. The power spectrum in figure 9a2 is calculated by using the Welch method with windows of length 15 years and overlapping across 12 years. The time series in figure 9a1 demonstrates characteristic features of the ENSO system, including the irregularity of very distinct large peaks, representing El Niño events, both in terms of amplitude and occurrence. There is clearly a high degree of variability and, as we checked, the large peaks tend to be seasonally locked. The influence of the seasons is evident in the power spectrum in figure 9a2 with a distinct peak centred at 1 year. There are also two distinct peaks centred near 1.5 and 3.5 years, representing the time scales of ENSO variability. Note that the peak centred near 3.5 years is particularly broad.
The time series of the thermocline deviation h of model (1.1) in figure 9b1 lies on a chaotic attractor that resides in the period-doubling cluster at . It features irregularity in the form of relatively large peaks that occur every 4–8 years, as well as seasonal locking with a corresponding dominant frequency in figure 9b2 centred at 1 year. The power spectrum also displays a broad peak centred near 2.5 years. Comparison between rows (a) and (b) of figure 9 demonstrates that these two characteristic features of El Niño events from observational data (and many other models [18,24,55,77,78]) are captured by model (1.1). We find that the seasonal locking is very robust with respect to the choice of parameters. Number and distribution of El Niño peaks, on the other hand, depend on the choice of the parameter points in regions of overlapping p:q resonance tongues. Furthermore, the relative strengths between the peaks in the power spectrum, representing both seasonal and El Niño time scales can be influenced by the choice of the seasonal forcing strength c. We remark with regard to the amplitude of the oscillations in figure 9b1, that one could calibrate the remaining model parameters b and κ in order to scale the variable h. Importantly, by using figure 8 as a guide, solutions demonstrating fundamental ENSO characteristics can be found readily in the appropriate regions of parameter space. Therefore, we conclude that state-dependent delay associated with upwelling and ocean adjustment must be considered as a contributing mechanism for the observed irregular but seasonally locked behaviour of the ENSO phenomenon.
In spite of the good agreement of the power spectra in figure 9, row (b) does not feature the same degree of variability found in the observational data in row (a). Of course, the example solution represented in row (b) is calculated for the thermocline depth anomaly with the conceptual model (1.1) with only a negative delayed feedback and forcing term. Exactly how the variable h translates to an observable such as Nino3 is not obvious, especially since Nino3 is a spatially averaged quantity and model (1.1) has no spatial component. Making a closer comparison between the dynamics of model (1.1) and climate observables is left for future work.
4. Summary and discussion
The processes controlling the irregularity of the temporal variability of ENSO are still not satisfactorily determined and this hampers skilful prediction [8,9]. Certainly atmospheric noise, in particular, westerly wind bursts, contribute to the irregularity of the ENSO variability and some literature is in support of the theory that ENSO is a stable, linear system subject to stochastic forcing [79–82].
However, there are also quite a few deterministic mechanisms that can give rise to chaotic behaviour, such as bursting phenomena [83,84], period-doubling [21,22,73] and nonlinear resonances with the seasonal cycle [17,24,74].
This motivates the study of conceptual ENSO models, which have been very useful in the past for investigating the interactions of known mechanisms and for discovering new physical mechanisms of variability.
The work presented here makes a contribution in this spirit by considering a quite simple conceptual ENSO model in the form of a DDE. Models of this type are known to efficiently represent physical processes, such as wave propagation and feedbacks [15]. Moreover, DDE models lend themselves to an investigation with advanced tools. As our study demonstrates, this is also the case when the delays that arise are no longer considered to be constant, which is an assumption in much of the previous literature in climate science and other fields. In particular, it is possible to ascertain whether and how state dependence of delays contributes to the overall dynamics of the feedback-driven system.
More specifically, we performed a study of the influence of additional and state-dependent delays arising from the interaction between the thermocline and the ocean surface. To this end, we incorporated the state dependence of delays into a well-studied ENSO DDE model with only negative feedback and periodic forcing [19,35,69]. We derived state-dependent delay terms for upwelling in the eastern Pacific and for ocean adjustment in the central Pacific, and then investigated their effect individually and also in combination. To capture important features of the overall dynamics, we considered the respective bifurcation sets in the -plane of the forcing strength c and the constant portion of the overall delay. We demonstrated how an increasing complexity of the bifurcation set arises with increasing state dependence. Relevant levels of state dependence of the delays were shown to severely enhance resonance phenomena, leading to considerably larger regions with overlapping resonance tongues and associated period doublings and chaotic dynamics. The latter generate times series with characteristic features of the ENSO phenomenon, specifically, its irregular occurrence every few years in combination with seasonal locking. We remark that chaotic behaviour arising from nonlinear resonances has always been considered to be a less important mechanism of ENSO variability because the observed windows of its existence in parameter space were relatively small [25]. The widening of chaotic regions in the parameter space, as found here due to state-dependent delays, is hence highly relevant for attributing irregular ENSO variability more to deterministic processes rather than noise.
It is worth noting that other mechanisms may play a role in the complexity of the model dynamics. For example, the authors of [21] showed that the inclusion of additional Rossby wave modes and an asymmetric coupling function each increase the complexity of the dynamics of their ENSO delay-difference model. Furthermore, including positive delayed feedback due to Kelvin waves and asymmetry of the coupling each increase the complexity of the bifurcation diagram of the GZT DDE model studied here [20]. It is an interesting observation that in the state-dependent DDE model the breaking of symmetry, discussed in §§3b, is achieved without the need of introducing an asymmetric coupling function, as was the case in [18,20,21].
The bifurcation sets presented here highlight the role of invariant tori in producing dynamical mechanisms resulting in complicated behaviour. The schematic bifurcation diagram in [85], based on simulation runs of a hybrid coupled general circulation model, also indicates the role of torus dynamics in ENSO. The origins of dynamics on a torus in an intermediate coupled ocean-atmosphere model for the ENSO system was identified more thoroughly with continuation methods in [86]. A bifurcation set for the DDE ENSO model introduced in [18] is presented in fig. 6 of [87]. Like the bifurcation sets presented here, the authors of [87] found resonance tongues emanating from a curve of torus bifurcations; however, these resonance tongues appear not to interact with each other (at least over the parameter ranges considered). The bifurcation sets for a conceptual ENSO model presented in [83,84] do not feature resonance tongues but still demonstrate the emergence of chaotic behaviour due to period-doubling. We remark that most of the dynamics we reported on involve two-tori and their further bifurcations and could, hence, potentially be observed in models with dimension as low as three. On the other hand, some more complicated dynamics might well be of higher dimension; a closer analysis of this is left for future work.
From a more general point of view, we demonstrated that state dependence of delays can have a serious impact on the dynamics, even when the constant-delay DDE already features complicated nonlinear behaviour. Our results follow on from the study by Calleja et al. [44], which showed that state dependence of delays alone can create surprisingly complicated nonlinear dynamics when the constant-delay DDE has only trivial, linear dynamics. Together, these results must be seen as a ‘health warning’ that delays should not simply be assumed to be constant as a matter of course. Given the simple form of the model we considered, with only delayed negative feedback and periodic forcing, we expect that our case study will be of interest to researchers from various other areas of application. The results presented here should provide motivation to investigate state-dependent delays more rigorously. A concrete example concerns the dynamics driven by state-dependent delays in lasers with frequency-dependent feedback [50].
In the specific context of climate modelling, there are also several avenues for future work. The work presented here is an initial step in the analysis of the effects of state dependence on ENSO. For simplicity, we have taken the state dependence of the delay to be linear. An interesting question from the modelling perspective is whether state-dependent delays can be derived explicitly from more complicated ENSO models of an intermediate complexity and, in particular, what mathematical form would they have. This paper might provide motivation for such a more quantitative analysis. One of the tools that can be used in this context is the Mori–Zwanzig formalism [88–90], which is a formal model reduction technique. Another focus for future work will be to study other physical mechanisms that may lead to state-dependent delays in the ENSO system. For example, a contributing factor to the delay times of the DAO paradigm is the position of the zone of ocean-atmosphere coupling in the central Pacific Ocean, which is closely associated with the western Pacific warm-pool. In fact, a state-dependent DDE for the position of the western Pacific warm-pool has already been suggested in [91], but has not yet been analysed in much detail. Finally, it will be of practical relevance to study the effect of atmospheric noise on the ENSO behaviour in DDE models with state-dependent delay.
Data accessibility
Not applicable.
Authors' contributions
A.K. carried out the numerical computations; all authors conceived of and designed the study, and drafted, read and approved the manuscript.
Competing interests
The authors declare that they have no competing interests.
Funding
H.A.D. acknowledges support by the Netherlands Earth System Science Centre (NESSC), financially supported by the Ministry of Education, Culture and Science (OCW), Grant no. 024.002.001.
Disclaimer
Not applicable.
Acknowledgements
We thank Jan Sieber for advice on implementing DDE-Biftool for the state-dependent delay case and Tony Humphries for advice on numerical integration of state-dependent DDEs. We are grateful to the Centre for eResearch at the University of Auckland for their help in facilitating the use of the New Zealand eScience Infrastructure high-performance computing cluster.
