Taming of the Hopf bifurcation in a driven El Niño model

1 Introduction

El Nino generally represents the extensive warming of the central and eastern Pacific Ocean. In 1986, Geoffrey K. Vallis studied climatic change occurs in east-central equatorial Pacific El Niño phenomenon. During El Niño phase, the trade winds deteriorate in the central and western Pacific leading to a depression of the thermocline in the eastern Pacific and cooling the surface, cutting off the supply of nutrient rich thermocline water. Hence, a rise in sea surface temperature (SST) and a severe decline in primary productivity happen. During non-El Niño phase, the west receives trade winds through the tropical Pacific piling up warm surface water in the west Pacific. The sea surface becomes higher at Indonesia than at Ecuador, and the SST reaches higher in the west. The cool temperatures off South America are because of an upwelling of cold water from deeper levels, which is nutrient-rich, supporting high levels of primary productivity, diverse marine ecosystems [18], [20], [21]. The El Niño phenomenon is a very complex and irregular one, modeling its intricate dynamical behavior is still a challenge for scientific community [12], [13], [14], [15].

Mathematical model that describes El Niño phenomena is the continuous-time Vallis model [1], which consists of a set of three autonomous first-order nonlinear ordinary differential equations. Borghezan and Rech [11] investigated chaos [9] and periodicity for this model, also reported the existence of periodic structures embedded in a chaotic region.

Vallis model [1], [2] is still considered as the simplest mathematical model for El Niño. The localization problem of compact invariant sets of nonlinear time-varying systems with the differentiable right-side is investigated using the Vallis model [3]. The sufficient conditions for the existence of periodic solutions and their stability nature are discussed in [4]. Most of the earlier discussions about chaos in Vallis model use the parameter conditions P = 0 and, in [5], authors numerically investigated the existence of chaotic solutions for P = 0 and P = 0.83. In [23] Kurgansky et al. studied the long-term climate variability using nonlinear chaotic model and emphasize the influences of parameter variations.

Driven oscillators are used widely to investigate the nonlinear behavior of interdisciplinary systems [10]. The solution for the driven oscillators holds two parts, a transient part and a steady-state part, which must be used together to fit the physical boundary conditions of the problem. In 2006, Moore et al. investigated optimal forcing patterns for coupled models of ENSO and revealed the significance of excitation terms in prediction of climatic variation [17]. Many biological rhythmic processes can be modeled by nonlinear differential equations exhibiting limit cycle behavior (that is they admit Hopf bifurcation). Hence, it is important to study such phenomenon in order to get deeper knowledge about the systems.

In the present work, we aim to study the effect of driven force on the Hopf bifurcation. For this purpose, we consider a Driven Vallis Model (DVM) and made the bifurcation analysis with respect to an amplitude of the external force. We find that in the absence of driven force, we can observe either supercritical or subcritical Hopf bifurcation whereas the presence of driven term leads to existence of both supercritical and subcritical Hopf bifurcation.

2 Driven Vallis Model (DVM)

The El Niño phenomenon [24] is about the band of warm ocean water that develops in the central and east-central equatorial Pacific which has great impact on the global climate.

The generalized mathematical model is defined as

where, u is the velocity of the ocean surface flow, T_w and T_e are the relative temperatures at the western and eastern edges of the ocean basin, respectively. A⁻¹ is the relaxation time of the temperature, B is the coupling factor between the temperature difference and the ocean surface flow C is the coefficient associated with internal friction in ocean water l is the ocean basin width T₀ is the relative temperature deep in the ocean u∗ is the velocity of the tradewind u∗ = u₀(1+sin ωt) T∗ is the steady-state ocean temperature for u∗ = 0 [16, 19].

There are many mathematical models available for El Niño [15]. Kuhlbrodt et al. [22] simplified the mathematical model into continuous-time Vallis model [1], [2] defined by the non-dimensional model.

There are many mathematical models available for El Niño and one of these simple models is the continuous-time Vallis model [1], [2] defined by the non-dimensional model,

For the parameter values, P = 0, b = 103, c = 3, system (1) shows chaotic oscillations [3], [4], [5], [6]. In this paper, we discuss the dynamical behavior of a driven Vallis model as in (2).

where f(t) = 1+f cos(ωt). Using a computer search algorithm, we find the parameters for chaotic solutions to be b = 90, P = 0.5, c = 1, ω = 4, f = 3.6 and the 3D phase portraits of system (2) for initial conditions [1.4, 0.1, 0.3] are shown in Figure 1a and the 3D Poincare sections for z = 0 are shown in Figure 1b.

Figure 1:

(a) 3D phase portraits of the Vallis model; (b) 3D Poincare sections where the red markers shows the phase trajectory crossing the ‘z = 0’ axis while section crossing ‘+’ to ‘−’ and black markers shows the phase trajectory crossing the ‘z = 0’ axis while section crossing ‘−’ to ‘+’.

Most of the dynamical analysis of the Vallis model is done assuming that the average effects of equatorial winds (P) are negligible or zero [3], [4], [5], [6], [7] until it was proved that the parameter P has a considerable effect on the dynamical behavior of the Vallis model (2). Hence, in this paper, we investigate the effect of the parameter P on the bifurcation of the system and also other dynamical properties.

3 Dynamical properties of driven Vallis model

3.1 Bifurcation

Bifurcation plots are derived and presented to show the impact of the system parameters on its dynamical behaviors. We derive the bifurcation plots for three different cases in Table 1:

Table 1:

Different cases for bifurcation.

Case name	Control parameter range	Other parameters
Case-A	0.5 ≤ ω≤ 5	b = 90, P = 0.5, c = 1, f = 3.6
Case-B	2 ≤ f≤ 13	b = 90, P = 0.5, c = 1, ω = 4
Case-C	0.5 ≤ c≤ 2.75	b = 90, P = 0.5, f = 3.6, ω = 4

3.1.1 Case-A

As described in Table 1, for deriving the bifurcation plots for case-A, we choose ω as the control parameter. The initial condition for the first iteration is taken as [1.4, 0.1, 0.3]. The ω is varied from 0.5 to 5, we could see that the system shows chaotic regions with both period doubling and period halving routes. The creation of period doubling followed by their annihilation via period-doubling bifurcation (i.e. antimonotonicity) is an important behavior which is useful for studying intricate behavior of the chaotic system. It should be noted that forward period-doubling bifurcation sequences followed by reverse period doubling sequences is a unusual behavior. The process of period doubling and period halving occurring in a bifurcation diagram of a system is termed as antimonotonicity.

To show the existence of multistability, we use a robust way to plot the bifurcation plots where the initial conditions are changed in every iteration to the end values of the state variables. Figure 2 (blue) shows the forward continuation where the parameter ω is increased from minimum to maximum and Figure 2 (red) shows the backward continuation where the parameter ω is decreased from maximum to minimum and the local maxima of the state variables are plotted. The respective finite time Lyapunov exponents (LEs) are calculated using Wolf algorithm [8] for run time of 40,000 s with the initial conditions changed as like for Figure 2b are shown in Figure 2c.

Figure 2:

(a) Maximum of driven Vallis model (1) with forward (red) and backward (black) continuation. (b,c): LEs for b-forward continuation, c-backward continuation.

3.1.2 Case-B

For deriving the bifurcation plots for case-B, the control parameter is chosen as the amplitude of the forcing term (f). We could observe pocket of chaotic behavior for the range of 2 < f < 13. Antimonotonicity behavior observed in bifurcation plot (Figure 3). We use only the forward continuation to plot the bifurcation and as shown in Figure 3a, the DVM system takes a period doubling route to chaos. The Lyapunov spectrum for the parameter f is presented as in Figure 3b to confirm the existence of chaotic oscillations.

Figure 3:

Bifurcation of Vallis model with parameter f (a) and the corresponding LEs (b).

3.1.3 Case-C

In this case, the control parameter for bifurcation is the parameter c which is varied between the range [0.5, 2.75]. We could observe chaotic behavior for the parameter range 0.8 < c < 1.18 and 1.5 < c < 1.95. Again, only forward continuation is used to derive the bifurcation plots and the system takes period halving exit from chaos as shown in Figure 4a. The corresponding Lyapunov spectrum is presented in Figure 4b.

Figure 4:

Bifurcation of Vallis model with parameter c (a) and the corresponding LEs (b).

4 Tailoring of a Hopf bifurcation region through the variation of a parameter

To illustrate the dynamical transition and their stability, we plotted one-parameter bifurcation diagram using XPPAUT software for the variable x by varying the forcing term f(f(t) = Asin(ωt)) in Figure 5. The black filled circles and red filled circles denote the stable and unstable steady state, respectively. The blue unfilled circles denote the unstable oscillatory state which is also co-exist with unstable steady state. Figure 5 clearly illustrates the dynamical transitions as a function of forcing. We noticed that the stabilization stable fixed point that is transition from unstable steady state to stable steady state arises through inverse subcritical Hopf bifurcation at f = −2.57. Further stable fixed point losses their stability (i.e. the transition from stable steady state to unstable steady state) through subcritical Hopf bifurcation when f = 2.77.

Figure 5:

One parameter bifurcation diagram for DVM model as a function forcing “f”. Here, HB is the Hopf bifurcation point. Red filled circles and black filled circles denote the unstable and stable steady states, respectively. Blue unfilled circles represent the unstable oscillatory state. Other system parameters are fixed as b = 90, P = 0.5, c = 1.0.

4.1 Two-parameter bifurcation diagram

Further, to understand the dynamical transition more clearly, as a function of forcing f and system parameter c, we have plotted two-parameter diagram in (f, c) space in Figure 6. The US and SS are the unstable and stable steady states, respectively. At lower values of c, we found broader region of the stable steady-state. The observed stable steady state region decreases with increasing the region of unstable steady state when increasing the values of c to higher ranges. At larger values of c, the stable steady-state suppresses completely and the entire region is accompanied by the unstable steady-state. From the analysis, we can observe that the higher range of control parameter c, the system becomes completely unstable. As we found the one parameter bifurcation plot, the unstable steady state gets stabilized (US to SS ) via inverse subcritical Hopf bifurcation as a function forcing ‘f’. Subsequently, the stable steady gets destabilized (SS to US) occurs through subcritical Hopf bifurcation.

Figure 6:

The two parameter diagram in (f, c) space. SS and US are Stable steady State, Unstable steady states and HB Hopf bifurcation, respectively. Other parameters are same as in Figure 5.

4.2 Phase portraits

To understand the oscillatory behaviors we have plotted time evolution and phase portraits of such state in different planes. Firstly, time evolution for x, y and z variables are plotted as function of time in Figures 7(a)–(c) and the emergence two attractors in each variables are evident. To validate this we also plotted the phase portraits in x-y, y-z and z-x planes which is shown in Figures 7(d)–(f). The phase portraits clearly shows the double scroll attractor, that is the system scrolling between the two attractors.

Figure 7:

Time evolution of oscillatory state for driven Vallis model in terms of (a) x variable, (b) y variable, and (c) z variable. Phase portraits for (d) x-y plane, (e) y-z plane, and (f) x-z plane. Other parameter are same as in Figure 5.

5 Conclusion

We have derived the dynamical properties of the driven Vallis model. A new method of taming the Hopf bifurcation points are presented and by tuning the driving force the system changes its state from stable steady state to oscillatory state and vice versa. We find that in the absence of driven force, we can observe either supercritical or subcritical Hopf bifurcation whereas the presence of driven term leads to existence of both supercritical and subcritical Hopf bifurcation.