Modeling and simulations for the mitigation of atmospheric carbon dioxide through forest management programs

1.   Introduction

Carbon dioxide concentrations in the environment have risen rapidly during the last several years. Global warming is primarily caused by excessive carbon dioxide levels, which have detrimental effects on air quality and human health. Climate change has several harmful consequences for people and ecosystems, such as melting glaciers and icecaps, tidal waves, floods, and siltation; increased chances of annihilation of endangered species of animals and plants; alterations in rainfall patterns modifying the global supply of water and food; a rise in food-borne, water-borne, and vector-borne diseases; and a bump in high-temperature diseases ^[1,2]. The atmospheric $CO_2$ concentration surpassed 418 ppm in 2022, nearly 49% higher than the value of 280 ppm at the start of the industrial age ^[3]. One of several primary human-caused origins of carbon dioxide emissions is deforestation. Deforestation has resulted in the loss of almost 420 million ha of forests globally since 1990, with an average pace of 10 million ha year$^{-1}$ between 2015 and 2020 ^[4].

Many nations worldwide have implemented reforestation and afforestation plans to compensate for the loss of forests driven by deforestation. It was predicted that 27 million acres of new trees would be planted in 2015, an increase of 1.57% yearly ^[5]. Large-scale plantations are necessary to reduce elevated $CO_2$ concentration in the air; however, reforestation is not possible on the needed scale due to a variety of demographic and economic factors ^[6]. Nations in this sort of situation are implementing sustainable forestry strategies that enhance forest biomass and minimize deforestation rates, helping to reduce $CO_2$ emissions in the atmosphere. Since being the first tropical nation to successfully control deforestation, Costa Rica increased its forest area from 24.4% in 1985 to more than 50% by 2011. This was accomplished via the use of forest management strategies ^[7].

Innovative methods of managing forests are being developed due to global concerns over forest degradation. One of several growth strategies is to enhance the amount of high biomass in the forest and diminish the atmospheric $CO_2$ level using genetically modified plants with longer roots, rapid growth rates, and biomass production ^[8,9]. Genetically manipulated eucalyptus plants have been shown to grow quicker and soak up more $CO_2$, leading to increased forest biomass output per unit area in Brazil ^[10]. Another technique that is frequently employed in modern times to boost productivity and forest cover is agroforestry. To increase production and the sustainability of the ecosystem, agroforestry defines the system of land use that incorporates trees into agricultural landscapes and farms ^[11]. Agroforestry is a crucial component of the broader regional plan for managing forests and reducing climate crises. Indian agroforestry can sequester 25 tonnes of carbon on average per hectare, according to estimates ^[12,13]. Forests provide a significant source of income, firewood for heating and cooking, and other necessities for people in rural communities. The population's reliance on forestry resources is frequently brought on by the absence of alternatives or the incapacity of the inhabitants to pay for them ^[14]. Some initiatives offering rural residents government subsidies for their way of life, such as biogas and fuel-efficient stoves, as well as fiber and tin, successfully reduce deforestation rates by promoting alternatives to forest resources ^[15].

Multiple approaches to mitigating climate change have been designed in the past few years, with a particular emphasis on reducing atmospheric $CO_2$ emissions from forest degradation and deforestation. Reducing emissions from deforestation and degradation plus (REDD+) system, initiated by the United Nations Framework Convention on Climate Change (UNFCCC) Conference of the Parties, is one of the primary worldwide initiatives in this field. With sustainable forest management, the REDD+ mechanism seeks to direct and financially support actions that minimize deforestation and boost carbon stocks already present in the forest ^[16,17]. It has been determined that a REDD+ program implemented at the state scale can considerably lower carbon emissions and forest cover loss. According to research, the Norway-Guyana REDD+ initiative prevented 12.8 million metric tons of $CO_2$ emissions and reduced tree cover loss by 35% between 2010 and 2015 ^[18]. Therefore, the actions involved in forest management play a vital role in lowering the pace of deforestation and increasing the biomass of the forest.

Recently, several mathematical models have been put out to investigate the behavior of different aspects of the human population ^{[19,20,21,22,23,24,25]}. The authors in ^[26] investigated the stability criteria of a mathematical model for the equilibrium points. It was revealed that the equilibrium between biomass and carbon dioxide becomes unstable, leading to a rise in atmospheric $CO_2$ at highly elevated deforestation rates. The authors in ^[27] used a mathematical model to determine how the population, biomass, and atmospheric $CO_2$ interacted. They concluded that the rate of deforestation has adversely affected the system's dynamics. Research on how different plant capacities to absorb $CO_2$ affect atmospheric $CO_2$ level is given in ^[28]. To investigate the impact of population on the dynamics of $CO_2$, the authors in ^[29] designed a dynamic model. According to their research, the pressure of population growth increases the rate of forest biomass loss, which raises the equilibrium level of $CO_2$. To preserve forest biomass, several researchers have proposed nonlinear mathematical models ^[30,31,32]. Depending on these studies, it is possible to preserve forest biomass using technical means such as the cultivation of genetically modified plants and by offering financial incentives to the population, which would eventually result in a decline in the pace of deforestation. These investigations, however, need to examine the consequences of forest management on $CO_2$ concentration. In this manuscript, we formulate a novel fractional mathematical model to explore the influence of forest management strategies on the emission of $CO_2$. According to our assumptions, sustainable management strategies aim to reduce deforestation by offering financial incentives and encouraging individuals to use different resources. These efforts aim to lower deforestation rates and enhance forest biomass via afforestation and the planting of genetically modified plants. The models stated above rely on classical derivatives and frequently fail to adequately represent the complexities of systems that exhibit unusual spreading and non-exponential response.

Modern research has highlighted the potential of using fractional differential equations to represent a variety of events across various disciplines, especially epidemiology ^{[33,34,35,36,37]}. The calculus of arbitrary-order derivatives and integrals is studied in the domain of mathematical analysis known as fractional calculus. The vital feature of utilizing fractional derivatives to describe mathematical models is that they have the quality of being nonlocal. The $m$-th derivative of a function $\Phi(\epsilon)$ at $\epsilon$ is just a local characteristic when $m$ is an integer. The fractional derivative of $\Phi(\epsilon)$ at $\epsilon = a$ depends on all values of $\Phi(\epsilon)$, including those farthest from $a$. Fractional-order derivatives incorporate derivatives and integrals of arbitrary orders; the current and past states determine the future predicament. Memory's inherited qualities and potency are key aspects across many biological systems. A dynamic model can assist in expressing these aspects using fractional-order derivatives in modeling ^[38].

The most famous operators are the Caputo fractional derivatives and integrals, which are traditionally employed for simulating several real-life problems. Additionally, there is a close connection between the Riemann-Lioville derivative operator and the Caputo fractional derivative. These operators could yield better findings when used to examine the structure of models. However, the singularity of the kernels of these operators ^[39,40] is the primary problem. To better comprehend the dynamics of models, researchers have emphasized the need for fractional operators with nonsingular kernels. Many researchers were successful in providing fractional operators with nonsingular kernels to this point, such as the Caputo-Fabrizio fractional derivative ^[41]. Recently, Atangana and Baleanu introduced a brand-new derivative with a single-parameter Mittag-Leffer (ML) kernel, which is the Atangana-Baleanu (ABC) operator ^[42]. One of its major achievements is that this operator has a non-local and non-singular kernel. It is advantageous for professionals who work in the dynamic simulation of real-life problems. Various conventional and fractional derivatives were taken into consideration when performing numerical computations in ^{[43,44,45,46,47]}. To approximate the ABC fractional derivative, the authors carried out an examination and numerical simulation of a model ^[45]. They applied a two-step family of Adams-Bash-Forth methodology.

To overcome these limitations, a unique use of the ABC fractional derivative operator is presented in this study. In contrast to earlier research, which mostly dealt with traditional fractional derivatives, this work makes use of the unique benefits of the ABC operator to improve complex system modeling. This paper's key contribution is its demonstration of how fractional derivatives of ABC can enhance the stability and precision of solutions to fractional differential equations. We demonstrate the superiority of the ABC operator over conventional techniques in terms of solution fidelity and computing economy through several in-depth evaluations and simulations. The suggested method has been widely used in many real-world problems due to its non singular and non-local properties ^[48,49], and has already been compared and proved to be superior to the method in the literature. To the best of our knowledge, this method has not yet been applied to problems related to climate change and sustainability. Thus, this study introduces a new methodology that not only links theoretical developments with real-world applications but also establishes new standards for precision and effectiveness in complex system modeling. Here is a list of the contributions of this research work:

● A novel mathematical fractional model for the dynamics of $CO_2$ concentration in the atmosphere is addressed and analyzed using the Atangana-Baleanu derivative.

● Non negativity of the solutions is proven, and the conditions for the existence and uniqueness of the solution of the considered system are established.

● The proposed fractional model is solved numerically using the Atangana-Toufik scheme.

● The effects of the fractional ABC derivative on the dynamics of the proposed model are discussed.

1.1.   Preliminary definitions

In this section, we recall some fundamental definitions from fractional calculus and a well-known theorem that will be used in the following sections.

Definition 1.1. Suppose that $\Lambda$ is an open subset of $\mathbb{R}$, and $\omega \in [1, \infty)$, the Sobolev space $H^\omega(\Lambda)$ is stated ^[35,45] as:

Definition 1.2. The ABC fractional operator in the Caputo sense ^[35,45] is defined for a differentiable function $\wp : [a, b] \rightarrow \mathbb{N}$, which is defined on $[a, b]$ such that $\wp \in H^1(a, b)$, $b > a$, and $\delta \in (0, 1]$ as:

where $E_\delta$ is a ML function of one parameter, and given by

where $\aleph(\delta) > 0$, is a normalized function, and $\aleph(0) = \aleph(1) = 1$.

Definition 1.3. ^[45] The associated ABC fractional integral with a non-local kernel is expressed as

Theorem 1.4. ^[35] The fractional-order differential equation

has the following unique solution

The remaining paper is organized as follows: Section 2 formulates the fractional model and its associated characteristics. In Section 3, the equilibrium states of the proposed fractional model are calculated. Section 4 examines the stability of the equilibrium points of the suggested model. Numerical analysis and simulations are carried out in Section 5. Finally, Section 6 is devoted to the conclusion.

2.   Model formulation

The growing study of climate change demands a thorough comprehension of greenhouse gas dynamics. The $CO_2$ gas is one of the main greenhouse gases causing global warming. Precisely estimating the gas' concentration and how it impacts the environment is essential in order to forecast future climate circumstances and develop practical mitigation plans. The complex mechanisms controlling $CO_2$ emissions are well represented by the integer-order model, which offers a stable framework. By integrating these processes into a logical system of differential equations, the model allows accurate simulations of time-varying $CO_2$ levels. For researchers and public officials to evaluate the effects of different emission reduction programs and to create well-informed environmental policies, this forecasting endurance is crucial.

In this study, we consider a geographical region where a growing human population affects the forest cover. Sustainable forest strategies are implemented to decrease deforestation and enhance forest biomass. As it acts as one of the fundamental $CO_2$ sinks, the conservation, and depletion of forest biomass will change the concentration of $CO_2$. To describe the model, we assumed four different state variables: the concentration of $CO_2$ in the atmosphere $\mathcal{C}({t})$ at any time $t$, the human population $\mathcal{N}(t)$ at any time $t$, the biomass of the forest $\mathcal{B}(t)$ at any time $t$, and the number of forest management programs $\mathcal{P}(t)$ at any time $t$.

The differential equation below describes how the concentration of carbon dioxide in the atmosphere changes over time:

here, the natural emission rate of $CO_2$ from processes including respiration, biological decomposition, and volcanic activity, all of which continuously contribute $CO_2$ to the atmosphere, is represented by the term $Q$. The concept of $\lambda\mathcal{N}$, which stands for emission rate coefficient and scales with the human population or activities that cause $CO_2$ emissions (such as burning fossil fuels and industrial operations), describes emissions of $CO_2$ that are caused by humans. The term $-\alpha\mathcal{C}$ refers to the amount of $CO_2$ absorbed by natural sinks other than forest biomass, such as soil and seas; this absorption is proportionate to the amount of $CO_2$, or $\mathcal{C}$, that is currently present in the atmosphere. Last, the absorption of $CO_2$ by forest biomass is denoted by the term $-\lambda_1\mathcal{B}\mathcal{C}$, where $B$ is the quantity of forest biomass and $\lambda_1$ is the rate of uptake coefficient. This uptake depends on both the amount of $CO_2$ and the quantity of forest biomass.

Given that population expansion is dependent on forest resources and causes deforestation it is assumed that both the population and biomass of the forest would expand logistically. A rise in biomass, on the other hand, promotes population expansion. The efficiency of forest management programs decreases with intensity, while they do prevent deforestation and increase forest biomass. A decrease in population is also a result of the negative climatic effects brought on by elevated atmospheric $CO_2$. With the consequences of deforestation, the management of forests, and $CO_2$-induced climate change taken into account, these dynamics are reflected in the following differential equations that describe population and forest biomass as:

We assumed that the rate at which forest management initiatives are carried out is proportionate to the gap between the carrying capacity of the forest and its present biomass density. Let $\nu$ and $\nu_0$ denote the implementation and declination rates of the forest management initiatives, respectively. The following is a description of the dynamics of the forest management programs:

Therefore, the following model incorporates the dynamic nature of the problem:

with the following initial conditions:

It is assumed that all the involved parameters and initial data of the model are non-negative. The model parameters used in this model are defined and interpreted in Table 1. Generally, fractional models of real-life phenomena are known to express the memory effect quite effectively; the system (2.6) is considered using the ABC fractional derivative. By applying the ABC fractional derivative of order $\delta \in (0, 1]$ of the ML kernel, the system is given as:

and the following short form can be used to express it as:

with the initial condition

where $\wp : [0, +\infty)\rightarrow \mathbb{R}^4$ and $G : \mathbb{R}^4 \rightarrow \mathbb{R}^4$ are vector-valued functions such that

and

respectively. Clearly, $G_j;$ $j = 1, 2, 3, 4,$ are functions of $t$, $\mathcal{C}$, $\mathcal{N}$, $\mathcal{B}$, and $\mathcal{P}$.

2.1.   Fundamentals of the model

In this part, we use nonnegative initial information to define the range of solutions to the non-linear system. Our fundamental aim is to present that the observed feasible region in $\mathbb{R}^{4}_+$ is positively invariant about the suggested model.

2.1.1.   Invariant region

Theorem 2.1. ^[46] The set

is positively invariant for the fractional model (2.6), where $\mathcal{C}_m = \frac{Q + \lambda \mathcal{N}_m}{\alpha}$, $\mathcal{N}_m = L + \frac{\xi LM}{s}$, and $\mathcal{P}_m = \frac{\nu M}{\nu_0}$.

2.2.   Non-negativity and boundedness of the solutions

Under this section, we explore the physical characteristics of solutions, i.e., non negativity and boundedness of solutions to the presented model. Here, we establish the necessary and sufficient criteria for the existence of positive solutions. Our objective is to demonstrate the qualitative behavior of solutions. Therefore, we suppose that for any initial conditions for the system, $\mathcal{C}(0) = \mathcal{C}_0$, $\mathcal{N}(0) = \mathcal{N}_0$, $\mathcal{B}(0) = \mathcal{B}_0$, and $\mathcal{P}(0) = \mathcal{P}_0$ $\in \Omega$, then there is a unique solution for any time $t$. Let us now define positivity criteria.

Theorem 2.2. ^[50] Assume that the initial condition

then, the solutions exist, and they are all positive for all $t\geq 0$.

Proof. Let's begin with the $\mathcal{C}(t)$ class. We define the following norm:

First equation of the system (2.6) is given by

Then, we get

Second equation of the system (2.6) is given by

Then, we have

Third equation of the system (2.6) is given as

Then, we get

Fourth equation of the system (2.6) is written as

Then, we have

Thus the suggested model has non-negative solutions $\forall t\geq0$. □

Theorem 2.3. ^[51] Along with the initial conditions, the solution of the fractional model considered in (2.6) is unique and bounded in $\mathbb{R}_+^4$.

Proof. We obtained the following equations:

If $\left({\mathcal{C}(0), \mathcal{N}(0), \mathcal{B}(0), \mathcal{P}(0)}\right) \in \mathbb{R}^4_+$, it is practically unattainable to move out of $\mathbb{R}^4_+$. Since the domain is positively invariant, any non negative vector field goes into $\mathbb{R}^4_+$. For the suggested system, this guarantees a positive and bounded solution. □

2.3.   Existence and uniqueness of solution

In this part, we employ theorems to show the existence and uniqueness of the solution for the suggested model (2.6).

Theorem 2.4. ^[47] Suppose that we have eight positive constants $\hat{\rho_1},$ $\hat{\rho_2},$ $\hat{\rho_3},$ $\hat{\rho_4}$, and $\bar{\rho_1},$ $\bar{\rho_2}$, $\bar{\rho_3}$, $\bar{\rho_4}$, then the following two conditions hold:

(1)

(2)

Proof. The system (2.6) has a unique solution if the stated requirements are fulfilled. Let us begin with the first equation of the system $G_1(\mathcal{C}, t)$. Then, we will try to prove the first condition for the following equation:

Define the following norm as:

Let us consider $\mathcal{C}$, $\mathcal{C}_1$ $\in \mathbb{R}^2$ and $t\in [0, T]$,

where

For $\mathcal{N}$, $\mathcal{N}_1$ $\in \mathbb{R}^2$ and $t\in [0, T]$,

where

For $\mathcal{B}$, $\mathcal{B}_1$ $\in \mathbb{R}^2$ and $t\in [0, T]$,

where

For $\mathcal{P}$, $\mathcal{P}_1$ $\in \mathbb{R}^2$ and $t\in[0, T]$,

where

Finally, we able to establish the given condition (1).

Now, we will prove the second criterion for the proposed system. $\forall(\mathcal{C}, t) \in \mathbb{R}^2 \times [t_0, T]$, then we will show that

where

and with under condition

Now $\forall(\mathcal{N}, t) \in \mathbb{R}^2 \times [t_0, T]$, then we will show that

where

and with under condition

Now $\forall(\mathcal{B}, t) \in \mathbb{R}^2 \times [t_0, T]$, then we will show that

where

and with under condition

$\forall(\mathcal{P}, t) \in \mathbb{R}^2 \times [t_0, T]$, then we will show that

where

and with under condition

which completes the proof. □

3.   Equilibrium points of the system

To determine the equilibrium points of the system (2.6), assume that the rate of change with regard to time is zero, then we obtain

By solving the above system, we get the following four steady states denoted by $E_1$, $E_2$, $E_3$, and $E_\ast$:

(1) $E_1\left(\frac{Q}{\alpha}, 0, 0, \frac{\nu M}{\nu_0}\right)$ always exists.

(2) $E_2\left(\frac{Q}{\alpha+\lambda_1 M}, 0, M, 0\right)$ always exists.

(3) $E_3\left(\mathcal{C}_3, \mathcal{N}_3, 0, \mathcal{P}_3\right)$ exists, if the following condition is fulfilled:

where $\mathcal{C}_3 = \frac{s(Q+\lambda L)}{s\alpha+\theta \lambda L}$, $\mathcal{N}_3 = \frac{L(s\alpha-\theta Q)}{s\alpha+\theta \lambda L}$, $\mathcal{P}_3 = \frac{\nu M}{\nu_0}$.

(4) $E_\ast \left(\mathcal{C}_\ast, \mathcal{N}_\ast, \mathcal{B}_\ast, \mathcal{P}_\ast \right)$ exists, if the following conditions are fulfilled:

4.   Stability analysis

A theoretical discussion on the local and global asymptotic stabilities at the equilibria of a fractional model (2.6) is presented in this section.

4.1.   Local asymptotic stability

Theorem 4.1. 1) The equilibrium point $E_1$ is always unstable.

2) The equilibrium point $E_2$ is unstable whenever $E_\ast$ exists.

3) The equilibrium point $E_3$ is unstable whenever $E_\ast$ exists.

4) The equilibrium point $E_\ast$ is locally asymptotically stable if and only if the given condition is satisfied:

Proof. For the proposed fractional system (2.6), the Jacobian matrix $\mathbf{\hat{J}}$ is evaluated as:

Let $\mathbf{\hat{J}_{E_1}}$, $\mathbf{\hat{J}_{E_2}}$, $\mathbf{\hat{J}_{E_3}}$, and $\mathbf{\hat{J}_{E_\ast}}$ represent the Jacobian matrices of the considered fractional model at $E_1$, $E_2$, $E_3$, and $E_\ast$, respectively. Then

1) The characteristics equation of the $\mathbf{\hat{J}_{E_1}}$ is given by $|\mathbf{\hat{J}_{E_1}} - \hat{\lambda} I| = 0$, where $I$ is the unit matrix. This characteristics equation gives the eigenvalues of $\mathbf{\hat{J}_{E_1}}$, which are

From Eq (4.2), $\hat{\lambda}_3 > 0$, as all the involved parameters are positive; therefore, ${E_1}$ is always unstable.

2) The characteristics equation of the $\mathbf{\hat{J}_{E_2}}$ is given by $|\mathbf{\hat{J}_{E_2}} - \hat{\lambda} I| = 0$. This characteristics equation gives the eigenvalues of $\mathbf{\hat{J}_{E_2}}$, which are

and $\hat{\lambda}_3$, $\hat{\lambda}_4$ are complex with negative real part.

Since, $\hat{\lambda}_1 < 0$, as $\alpha$, $\lambda_1$, and $M$ are positive. From Eq (4.3), $\hat{\lambda}_2 = s + \xi M - \frac{\theta Q}{\alpha + \lambda_1 M} > 0$, if condition (3.2) is satisfied. Hence, ${E_2}$ is unstable whenever ${E_\ast}$ exists.

3) The characteristics equation of the $\mathbf{\hat{J}_{E_3}}$ is given by $|\mathbf{\hat{J}_{E_3}} - \hat{\lambda} I| = 0$. Then, $\mathbf{\hat{J}_{E_3}}$ has the following eigenvalues:

and $\hat{\lambda}_3$, $\hat{\lambda}_4$ are roots of the following equation

which are complex with negative real part.

From Eq (4.4), $\hat{\lambda}_2 = \mu - \mathcal{N}\left(\phi - \frac{\phi_1 \mathcal{P}}{k_1 + \mathcal{P}}\right) + \frac{\eta_1 \mathcal{P}}{l_1 + \mathcal{P}} > 0$, if condition (3.1) is satisfied. Therefore, the considered system at ${E_3}$ is unstable whenever ${E_\ast}$ exists.

4) The characteristics polynomial of the $\mathbf{\hat{J}_{E_\ast}}$ is given by $|\mathbf{\hat{J}_{E_\ast}} - \hat{\lambda} I| = 0$, and can be written as:

where

On inspection, it is seen that all $\Delta_j{'s} > 0, j = 1, 2, 3, 4$. Following the Routh-Hurwitz rule, if the condition (4.1) is achieved, all the solutions of Eq (4.5) will be found in the negative half of the plane, which gives that all of the eigenvalues are negative. Therefore, the equilibrium point ${E_\ast}$ is locally asymptotically stable.

□

4.2.   Global asymptotic stability

Theorem 4.2. If the equilibrium $E_\ast$ of the presented fractional model (2.6) exists, then it is globally asymptotically stable in $\Omega$, assuming the following conditions are fulfilled:

Proof. To show the global asymptotic stability of the $E_\ast$, we define the following Lyapunov function:

where $\sigma_1$, $\sigma_2$, and $\sigma_3$ are positive constants that will be opted later.

Applying the ABC derivative with respect to $t$ on Eq (4.8), we have

Choosing $\sigma_1 = \frac{\lambda}{\theta}$, $\sigma_2 = \sigma_1 \frac{\xi}{\left(\phi - \frac{\phi_1 \mathcal{P}_\ast}{k_1 + \mathcal{P}_\ast}\right)} = \frac{\xi \lambda}{\theta \left(\phi - \frac{\phi_1 \mathcal{P}_\ast}{k_1 + \mathcal{P}_\ast}\right)}$, we have

Therefore, $_0^{ABC}D_{t}^{\delta} \mathbf{U_\ast}$ is negative definite in $\Omega$ if the following conditions are hold:

From inequalities (4.9)–(4.12), and stated condition (4.7) holds, then we can choose $\sigma_3 > 0$. Therefore, the conditions (4.6) and (4.7) are achieved, and $_0^{ABC}D_{t}^{\delta}\mathbf{U_\ast}$ is negative definite, which implies that $E_\ast$ is globally asymptotically stable in $\Omega$. □

5.   Numerical analysis

In this segment, we apply a well-known numerical scheme known as the Atangana-Toufik scheme ^{[35,43,48,49]} to find the numerical results of the proposed model for different values of $\delta$. Parameter values used in the simulation are provided in Table 2.

5.1.   Atangana-Toufik scheme

For the solution of a fractional model, we construct a numerical scheme using the Atangana-Toufik scheme. With this goal, we will consider that we get the following results from the system (2.6):

For a given point $t = t_{m+1}$, $m = 0, 1, 2 \ldots,$ then Eqs (5.1)–(5.4) are written as

Also, we have

The function $G_j(\tau, \wp(\tau))$ can be approximated applying two-step Lagrange polynomial interpolation in the interval $[t_k, t_{k+1}]$ as follows:

Therefore, the previous approximation can be added to Eqs (5.5)–(5.8) to give

where

and

Substituting the above integrals into Eqs (5.5)–(5.8), we have the following iterative schemes for the model equations:

Figures 1–4 depict the behavior of all the compartments within the system when the fractional parameter is changed from 0.80 to 1.00. This study sheds light on the advantages of applying the ABC fractional derivative operator by demonstrating the fractional parameter's substantial influence on the dynamics of the system.

5.2.   Simulations

Numerical simulations are carried out for the parameter values listed in Table 2 to examine the consequences of forest management policies and deforestation on the dynamics of atmospheric $CO_2$. Thus, we set the natural $CO_2$ emission rate, $Q = 1$ ppm month$^{-1}$. Since the population's carrying capacity is $L = 1000$ persons, we set the population growth rate coefficient induced by forest biomass, $\xi$, to be $0.0000002$ (ton month)$^{-1}$. For the simulations, the following initial conditions are considered: $\mathcal{C}(0) = 300$, $\mathcal{N}(0) = 500$, $\mathcal{B}(0) = 500$, and $\mathcal{P}(0) = 100$ by setting $\phi = 0.0003$ (person month)$^{-1}$, $\mu = 0.2$ month$^{-1}$, and $\theta = 0.00001$ (ppm month)$^{-1}$.

Figure 1 shows that there is a discernible shift in the response of $CO_2$ concentration when the fractional parameter increases from 0.80 to 1.00. The behavior of atmospheric carbon dioxide is reduced with smaller fractional levels. It is evident from Figure 2 that higher values in the fractional parameter lead to a rise in the human population. Figure 3 illustrates an inverse connection between a fractional parameter and the human population. There is a noticeable tendency at lower fractional orders that progressively decreases as the parameter rises. This implies that more complicated, real-world dynamics that are not as well characterized by integer-order derivatives are captured by fractional-order derivatives. Figure 4 illustrates how fractional parameters affect forest management initiatives and exhibit behavior that is comparable to Figure 1. The behavior is more consistent with integer-order derivatives as the parameter gets closer to 1.00, indicating that it is converging to the right steady-state ^[46]. As seen in these diagrams, the fractional derivative has a significant impact on the concentration of $CO_2$ throughout time. The simulations performed using fractional derivatives produce the reported gap due to the memory effect, which is missing in the classic integer-order derivative. Future $CO_2$ dynamics will be influenced by memory effects on previous data, making it easier to control them.

Figures 5 and 6 show how deforestation and forest management strategies affect atmospheric $CO_2$ dynamics. These statistics emphasize how important it is to implement effective mitigation techniques and switch to low-carbon energy sources to lower atmospheric $CO_2$ levels.

The relationship between deforestation activities and rising $CO_2$ concentrations in the atmosphere is seen in Figure 5. The findings show that regions with considerable deforestation have greater $CO_2$ concentrations. This is explained by the fact that less forest cover means less carbon being sequestered. As essential carbon sinks, forests take up carbon dioxide from the atmosphere. Deforestation, which removes trees from the forest, causes the atmosphere to re-release stored carbon, intensifying the greenhouse effect and accelerating global warming. The results of using different forest management strategies targeted at reducing $CO_2$ emissions are shown in Figure 6. The picture shows how $CO_2$ levels gradually drop over time in areas with strict forest management strategies. Reforestation, afforestation, and sustainable logging techniques are examples of effective policies that guarantee the maintenance and improvement of forest carbon sinks. By improving tree's capacity to store carbon, these methods can lower the atmospheric concentration of $CO_2$ globally.

6.   Conclusions

The present study developed a fractional-order model to analyze the dynamics of novel atmospheric $CO_2$ gas through forest management programs. The Atangana-Baleanu derivative operator and generalized Mittag-Leffler function were employed to develop this novel fractional model. Carbon dioxide, the leading greenhouse gas, was mitigated in the atmosphere using sustainable forest management. Our study aimed to examine the memory effect and the impact of forest management initiatives on stabilizing elevated $CO_2$ levels in the atmosphere. The existence and uniqueness of the solution for the proposed model were proven using theorems. Steady states were computed and investigated for stability analysis. The equilibrium point $E_\ast$ of the model is locally and globally asymptotically stable in the region $\Omega$ using the Lyapunov function theory. The Atangana-Toufik scheme was utilized to find the approximate solution to the proposed model. The results obtained are represented graphically and are consistent with the system. Simulations were conducted for different values of fractional order $\delta$, revealing that the fractional order significantly impacts the concentration of $CO_2$. Analysis of deforestation and forest management programs indicates, that to enhance forest cover and reduce atmospheric $CO_2$ levels, it is more advantageous to implement forest management strategies that reduce the rate of deforestation and increase the growth rate of forest biomass at low implementation costs. A constraint of this research is the assumption of uniformity in forest management methodologies, which does not accurately represent actual diversity in the field. In the future, we will take into account different forest kinds and management techniques. We will also add climate change aspects to the model, which might improve its forecasting power.

Author contributions

Muhammad Bilal Riaz: Conceptualization, methdology, validation, formal analysis, investigation, writing—review and editing, writing—review and editing; Nauman Raza: Software, formal analysis, visualization; Jan Martinovic: Software, validation, writing—review and editing, project administration; Abu Bakar: Software, investigation, writing—original draft preparation; Osman Tunç: Data curation, writing—review and editing, visualization. The authors declare that this study was accomplished in collaboration with the same responsibility. All the authors read and approved the submitted Manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This article has been produced with the financial support of the European Union under the REFRESH – Research Excellence For Region Sustainability and High-tech Industries project number CZ$.10.03.01/00/22\_003/0000048$ via the Operational Programme Just Transition.

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this manuscript.