On the existence, stability and chaos analysis of a novel 4D atmospheric dynamical system in the context of the Caputo fractional derivatives

1.   Introduction

Over the last three decades, a great deal of research has been done on the fascinating nonlinear phenomenon known as chaos. Mathematical modeling, which promotes deeper insights and well-informed decision-making, is regarded as a very effective method of formulating and addressing real-world problems. It has great potential in many fields, including biomedical engineering applications to the human brain and heart ^[1], comprehensive liquid mixing with low power consumption, circuit systems ^[2], collapse prevention of power systems, chemistry, physics ^[2,3], engineering, and climate science. Within the field of climate research, mathematical models facilitate the simulation and comprehension of intricate climate systems, empowering scientists to forecast and examine climate trends. For example, the Lorenz-84 system illustrates the irregularity of the atmosphere ^[4], and Yong-Guang's study examined a low-order model for Hadley circulation ^[5]. The complexity and intrinsic unpredictability of weather and climate systems have been acknowledged by chaotic models, which have revolutionized our understanding of atmospheric dynamics. Although challenges remain advancements in computer methods and multidisciplinary studies have the potential to improve our capacity to simulate and forecast atmospheric events. Future atmospheric research will be guided by the insights provided by chaos theory, which provides a more profound understanding of the complex and frequently unpredictable nature of the Earth's atmosphere. The surface of the Earth and the lower atmosphere are warmed by the absorption and reradiation of heat by atmospheric gases such as carbon dioxide ($CO_2$), methane ($CH_4$), and water vapor ($H_2O$). The term "greenhouse effect" refers to this mechanism of heat absorption and radiation. The modeling and control of greenhouse gases (GHGs) have attracted the attention of numerous scholars in recent years, leading to the development of various control strategies.

The concept of fractional differentiation ^[6] is 300 years old; however, in the last 50 years, it has gained a lot of attention and discussion. One of the most popular trends in applied mathematics has been the investigation of fractional calculus ^[7,8]. Due to its more accurate explanation of physical occurrences. It has been recently shown that approaches such as integer-order differential operators are not necessarily appropriate for describing complicated and nonlinear phenomena. Fractional calculus has increasingly been applied to describe a wide range of real-world events in fields like economics, biology ^[9], chemistry, control theory, mechanics, signal and image processing, and electricity. Due to the limitations of existing methods, scientists and mathematicians are aiming to create complex and rigorous mathematical operators that can accurately mimic and capture these natural phenomena. In this sense, the power law and the definition of local differential operators have been employed.

In recent years, researchers have shown interest in the synchronization and control of chaotic systems ^[10] with fractional-order dynamics. This led to numerous notable contributions to the Chua and Liu's circuit of chaos synchronization ^[11,12,13]. Many fractional-order chaotic systems have been effectively controlled using a variety of techniques such as adaptive controller ^[14,15], optimal controller ^[16,17], or sliding mode controller (SMC) ^[18,19,20], etc. The SMC is one of the most popular methods to control the chaos. In ^[21], authors applied adaptive SMC for uncertain fractional reaction-diffusion systems. In ^[22,23], SMC was used for fractional reaction-diffusion neural networks.

Understanding and projecting Earth's climate system requires the use of climate modeling, which also offers important insights for adaptation and mitigation strategies related to climate change. Scientists can estimate future climate scenarios and evaluate the possible implications of human activity using these models, which mimic and analyze the intricate interactions between numerous atmospheric, marine, and land surface systems ^[5,24]. Inspired by meteorological models and other forms of abstract models such as a chaotic atmosphere propagation model ^[25], atmospheric circulation systems ^[26], a model of atmospheric dynamics of carbon dioxide gas ^[27], post-quantum model of Navier problem ^[28], a fractal-fractional model of polluted lakes ^[29], financial systems in the sliding modes ^[30], memristor-based chaotic model of circuits ^[31], quantum-based model on integro-differential system ^[32], hybrid comformable initial value problems ^[33], and fractal-fractional model of HIV-1 infection ^[34], we have developed an atmospheric model describing the interaction between westerly currents, temperature, and GHGs in the Caputo fractional domain. To conduct the study, we analyze chaos using fractional-order commensurate and non-commensurate derivatives through analytical techniques and numerical simulations. A fractional derivative is introduced to obtain more details about the system's chaotic behavior. In the presence or absence of uncertainties and outside disturbances, we formulate sliding mode control (SMC) laws to manage chaos in the proposed atmospheric fractional-order system. We evaluate the controllers' efficiency by integrating them into the westerly current dynamics.

The rest of the paper is organized as follows: In Section 2, we give some basic definitions, lemmas, theories, and notations related to fractional calculus. Section 3 provides an extensive overview of the proposed fractional-order atmospheric model. In Section 4, we discuss the existence and boundedness of solutions, equilibrium points and their stability, and the conservative behavior of the system. Section 5 presents the detection of the chaos behavior of the system with the help of the Lyapunov exponent with graphical representations of the Poincaré maps and bifurcation analysis. Section 6 focuses on designing the SMC for chaotic systems and analyzing their stability. In Section 7, we provide a numerical example of the suggested work. Section 8 draws a conclusion for the current work.

2.   Fundamental concepts

Definition 2.1 (^[6]). The Caputo derivative of the continuous function $s(t)$ of fractional order $m-1 < \gamma < m$ is defined as

where $\Gamma(\cdot)$ is the Gamma function.

Definition 2.2 (^[6]). For the integrable function $s(t)$, the Riemann-Liouville fractional-order integral is defined as

Lemma 2.3 (^[35]). Let $s(t)$ be a continuous function on $[t_0, \infty]$ that satisfies

here, $0 < \gamma\le$1, ($\alpha$, $\zeta$) $\in\mathbb{R}^2$, and $\alpha\neq0$, then

Lemma 2.4 (^[6]). Consider the system

with $\gamma\in(0, 1]$, $s : [t_0, \infty)\times \Gamma\rightarrow \mathbb{R}^n$. When $s(t, h)$ holds the local Lipschitz criteria regarding $h$ on $v$, the system has a unique solution.

3.   Model formulation

A mathematical representation of atmospheric circulation that incorporates a condensed version of Earth's topography and climatic patterns was proposed by Edward N. Lorenz in 1984, becoming popularized as the Lorenz-84 climate model. The globe-encircling wind, also known as the Hadley cell circulation ^[5], is an important atmospheric circulation pattern that transfers heat from the equator toward the poles, being propelled by their temperature differential. This circulation can be observed in the Lorenz-84 model using a global wind component. In this model, the globe-encircling wind interacts with the thermal differential between the equator and the poles to provide poleward heat movement. As heat gets distributed with the wind, the westerlies get stronger due to greenhouse effect-induced global warming. In ^[25], authors demonstrated chaos for the fractional order model (3.1).

Figure 1 demonstrates the dynamics of the model (3.1) for $\beta = 1, \, \beta = 0.95, \, \beta = 0.85,$ and $\beta = 0.75$ for the parameter values $a = 0.24, \, b = 11$, $P = 3$, $M = 8,$ and the constant $CO_2$ input $F = 1.66$.

In the present work, instead of using a constant input of $CO_2$, we have incorporated the radiative forcing of GHGs as a new equation to investigate their continuous impact on the Lorenz-84 climate circulation system ^[36]. We aim to examine the change in the strength of the westerly current with the continuous change in growth rate of GHGs in the atmosphere. Therefore, we have installed an additional equation representing the radiative force of GHGs; its effect on westerly current is represented by the incorporation of a term representing the continuous impact of GHGs in the first equation. The Caputo fractional derivative framework is utilized in this study to further examine the atmospheric system. The proposed fractional differential equation (FDE) is:

where $\nu$ is the force of the symmetrically circling westerly wind current around the Earth. Heat is transported poleward by a series of superposed large-scale currents, the cosine and sine phases, that are represented by $\rho$ and $\psi$. $\mu$ represents the GHGs, $M$ is the symmetric thermal forcing brought on by the temperature differential across latitudes, and $P$ is the difference in temperature between the ocean and the continent. The wave advection force that results from the westerly current is expressed as $b$. $n$ is the rate of GHGs added to the westerly current, and $r$ is the natural depletion rate of GHGs that may happen due to chemical reactions.

4.   Theoretical dynamics

4.1.   Existence and uniqueness of the solution

Our investigation aims to establish whether the solutions to the fractional four-dimensional chaotic atmospheric model exist and are unique under the Caputo fractional derivative framework. By careful exploitation of mathematical details, we have confirmed that the solutions not only exist but are unique, thus ensuring the validity of model (3.2) for further analysis and applications.

Theorem 4.1. In the region $\mathfrak{B} \times [0, T],$ where $\mathfrak{B} = \{(\nu, \rho, \psi, \mu) \in R^4:max\{|\nu|, |\rho|, |\psi|, |\mu| \}\le F\},$ and $T < + \infty$, there exists a unique solution of the system (3.2).

Proof. Let us consider

and a function

where

Here, $V(t, U)$ is defined on $B\times[0, T]$ and let $N = sup_B ||V(t, U)||.$

We interpret the norm as $||U(t)|| = sup_{t\in[0, T]}|U(t)|.$ We shall demonstrate the existence of some $\Phi$ such that

Consider that

This implies

where

Consider

This gives

Let $\Delta$ represent the Picard's operator that is being constructed through the utilization of the fractional integral and function V. Consequently, we get the following equation

It is necessary to demonstrate that this operator is a contraction mapping and that it maps a whole non-empty metric space into itself.

Consider

where $m$ is a constant. Using (4.5) as the norm, we obtain

The above inequality in (4.6) holds if

We now isolate a condition in which the operator $\Delta$ is contraction. We take the following actions to fulfill this requirement:

The preceding formula demonstrates that when the relationship is defined by the Picard's operator $\Delta$, it turns into a contraction, and

In this way, it is proven that the Picard's operator $\Delta$ is a contraction. We determine that the operator $\Delta$ has a unique fixed point by applying the Banach fixed-point theorem. As a result, there is only one solution for the FDE (3.2). □

Here, note that

Based on the model's parameters, we can interpret $\Phi_i$, where i = 1, 2, 3, 4, as follows:

$\Phi_1$ = Strength of the westerly current due to the impact of sinusoidal waves;

$\Phi_2$ = Strength of the sine waves;

$\Phi_3$ = Strength of the cosine waves;

$\Phi_4$ = Total amount of greenhouse gas present in the armosphere;

All these $\Phi_i$, i = 2(1)4, indicate the strength of the westerly current.

(1) The combined effect of the sine and cosine waves maximizes the strength of the westerly current if $\Phi_1$ is at its maximum.

(2) Suppose that $\Phi_4$ represents the maximum. If the depletion rate increases in the atmosphere, then the rate of GHGs reduces, and due to this effect, the strength of the westerly waves decreases.

4.2.   Boundedness

In this section, we establish that the solutions of system (3.2) are bounded.

Theorem 4.2. The solutions of system (3.2) are bounded uniformly.

Proof. Define a function,

Utilizing the Caputo fractional derivative, we obtain

There is a unique and existing solution in

The inequality stated above gives that

By Lemma (2.3), we have

Clearly,

as t$\rightarrow$ $\infty$.

Thus, the system (3.2) that initiates in $\mho$ and its whole solution remained bounded in

□

4.3.   Stability of the equilibrium points

The stability of the system (3.2) is examined here at each equilibrium point. The following system of equations must be solved in order to determine the equilibrium points:

Solving the aforementioned system, we obtain more than one equilibrium point. Let $(\nu^*, \rho^*, \psi^*, \mu^*)$ be one such equilibrium point. From Eq (4.11), we get

Solving (4.10) for $\psi^*$, we get

Solving (4.9) for $\rho^*$ and then substituting the value of $\psi^*$, it gives

From Eqs (4.12) and (4.14),

Substituting the values of $\rho^*, \psi^*, \mu^*$ in (4.8), we obtain an equation in terms of $\nu^*$.

where,

Here, two cases arise, because $\mu$ has two values:

ⅰ) If $\mu^*$ = 0, then

Since $B_1 < 0$, $B_2 > 0,$ $B_3 < 0$, if $aM < P^2$, by Descartes' rule of signs, Eq (4.16) has at most two real roots. If $aM > P^2$, then Eq (4.16) has at most three real roots.

ⅱ) If $\mu^* = \frac{l}{r},$ then

Since $B_1 < 0$, $B_2 > 0,$ $B_3 < 0$, if $(aM+\frac{l}{r}) < P^2$ then by Descartes' rule of signs, Eq (4.15) has at most two real roots. If $(aM+\frac{l}{r}) > P^2$, then Eq (4.15) has at most three real roots. Consequently, there is a coexistence point of equilibrium.

The Jacobian matrix of system (3.2) is

Theorem 4.3. The axial equilibrium point, $E_0 = (0, 0, 0)$, is always present, but it is unstable.

Proof. The eigenvalues of the Jacobian matrix $J$ at $E_0$ are

It is clear that $E_0$ is unstable. □

Theorem 4.4. When $P\neq0$ and $b = 0$, there are two unstable equilibrium points $E_1 = (\nu_1, \rho_1, 0, 0)$ = $E_1(\nu_1, \frac{-P}{\nu_1-1}, 0, 0)$ and $E_2 = (\nu_1, \rho_1, 0, \mu_1)$ = $E_2(\nu_1, \frac{-P}{\nu_1-1}, 0, \frac{l}{r}).$

Proof. (ⅰ) At $E_1$, the Jacobian matrix is

Eigenvalues of $J$ at $E_1$ are

$\lambda_{12} > 0$ indicates that the equilibrium point $E_1$ is unstable.

(ⅱ) At $E_2$, the Jacobian matrix is

The corresponding eigenvalues are

When b = 0, system (3.2) yields

If $a > 0$, then Eq (4.18) has a solution $\nu_1 < 1$, then $\lambda_{21} < 0.$ Since $\frac{8P^2}{(\nu_1-1)^2} > 0$ and $(\nu_1-1-a) > 0$, $\nu_1 > 1+a$ and it clearly shows that $Re(\lambda_{23}) < 0$ and $Re(\lambda_{24}) < 0$. Hence, $E_2$ is a stable equilibrium point.

□

4.4.   Conservative behavior

The system (3.2) is expressed in vector form as $^CD^\beta _{t_0} U = V(U)$ such that

where $V_1$–$V_4$ are defined in Eq (4.2). In this paper, we consider the values of the parameters as:

The vector region $f$ on $\mathbb{R}^4$ has a divergence ^[37] that can be expressed as

Consider any smooth-boundary region in $\mathbb{R}^4$ denoted by $\mathfrak{B}$. Let $\mathfrak{B}(t) = \mathfrak{P}_t(\mathfrak{B})$, where $\mathfrak{P}_t$ is the path of the vector field $f$. Let $\mathfrak{V}(t)$ represent $\mathfrak{B}$ hyper-volume. According to Liouville's theorem, it implies that

The $(div V)$ value can be substituted in (4.19) to get

When the linear differential inequation (4.20) is integrated, we get

We can conclude that $-1.94+2F < 0$ that is $F < 0.97$ is essential for the dissipativeness of the system (3.2). The volume $\mathfrak{V}(t)$ decreases exponentially to zero as $t\rightarrow \infty$, according to Eq (4.21). As a result, the novel chaotic system (3.2) is dissipative.

5.   Existence of chaos

5.1.   Lyapunove experiment

In this section, we explain the chaos that results from different fractional-order derivatives concerning time. $\beta = 1$ and $\beta = 0.95$ are chosen for this purpose. Values of the Lyapunov exponent for $\beta = 1$ and $\beta = 0.95$ are presented in Tables 1 and 2. Since at least two Lyapunov exponents are positive, for $\beta = 1$, $\beta = 0.95$, the system exhibits hyperchaotic behavior as shown in Figure 2. Hence, system (3.2) is a hyper-chaotic model. Authors of ^[7,38,39] have uncovered certain fractional-order hyperchaotic systems with a one-positive Lyapunov exponent.

5.2.   Poincaré map

An additional technique to describe the chaotic reaction is a Poincaré map or a return map. It forecasts the point on the trajectory where the current Poincaré section ^[40] crosses over. This mapping aims to achieve a minimum of two objectives:

(1) It helps to reveal any sequential order from the recurrent visit of paths to a phase space cross-section in graphical form.

(2) It aids in the prediction of the system's subsequent states in that Poincaré section.

In this section, we have considered two Poincaré sections: $U = \{\rho(t), \psi(t)\in \mathbb{R}^2| \nu(t) = -0.1\}$ and $V = \{\nu(t), \psi(t)\in \mathbb{R}^2| \rho(t) = 1\}.$

While the Poincaré sections across chaotic attractors of the system (3.2) are relatively complex, they typically exhibit some form of structure. We aim to observe the Poincaré section for changes concerning the growth rate $l$ and depletion rate $r$ of GHGs in the atmosphere. $r = 0.001, \, l = 0.3$ leads to Figure 3a, while $r = 0.0001, \, l = 0.0037$ leads to Figure 3b. Both figures represent the Poincaré sections parallel to the $\rho-\psi$ plane sectioned by $\nu = -0.1$. On the other hand, $r = 0.001, l = 0.0037$ lead to the Poincaré section of Figure 3c parallel to $\nu-\psi$ plane sectioned by $\rho = 1$. Moreover, a Poincaré section traversing a chaotic attractor typically consists of an infinite number of points. It is significant to remember that each surface in the phase plane's Poincaré section refers to a chaotic attractor. Thus, for $\beta = 1$, it is expected that system (3.2) will produce chaotic attractors. The outcomes of the modeling shown in Figure 4a are validated in this section via the Poincaré section parallel to $\rho-\psi$ and $\nu-\psi$ plane.

5.3.   Bifurcation

Quantitative bifurcation analysis is a tool to observe the existence of chaos in a system. Here, we examine the dynamics of the system (3.2) for the change of the values of the parameter $a$. The presence of chaos is shown in Figure 5 for $b = 11, \, P = 3, \, n = 0.25, \, l = 0.3, \, r = 0.0001$, and $M = 8$. The $a-\nu, \, a-\rho, \, a-\psi$ phase portraits for a range of parameter $a$ values can be seen in Figure 5. These figures indicate the existence of chaos.

6.   Design of sliding mode controller (SMC)

The fundamental principle of the SMC rule is to drive the system state trajectories to certain, predefined sliding surfaces where stability and other desirable system properties are present by use of a discontinuous control. This is accomplished by establishing the sliding surface, which facilitates the assessment of the fractional-order system's efficacy in sliding mode. The sliding mode controller is created after the sliding surface is determined. It is composed of two parts:

● We need to construct a sliding surface that depicts the intended system dynamics prior to creating a sliding mode controller.

● We need to create a switching control rule such that, at each point on the sliding surface, it will result in a sliding mode and drive any states that are outside of it to the surface in a finite amount of time.

A useful control method that guarantees both the occurrence of sliding motion and chaotic system control is the SMC law. Among its many fascinating characteristics are the resistance to parameters such as uncertainty and the desensitization to outside disruptions.

With the aid of an initial numerical analysis, we found that controlling the first equation is crucial to control the chaos within the system. The SMC law is a practical control strategy that ensures both chaotic system control and sliding motion. Resistance to constraints like uncertainty and desensitization to outside disruptions are only two of its many intriguing traits. Let us assume that $A(t)$ represents a control input that has to be applied in the following manner to the first state equation of the fractional-order system (3.2)

We propose a sliding surface $\lambda(U(t))$ as a function of the system states $U(t) = (\nu(t), \rho(t), \psi(t), \mu(t))$ as follows:

in which the function $\delta$(U(t)) is described as

such that it involves all the state variables. Here, the constant $s > 0$ is arbitrary.

The sliding surface and its derivative in a sliding mode must satisfy

From (6.2) and (6.4), we get $\lambda^\prime(U(t)) = ^CD^{\beta_1} \nu(t)+\delta(U(t)) = 0.$ Hence,

Hence, from the SMC theory, the first equation of systems (6.1) and (6.5) is applied to determine the corresponding control law

With respect to the sliding criterion, the reaching law can be selected as follows:

where

and the controller gain is represented as $W_r$.

Here, the total control law can be defined as follows:

Theorem 6.1. The fractional system (6.1) with the influence of the control law (6.7) is globally asymptotically stable if the controller gains $W_r < 0$.

Proof. Let us choose the Lyapunov candidate

Its derivative is computed by

Thus, the Lyapunov function X meets the conditions of the Lyapunov theorem, that is $X > 0$ and $\dot{X} < 0$. Hence, the controller system (6.1) with the impact of the control law (6.7) is globally asymptotically stable when $W_r < 0$. □

Theorem 6.2. The finite-time interval for fractional system (6.1) to reach the sliding surface is $[0, \tau_r]$ with $\tau_r = -\frac{1}{2W_r}(\lambda(U(0)))^2$.

Proof. For stability and finite-time convergence, the control gain $W_r$ must be chosen such that $\dot{X}(U(t))$ is negative definite. From Theorem 6.1, we have

where $W_r < 0$. We consider the time to reach the sliding surface to be $\tau_r$. Integrating both sides of Eq (6.9), from the initial time $t = 0$ to the time $\tau_r$, when $\lambda(U(t))$ tends to zero, we get

We have $X(U(\tau_r)) = 0$ and $X(U(0)) = \frac{1}{2}(\lambda(U(0)))^2$. This yields that

where $\lambda(U(0))$ is the value of the sliding surface at the initial state. By decreasing $W_r$, one can reduce the finite-time interval, ensuring the system reaches the sliding surface faster. □

Theorem 6.3. When there are uncertainities and an external disturbance, the fractional system (6.5) takes on the following form:

Where $\Delta h (\nu, \rho, \psi, \mu)$ and $f(t)$ are considered to be bounded such that $\Delta h(\nu, \rho, \psi, \mu) < g_1$, and $f(t) < g_2$, where $g_1, g_2 > 0$. Then, the above system (6.11) with the impact of SMC law (6.7) is globally asymptotically stable if $W_r < -(g_1+g_2)$.

Proof. We consider here the Lyapunov function as

We have

Hence, $\dot{X} < 0$ and $W_r < -(g_1+g_2).$ □

7.   Numerical simulations

Using the Adams-Bashforth-Moulton approach, which is compatible with FDEs, numerical simulations are performed ^[41]. We consider the parameter values at which the projected atmospheric model (3.2) becomes chaotic, namely $a = 0.24, \, b = 11, \, P = 3, \, n = 0.25, \, l = 0.3, \, r = 0.0001,$ and $M = 8.$ $(\nu_0, \, \rho_0, \, \psi_0, \, \mu_0) = (-0.1, \, -0.1, \, -0.1, \, 0.1)$ represents the initial approximation of state variable of the atmospheric model (3.2).

7.1.   Numerical validation of stability profile

The equilibrium points corresponding to the model (3.2) are $E_1(7.995, 0.002, 0.033, 0),$ $E_2(0.184, 0.509, 1.271, 0),$ $E_3(-0.163, 0.759, -1.175, 0),$ $E_4(3133, 0, 0.00008, 3000),$ $E_5(0.008-0.089i, 248.584+22.469i, 22.605-247.08i, 3000),$ and $E_6(0.008+0.089i, 248.584-22.469i, 22.605+247.08i, 3000).$

Here, equilibrium points $E_1$–$E_4$ are real, whereas $E_5$, and $E_6$ are complex conjugate. As atmospheric events are real phenomena, we shall consider only real equilibrium points.

$(1)$ Eigenvalues of the Jacobian matrix at $E_1$ are

From the Matignon criterion provided by ^[42],

As $|arg(\lambda_{1, 3})| < \frac{\beta\pi}{2}$ holds for any values of fractional derivative $\beta$, the equilibrium point $E_1$ is unstable.

(2) Eigenvalues of the Jacobian matrix at $E_2$ are

Since $|arg(\lambda_{2, j})| < \frac{\beta\pi}{2}, \, 3\le j\le4,$ for $0 < \beta\le1,$ the equilibrium point $E_2$ is unstable.

(3) At $E_3$, the eigenvalues of the Jacobian matrix are

Since $|arg(\lambda_{3, 4})|\le\frac{\beta\pi}{2}$, for $0 < \beta\le 1,$ the equilibrium point $E_3$ is always unstable.

(4) Eigenvalues of the Jacobian matrix at $E_4$ are

Since $|arg(\lambda_{4, i})| < \frac{\beta\pi}{2}, \, 1\le i\le2$, $0.945 < \beta\le1$, the system (3.2) is unstable at $E_4$ for $0.945 < \beta\le\; 1$.

7.2.   Dynamics of system (3.2)

We consider $\beta_1 = \beta_2 = \beta_3 = \beta_4 = \beta$. Figures 1 and 6 display the system trajectories of Eq (3.2) for commensurate orders $\beta = 1, \beta = 0.95, \beta = 0.85,$ and $\beta = 0.75.$ Figure 1 represents the dynamics of the system (3.1) illustrating the impact of a fixed input of $CO_2$. On the other hand, Figure 6 depicts the dynamics of the system (3.2) illustrating a continuous input of GHGs.

The Southern Hemisphere westerlies (SHWs) are crucial for regulating global climate and ocean circulation, but their future changes under varying levels of greenhouse gases (GHG) remain uncertain. According to this study, GHG-induced surface warming is unevenly distributed in the Southern Hemisphere under moderate to high GHG forcing, with the greatest warming occurring in sub-Antarctic coastal regions. Additionally, land areas on the continent are expected to warm more rapidly than the surrounding oceans. As a result, the equator-to-pole temperature gradient is likely to increase south of $60^0S$ latitude and decrease north of $60^0S$ latitude. This could lead to a strengthening of westerly wind speeds and a poleward shift of the SHWs ^[43]. The fourth equation of system (3.2) represents the changes in GHGs' intake in the atmosphere, highlighting the fact that, along with the growth of GHGs, there are depletions of GHGs due to their internal reactions. Chaotic behaviors are strongly observed in both Figures 1 and 6. However, in Figure 1, for various values of $\beta$, we observe the identical pattern of the strength of the westerly current, cosine, and sine waves. Whereas in Figure 6, an increasing strength of sine and cosine waves concerning time is noticeable. Also, the presence of continuous input of GHGs delays the amplitude of westerlies. This suggests that stronger westerlies result in more intense westerly currents compared to the new 4D atmospheric model (3.2) for the fractional order shown in Figure 6. 3D and 2D portfolios are shown in Figure 4. The increase in the strength of the westerlies can influence sea surface temperatures and marine ecosystems by affecting ocean circulation. Therefore, controlling these westerlies is critical.

A preliminary numerical investigation reveals that the installation of the control variable in the first equation of (3.2) will induce control in the underlying chaos. As a result, we obtained the resulting control system (6.1). The dynamics of the system's solution under the impact of SMC are projected in Figure 7, demonstrating that the resultant solution shows that chaotic behavior becomes controlled. The strength of the westerlies in the new 4D model increases as the fractional differential order decreases from $\beta = 1$ to $\beta = 0.75$. Moreover, we observe that when the fractional order decreases, the corresponding GHGs also decrease.

Figure 8a–d shows the dynamics of the non-commensurate fractional-order system (3.2) without controller, and Figure 8e–h shows the dynamics of the non-commensurate fractional-order system (6.1) with controller, for $\beta_1 = 0.85, \beta_2 = 0.95, \beta_3 = 0.9,$ and $\beta_4 = 1$.

We observed the dynamics of the controlled system (6.1) under the influence of uncertainties and external disturbances, namely $\Delta h(\nu, \rho, \psi, \mu) = 0.11 \sin(\nu^3+\rho^3)$ and $f(t) = 0.04 \sin(\pi t)$, respectively. Here, $|\Delta h(\nu, \rho, \psi, \mu)|\le g_1 = 0.11$ and $|f(t)|\le g_2 = 0.04$. By incorporating these parameters, we obtain the system (6.11). Figure 9 represent the dynamics of the non-commensurate ordered system (6.11) in the presence of assigned uncertainties and external disturbances. Westerlies, sine, and cosine waves remain in control even when external disturbances and uncertainties are applied. The system converges toward stability as pressure over the poles increases and the strong westerlies weaken. Figure 10 illustrates that even in the face of uncertainty, the controlled system (6.11) is efficient in controlling the chaos. In Figure 10, we show the power of the controller in system (6.11) for non-commensurate fractional derivatives $\beta_1 = 0.85, \beta_2 = 0.95, \beta_3 = 0.9,$ and $\beta_4 = 1$. The controller has forced the states of the system onto the sliding plane, where they remain for all following times. The simulation results demonstrate that in all three scenarios, the system responds to the derived control law (6.7). Put differently, the findings of the simulation demonstrate that the theoretical conclusions drawn are workable and effective for regulating the fractional-order atmospheric system.

7.3.   Validation of existence of Theorem 4.1

Using the following set of parameter values, we validated the existence of the solution for the projected model: $a = 0.24, \, b = 11, \, P = 3, \, n = 0.25, \, l = 0.3, \, r = 0.0001, \, M = 8.$ Furthermore, taking the guidance from the bounds shown in Figure 6, we considered $\eta_1 = 4$, $\eta_2 = 6.5$, $\eta_3 = 6.5$, and $\eta_4 = 170.$

(1) For a maximum $\Phi_1$, we examined the values of fixed parameters. Thus, we obtained $\Phi_1 = 156.24, \Phi_2 = 57.2, \Phi_3 = 57.2, \Phi_4 = 0.584.$ In this case, the strength of the westerly current is maximized due to the combined impact of sine and cosine waves; these observations are shown in Figures 4a, 6a, 6e, and 6i.

(2) To get the $\Phi_4$ maximum in comparison to others, we increased the depletion rate to $r = 0.5,$ maintaining the same values for the other parameters. Consequently, we obtained $\Phi_1 = 156.24, \Phi_2 = 57.2, \Phi_3 = 57.2, \Phi_4 = 170.55.$ When the depletion rate of GHGs increases, the exponential growth of GHGs (Figure 6m–p) converts to steady growth (Figure 11d), and this leads to depleting the amplitude of the strength of the westerly current. These observations are shown in Figure 11.

8.   Conclusions

In summary, our research on the fractional four-dimensional chaotic atmospheric model has yielded significant insights into the existence and uniqueness of solutions, the system's dissipative nature, and methods for controlling chaos. We calculated the Lyapunov exponents, which are critical for understanding the system's stability and chaotic behavior. Poincaré map and bifurcation analysis strengthened the existence of chaos in the projected model. To address the chaos in the atmospheric system, we derived a sliding-mode controller. This control mechanism is designed to globally stabilize the system, ensuring it can operate under the influence of external disturbances. The numerical simulations illustrate that the sliding-mode controller successfully stabilizes the chaotic system, thereby proving its practical applicability. The control law we proposed effectively mitigates chaos for both commensurate and non-commensurate orders of the system. All theoretical findings are well-supported by numerical simulations and graphical representations, enhancing our understanding of complex atmospheric models.

Author contributions

Conceptualization: A.J.R., C.B., R.G.; Methodology: A.J.R., R.G., Z.M.Y.; Software: C.B., S.E.; Validation: C.B., S.E., Z.M.Y.; Formal analysis: C.B., S.E., Z.M.Y.; Investigation: R.G., S.E., Z.M.Y.; Writing-original draft preparation: A.J.R., C.B., S.E.; Writing-review and editing: C.B., R.G., Z.M.Y.; Supervision: C.B., R.G.; Project administration: R.G. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

The authors extend their appreciation to Prince Sattam bin Abdulaziz University, Saudi Arabia for funding this research work through the project number (2024/01/922572).

Conflict of interest

The authors declare no conflicts of interest.