Integrated Modelling of Carbon Emissions under Population Growth, Forest Management Policies, and Tree Plantation Dynamics

The dynamics of carbon dioxide (CO$_2$) in the atmosphere are partially due to a complex interplay of natural processes and human activities. As global environmental concerns increase, it is important to understand how various factors such as human population growth, land use changes, forest area, and tree plantations mark the concentration of carbon dioxide and, in turn, contribute to climate change. This model provides a framework for analyzing these interrelated components and their impact on carbon dynamics in the environment [1], [2].

Carbon concentration in the atmosphere, often regarded as a major driver of global warming and climate change, is largely molded by both natural and anthropogenic factors [3].

In contrast to the release of carbon dioxide, natural processes such as photosynthesis in plants, carbon sequestration in soil, and the absorption of CO$_2$ by oceans and forests play a significant role in balancing atmospheric carbon levels. Forests, in particular, are crucial in mitigating the rise of CO$_2$, as they absorb large quantities of carbon over the process of photosynthesis. However, human activities such as deforestation and land conversion often disrupt this balance, releasing stored carbon into the atmosphere and exacerbating the greenhouse effect [4].

Globally, land use and land cover change can increase or decrease carbon sequestration potential (IPCC). Forest conversion to other land uses can reduce carbon significantly; however, plantation establishment can boost the rate of carbon sequestration. Agricultural expansion and increasing demand for wood as fuel and timber, associated with population growth, are reducing the carbon sequestration potential. Due to the woody and long life of trees, the Kyoto Protocol of the Intergovernmental Penal on Climate Change (IPCC) declares forest plantation as a cost-effective approach for carbon absorption. Selection of the best stock for plantation can enhance the carbon absorption because of their fast growth and good health. Plantation consisting of the best stock will have the same investment cost, but will absorb more carbon [5], [6], [7], [8], [9], [10].

The human population is a crucial factor in the model, as the size of the population directly affects both the amount of carbon emissions and the demand for land and resources. As the population grows, so does the pressure on natural systems, resulting in increased emissions, land use changes, and greater consumption of energy. The rising population is a primary driver of industrial activities, agricultural expansion, and urbanization, all of which contribute to higher stages of carbon emissions [11], [12]. Additionally, population growth often leads to a reduction in forest area due to land being cleared for agriculture, settlements, and infrastructure development, further contributing to the depletion of carbon sinks [13], [14].

Forest area is another critical variable in this model. Forests help as vital carbon sinks, absorbing significant amounts of carbon dioxide from the atmosphere and storing it in biomass and soil. However, deforestation for agricultural increase, logging, and urban development decreases the volume of forests to act as carbon sinks [15].

Figure 1 describe the global warming and CO$_2$ the climate of the globe is greatly impacted by global warming, that is mainly due to the rising amount of greenhouse gases such as CO$_2$. In recent decades, there has been a substantial rise in world average temperatures, which has led to in numerous ecological and socioeconomic repercussions. The Figure 1 contains a concise summary of significant statistics regarding to the effects due to global warming between the years 1960 and 2020.

Global warming this model aims to provide a comprehensive understanding of how the carbon concentration in the atmosphere is affected by the growth of the human population, changes in forest area, and efforts to increase tree plantations [16]. Human activities, including industrial production, agriculture, deforestation, and transportation, release significant amounts of carbon dioxide into the atmosphere. These productions, driven by the growth of the human population and the increasing demand for energy and resources, lead to an accumulation of carbon in the air, raising concerns about the long-term sustainability of Earth’s climate systems [17].

By capturing the dynamic interactions between these variables, the model highlights the feedback loops that exist between population growth, carbon emissions, forest management, and tree plantation efforts [18], [19]. It underscores the importance of balancing human development with environmental conservation, particularly in the environment of combating climate change and sustaining ecological sustainability [20]. The key takeaway from this model is that achieving a sustainable future requires understanding and managing the complex relationship between carbon emissions, population growth, land use, and forest cover. The model provides insights into how human-driven factors such as population expansion and land use changes influence the environment, while also illustrating how natural systems like forests and tree plantations can act as mitigating forces. By examining the interactions between these variables, policymakers, environmentalists, and researchers can develop strategies to curb carbon emissions, restore forest ecosystems, and promote sustainable practices that will help mitigate the impacts of climate change [21], [22]. As forest area decreases, the ability of the natural environment to sequester CO$_2$ diminishes, leading to an increase in carbon concentrations. On the other hand, reforestation and afforestation efforts, particularly through tree plantations, can help alleviate this effect by increasing forest cover and enhancing carbon sequestration. Tree plantations, which involve the intentional planting of trees in areas where forests have been depleted or in land that was previously non-forested, provide a valuable strategy to offset carbon emissions [23].

In the following section, the mathematical formulation of the model is presented. The system is compartmentalized into four distinct classes, collectively denoted as CNFT, where C represents the atmospheric concentration of CO$_2$, N denotes the human population within the considered region, F corresponds to the forested land area, and T represents the area under active tree plantation. This classification forms the basis for the dynamic interactions captured within the model. The system of equations for the proposed model is follows:

with initial conditions:

$C(0)=C_0>0, N(0)=N_0>0, F(0)=F_0>0, T(0)=T_0>0$

By incorporating these properties, the proposed model is rendered analytically meaningful and suitable for capturing the essential dynamics of the coupled human–environment system.

Theorem 1. If the initial conditions of the model equations are non-negative, then the future solutions are also non-negative, i.e.: $C(t)>0, \quad N(t)>0, \quad F(t)>0, \quad T t)>0$.

Proof: Consider the differential equation (Eq. (1)), to demonstrate that $C(t)$ remains positive, there is:

Let $J=\lambda N-\alpha C-\left(\lambda_1+\frac{\gamma_1 T}{K_1+T}\right) F C$, then the differential Eq. (2) becomes:

in Eq. (4), where $K=e^{-J C_1}$ is a constant determined by the initial condition:

at $t=0$ (Eq. (5)) becomes $C(0)=\frac{Q-K}{J} \Rightarrow Q-K e^{-J t}>0$. This will hold true if $K \lt 0$, since $e^{J t}$ approaches 0 as $t \rightarrow \infty$. As $t \rightarrow \infty, C(t)$ approaches: $C(t) \rightarrow \frac{1}{J}$; this means that $C(t)$ stabilizes at a positive value as long as the parameters keep $J>0$. Using the method of separable variables, we have shown that $J(t)$ can be expressed in a way that ensures it remains positive for all $t \geq 0$, assuming $K\lt Q$ and $J>0$. Thus, $C(t)$ is positively invariant under these conditions.

Now to demonstrate that $N(t)$ remains positive, let's consider:

Let $\mathcal{D}_1=-\theta C N$. Thus, the differential Eq. (6) becomes:

$N(t) \geq e^{-D_1 t+C_1}$ at $t=0$, we have $N(0)=N_0$. Thus, in Eq. (7), $N(t) \geq N_0 e^{-D_1 t}$, we can conclude that $N(t)>0$ for all $t \geq 0$. Similarly, $F(t) \geq 0$ and $T(t) \geq 0$ for all $t \geq 0$. This means that if the initial conditions are non-negative, then the future solutions will also be non-negative. Thus, the proof is complete.

Theorem 2. All the solutions of the model equations are uniformly bounded and contained in a feasible region for all $t \geq 0$.

Proof: Let $\Gamma=\left\{(C(t), N(t), F(t), T(t)) \in R_+^4 \left\lvert\, 0 \leq N^*(t) \leq \frac{\Lambda}{k}\right.\right\}$ be the positive invariant set and let $N^*=C+N+F+T$ be the total population, $C, N, F$, and $T$ represent different compartments. Then:

To prove that $N^* (t)$ is bounded, we will solve the ordinary differential equation (ODE) (8):

$\frac{d N^*}{d t}=Q-k N^* \Rightarrow \frac{d N^*}{d t} \geq-k N^*$ with initial condition $N^*(0)=N_0^*$

First, we solve the corresponding ODE by using the basic method of integration:

where, $C_{!}$ is any constant of integration at $t=0$ then $N^*(t)$ becomes $N^*(t)=N_0^* e^{-k t}$ so that $N^*(t) \geq$ $N_0^* e^{-k t}$, which shows that $N^*(t)$ is a decreasing function, so that $N^*(t)$ bounded below. Now to show that $N^*(t)$ is also bounded above, consider the inequality $\frac{d N^*}{d t} \geq-k N^*$. By multiplying $e^{k t}$ on both sides and integrating from 0 to "$t$", we obtain $e^{k t} \frac{d N^*}{d t}+k e^{k t} N^* \geq 0$, which shows that $N^*(t) \leq N_0^* e^{k t}$ is bounded above. Combining both bounds, we conclude that $N_0^* e^{-k t} \leq N^*(t) \leq N_0^* e^{k t}$. Therefore, $N^*(t)$ is bounded for all $t \geq 0$.

Figure 2 presents a graphical validation of the theorem, which states that all solutions of the model equations are uniformly bounded and remain within a feasible region for all $t \geq 0$. The graph illustrates the contributions of each compartment $C(t), N(t), F(t)$ and $T(t)$ demonstrating the constraints $N_0 e^{-k t} \leq$ $N(t) \leq N_0 e^{k t}$ for $t \geq 0$.

In the framework of the model examining the impact of human population, forest management, and tree plantations on carbon emissions, the threshold parameter $R_0$ can be represented as a critical value that shows the balance between carbon emissions and carbon sequestration. When $R_0<1$, it suggests that carbon sequestration through natural forest processes and tree plantations is sufficient to offset emissions, leading to a sustainable system where carbon levels gradually decline. Conversely, if $R_0>1$, it indicates that carbon emissions exceed the sequestration capacity, causing a rise in atmospheric carbon and rendering the system unsustainable. This threshold also reflects forest sustainability: when human activities such as deforestation and poor land management are minimal ($R_0<1$), forests have the probable capacity to regenerate and continue absorbing carbon. However, if these harmful activities dominate $R_0>1$, forest degradation surpasses regrowth, resulting in increased emissions. Therefore, effective forest management and the implementation of tree plantations can play a vital role in reducing $R_0$ by enhancing carbon absorption and moving the system toward sustainability. To evaluate the threshold parameter $R_0$ for the proposed model using the next generation matrix method: $\frac{d X}{d t}=F(X)-V(X),$ where $ X=(C, N, F, T)$.

Now, the Jacobean matrices of $\mathcal{F}$ and $\mathcal{V}$, evaluated at the disease-free equilibrium $E_0$, are:

The threshold parameter $R_0$ is given by $R_0=\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)$. The equilibrium points $C_0, N_0, F_0$, and $T_0$ are essential for understanding the long-term behavior of the system. They represent the interactions between carbon emissions, human population growth, forest management, and tree plantation efforts. By analyzing and managing these equilibrium points, it is possible to guide the system toward sustainability, reduce carbon emissions, and promote effective carbon sequestration. These values also serve as a foundation for developing policies and strategies that address climate change, deforestation, and ecological resilience. To calculate the equilibrium points, we set the system derivatives to zero: $\frac{d C(t)}{d t}=\frac{d N(t)}{d t}=\frac{d F(t)}{d t}=\frac{d T(t)}{d t}=0$ Assuming that $C=N=F=T=0$, and after some simplification, we obtain the following equilibrium points for the proposed model:

\begin{align*} \left\{C_0, N_0 . F_0, T_0\right\}=\left\{\frac{Q \lambda N_0}{\alpha+\left(\lambda_1+\frac{\gamma_{1 T_0}}{k_{1+T_0}}\right) F_0}, K\left(1-\frac{\gamma \pi_{1 \theta} F_0}{\gamma}\right), \mu\left(1-\frac{F_0}{M}\right)-\phi N_0+\beta \frac{v\left(M-F_0\right)}{v_0 b}, \frac{v\left(M-F_0\right)}{v_0 b}\right\} \end{align*}

Figure 3 shows data points in three-dimensional space, allowing us to observe relationships between three variables. It uses a Cartesian coordinate system with x, y, and z axes to represent the data, making it possible to represent the relationships in a three-dimensional dataset.

In the section we check the stability of the model regardless of initial conditions, that all the solutions trajectories converges to equilibrium point as t tends to infinity. This analysis involve proof whether the system equilibrium point attracts all the solutions trajectories over time, it is used for long term behavior of the model. In mathematics point of view the global stability of an equilibrium point check by various mathematical approaches, such as geometric method, LaSalle’s invariance principle and Lyapunov functions, here we use Lyapunov functions to study the global stability of the corresponding equilibrium points of our model.

To study the stability, we linearize the system at the equilibrium point by computing the Jacobean matrix of the system. The Jacobean matrix $J$ describes how small changes in the state variables affect the time derivatives of the state variables. The Jacobean matrix is computed as follows:

$J=\left[\begin{array}{cccc} -\alpha-\left(\lambda_1+\frac{\gamma_1 T}{K_1+T}\right) F & \lambda & -\left(\lambda_1+\frac{\gamma_1 T}{K_1+T}\right) C & -C F \cdot \frac{\gamma_1 K_1}{\left(K_1+T\right)^2} \\ -\theta N & r\left(1-\frac{2 N}{K}\right)+\pi_1 \phi F-\theta C & \pi_1 \phi N & 0 \\ 0 & -\phi F & u\left(1-\frac{2 F}{M}\right)-\phi N+\beta T & \beta F \\ 0 & 0 & -v & -v_0 \end{array}\right]$

The resulting simplified Jacobean matrix at the equilibrium $C=N=F=T=0$ is:

We solve the characteristic equation:

$\operatorname{det}\left(J-\lambda^* I\right)=\left[\begin{array}{cccc}-\alpha-\lambda^* & \lambda & 0 & 0 \\ 0 & r-\lambda^* & 0 & 0 \\ 0 & 0 & u-\lambda^* & 0 \\ 0 & 0 & -v & -v_0-\lambda^*\end{array}\right], \lambda_1^*=-\alpha, \quad \lambda_2^*=r, \quad \lambda_3^*=u, \quad \lambda_4^*=-v_0$

To analyze the stability of the system, we examine the signs of the eigenvalues of the Jacobean matrix: $\lambda_1=-\alpha$ This is negative, which suggests that the system is stable with respect to changes in carbon concentration $C$. $\lambda_2=r$: This is positive (since $r>0$), which indicates that the population $N$ will grow exponentially and is unstable in that dimension. $\lambda_3=u$: This is also positive (since $u>0$), which means the forest $F$ will grow exponentially and is unstable in that dimension. $\lambda_4=-v_0$: This is negative, which suggests that the system is stable with respect to changes in tree plantation $T$.

Theorem 3. The proposed model is globally asymptotically stable at the disease-free equilibrium point if $R_0 \leq 1$; otherwise, it is unstable.

Proof: To show that the model is globally asymptotically stable when $R_0 \leq 1$, we use the following Lyapunov function:

$H(C, N, F, T)=\frac{1}{2}\left(\left(C-C_0\right)+\left(N-N_0\right)+\left(F-F_0\right)+\left(T-T_0\right)\right)^2$

To compute the time derivative of $H$, we need:

$\frac{d}{d t}\left(\left(C-C_0\right)+\left(N-N_0\right)+\left(F-F_0\right)+\left(T-T_0\right)\right)$

$ \frac{d H}{d t}=\left[Q+\lambda N-\alpha C-\left(\lambda_1+\frac{\gamma_1 T}{K_1+T}\right) F C\right]+\left[r N\left(1-\frac{N}{K}\right)+\pi_1 \phi N F-\theta C N\right]+\left[u F\left(1-\frac{F}{M}\right)-\phi N F+\beta F T\right]$

$\frac{d H}{d t}=(p-q)<0$ Thus, $\frac{d H}{d t}<0$ when $R_0 \leq 1$, which implies that $H$ is decreasing over time. Hence, the system is globally asymptotically stable at the disease-free equilibrium point if $R_0 \leq 1$. Otherwise, if $R_0>1$, the system is unstable.

Theorem 4. The proposed model is globally asymptotically stable at the endemic equilibrium point if $R_0>1$.

Proof: To show that the model is globally asymptotically stable, we will consider the following Lyapunov function:

$G\left(x_1, x_2, \ldots, x_n\right)=\frac{1}{2} \sum_{i=1}^n\left(x_i-x_i^*\right)^2$

where, $x_i=(C, N, F, T)$ and $x_i^*=\left(C^*, N^*, F^*, T^*\right)$.

First, we compute the time derivative of $G$: $\frac{d G}{d t}=\sum_{i=1}^n\left(x_i-x_i^*\right) \frac{d x_i}{d t}$. Let $N^*(t)=C+N+F+T$, and suppose:

$ \begin{aligned}& \frac{d G}{d t}=\left[N^*(t)-\left(C^*+N^*+F^*+T^*\right)\right] \frac{d G}{d t}(C+N+F+T) \cdot \frac{d G}{d t}(C+N+F+T)=Q-\theta N=0 \\ & \left(C^*+N^*+F^*+T^*\right)=\frac{\theta C^*-Q}{k_1} \cdot \frac{d G}{d t}=\left[N^*(t)-\frac{k_1 F^*-Q}{k_1}\right]\left[-\delta N(t)-\frac{k_1 C^*-Q}{k_1}\right] \\ & \frac{d G}{d t}=\left[N^*(t)-\frac{k_1 F^*-Q}{\theta}\right]\left[-\delta N^*(t)-\frac{\theta F^*-Q}{\theta}\right] \cdot \frac{d G}{d t}=-\theta\left[N^*(t)+\frac{Q}{\theta}\right]^2\end{aligned}$

Hence: $\frac{d G}{d t} \leq-\theta\left[N^*(t)+\frac{Q}{\theta}\right]^2<0$. Since $\frac{d G}{d t}<0$, all conditions of the Lyapunov function are satisfied, the model is globally asymptotically stable. The main idea of using a Lyapunov function in this model is to prove that the system naturally tends toward a stable equilibrium over time. By constructing a mathematical function that always decreases along system trajectories similar to how energy dissipates in a physical system, we can show that key variables like CO$_2$ concentration or forest area will not grow without bound but instead stabilize. This approach provides a rigorous way to confirm global stability without solving the system explicitly, makng it a powerful tool for understanding long-term environmental dynamics.

In this framework, sensitivity analysis plays a crucial role by examining how changes in these factors (e.g., population size, forest area, and tree plantation efforts) affect total carbon emissions. On occasion, an increase in population typically leads to higher emissions due to greater energy consumption and transportation, while enhanced forest management and development of tree plantations can reduce emissions by sequestering carbon. The sensitivity rule in this model raises the mathematical relationships that describe how sensitive carbon emissions are to changes in these variables. By calculating partial derivatives, the model can determine how small changes in population, forest management, and tree plantation areas influence emissions. This helps in identifying the most impactful factors for reducing emissions and informing policy decisions. The standard formula for the sensitivity index of $R_0$ is given by:

$\chi_{\xi}^{R_0}=\frac{\partial R_0}{\partial \xi} \times \frac{\xi}{R_0}$

where, $\xi$ can be any of the parameters in the set:

$\xi \in\left\{Q, \lambda, \alpha, \lambda_1, \gamma_1, K_1, r, K, \pi_1, \phi, \theta, u, M, \beta, v, v_0\right\}$

$\frac{\partial R_0}{\partial Q} \times \frac{Q}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial Q} \times \frac{Q}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}=1$

$\frac{\partial R_0}{\partial \lambda} \times \frac{\lambda}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial \lambda} \times \frac{\lambda}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx 0.1263$

$\frac{\partial R_0}{\partial \alpha} \times \frac{\alpha}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial \alpha} \times \frac{\alpha}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx 0.06731$

$\frac{\partial R_0}{\partial \gamma} \times \frac{\gamma}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial \gamma} \times \frac{\gamma}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx-0.4522$

$\frac{\partial R_0}{\partial K_1} \times \frac{K_1}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial K_1} \times \frac{K_1}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx-0.5234$

$\frac{\partial R_0}{\partial \mu} \times \frac{\mu}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial \mu} \times \frac{\mu}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx-0.3416$

$\frac{\partial R_0}{\partial \phi} \times \frac{\phi}{R_0}=\frac{\phi\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial \phi} \times \frac{\phi}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx 0.1291$

$\frac{\partial R_0}{\partial \theta} \times \frac{\theta}{R_0}=\frac{\partial\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)}{\partial \theta} \times \frac{\theta}{\left(\frac{Q+\lambda N_0}{\alpha+\lambda_1 F_0+\gamma_1 F_0 T_0}\right)\left(K_1+T_0\right)} \approx 0.3412$

From Figure 4, it is clear that the sensitivity index of $R_0$ with respect to $Q, \lambda, \alpha$ and $\phi$ indicates a direct relationship with $R_0$. On the other hand, the sensitivity index of $R_0$ with respect to $\gamma, K_1, \theta, \mu$ and $v$ shows an inverse relationship with $R_0$. The parameter $Q$ is more sensitive compared to the other parameters.

The aim of using optimal control in the model "Modeling the Impact of Human Population, Forest Management, and Tree Plantations on Carbon Emissions" is to progress a time-based strategy that successfully manages forest resources and human population pressures. By applying mathematical optimization techniques, the model purposes to identify the best combination of actions, such as reforestation, forest protection, and population-related efforts, that will minimize carbon emissions and support long-term ecological balance. This approach helps policymakers and environmental organizers make informed, efficient decisions to combat climate change and endorse sustainable development. In this section, we consider the rod accidents model to develop our control strategy. By implementing the two control variables, $u_1^*(t)$ is the control variable that represents the education/media campaign and $u_2^*(t)$ is the control variable that minimizes the Rod Accident.

With initial conditions:

$C(0)=C_0>0, N(0)=N_0>0, F(0)=F_0>0, T(0)=T_0>0$

In order to reduce the above-mentioned factors for this, the objective functional can be written as:

$J\left(u_1^*(t), u_2^*(t), v_3^*(t)\right)=\int_0^{T_f}\left(C_1^* C+C_2^* N+C_3^* F+C_4^* T+\frac{1}{2}\left(C_5^* u_1^2(t)+C_6 u_2^2(t)+C_7 u_3(t)\right)\right) d t$

The objective function is used for reducing and minimizing the number of accidents. The constants $C_m$ for $m$=1,2,3,…,4 are used as balancing cost factors, and $T_f$ represents the final time. Due to the nonlinear intervention among the population, we find an optimal control triplet as follows:

$J\left(u_1^*(t), u_2^*(t), u_3^*\right)=\min \left\{J\left(u_1^*(t), u_2^*(t)\right) \mid\left(u_1^*(t), u_2^*(t), u_3^*(t)\right) \in U\right\}$

In this section, we give proof for the existence of the control problem. We define the Hamiltonian $\widehat{H}$ for the optimal control problem as follows:

$\widehat{H}=L\left(C(t), N(t), F(t), T(t), u_1^*(t), u_2^*(t)\right)+\lambda_1^* \frac{d C(t)}{d t}+\lambda_2^* \frac{d N(t)}{d t}+\lambda_3^* \frac{d F(t)}{d t}+\lambda_4^* \frac{d T(t)}{d t}$

Theorem 5. For the control problem there exists $\mathbf{u}^*(t)=\left(u_1^*(t), u_2^*(t)\right) \in U$ such that:

$\min_{(u_1^*(t), u_2^*(t)) \in U} J\left(u_1^*(t), u_2^*(t)\right)=J\left(u_1^*(t), u_2^*(t),\right)$.

We brought several techniques into use for validation of the optimal control prevalence. Thus, all the control and state variables are nonnegative. That is why, in the process of reducing the problem, the required convexity of the objective functional is elaborated in equation, in $u_1^*(t)$ and $u_2^*(t)$ is gratified. The set of control variables $\left.u_1^*(t), u_2^*(t)\right) \in U$ is also convex and closed by definition. This optimal system is delineated, and it provides surety about the solidity that is needed for the validation of the optimal control system. Moreover, the integrand in the objective functional $C_1^* A+C_2^* B+C_3^* D+C_4^* C+\frac{1}{2}\left(C_5^* u_1^2(t)+C_6 u_2^2(t)\right)$ is convex on the control set $U$, which certifies the proof.

Theorem 6. Given optimal controls $u_1^*(t), u_2^*(t)$ and solutions $C^*, N^*, F^*, T^*$ of the corresponding state system, there exist adjoint variables $\lambda_m^*(t)$, for $m=1, \ldots, 4$, satisfying the following equations:

$\frac{d \lambda_1^*}{d t}=b A\left(\lambda_1^*-\lambda_2^*\right)+\left(\lambda_1^*-\lambda_4^*\right)\left(1-u_1^*(t)\right)+m \lambda_1^*$

$\frac{d \lambda_2^*}{d t}=-C_1^*+\left(\lambda_2^*-\lambda_3^*\right) a_1+\left(\lambda_2^*-\lambda_4^*\right)\left(1-u_1^*(t)\right)+m \lambda_2^*$

$\frac{d \lambda_3^*}{d t}=-C_2^*+\left(\lambda_1^*-\lambda_2^*\right) b S+\lambda_3^*(m+g+r+h)+\left(\lambda_3^*-\lambda_4^*\right)\left(1-u_2^*(t)\right)$

$\frac{d \lambda_4^*}{d t}=-C_3^*+\lambda_4^*(m+t+f+h)+\left(\lambda_4^*-\lambda_4^*\right)\left(1-u_2^*(t)\right)+\lambda_3^* t+\lambda_4^* f$

Proof: If we take the values as $A(t)=A^*, B(t)=B^*, D(t)=D^*, C(t)=C^*$ and differentiate the Hamiltonian with respect to state variables $A(t), B(t), D(t), C(t)$, respectively, we get the adjoint system:

$\frac{d \lambda_m^*(t)}{d t}=-\frac{\partial H}{\partial x_m}$, for $m=1,2, \ldots, 4$

where, $\lambda_m^*(t)$ represents the adjoint variables associated with the state variables $x_m(t)=$ $A(t), B(t), D(t), C(t)$. The transversal conditions are given by: $\lambda_m^*\left(t_f\right)=0, \quad m=1,2, \ldots, 4$, where $t_f$ is the final time.

Theorem 7. The control pair $\left(u_1^*(t), u_2^*(t), u_3^*(t)\right)$, which minimize the objective functional $J$ over the region $U$, is given by:

$\mathrm{u}_1^*(\mathrm{t})=\min \left\{1, \max \left(\frac{\left(\lambda_6^*-\lambda_1^*\right) \mathrm{C}+\left(\lambda_6^*-\lambda_2^*\right) \mathrm{N}}{\mathrm{C}_5^*}, 0\right)\right\}$

$\mathrm{u}_2^*(\mathrm{t})=\min \left\{1, \max \left(\frac{\left(\lambda_6^*-\lambda_3^*\right) \mathrm{F}+\left(\lambda_6^*-\lambda_5^*\right) \mathrm{T}}{\mathrm{C}_6}, 0\right)\right\}$

Proof: By using the optimality condition, we get:

$\frac{\partial \mathrm{C}}{\partial \mathrm{u}_1}=\mathrm{C}_5^* \mathrm{u}_1^*(\mathrm{t})+\left(\lambda_6^*-\lambda_1^*\right) \mathrm{A}+\left(\lambda_6^*-\lambda_2^*\right) \mathrm{N}$

$\frac{\partial \mathrm{C}}{\partial \mathrm{u}_2}=\mathrm{C}_6 \mathrm{u}_2^*(\mathrm{t})+\left(\lambda_3^*-\lambda_6^*\right)+\left(\lambda_5^*-\lambda_6^*\right) \mathrm{F}, \frac{\partial \mathrm{C}}{\partial \mathrm{u}_3}=\mathrm{C}_7 \mathrm{u}_3(\mathrm{t})+\left(\lambda_4^*-\lambda_6^*\right) \mathrm{T}$

The optimal control variables $\mathrm{u}_1^*(\mathrm{t}), \mathrm{u}_2^*(\mathrm{t}), \mathrm{u}_3^*(\mathrm{t})$ used to solve Eq. (12) are:

$\mathrm{u}_1^*(\mathrm{t})=\frac{\left(\lambda_6^*-\lambda_1^*\right) \mathrm{C}+\left(\lambda_6^*-\lambda_2^*\right) \mathrm{N}}{\mathrm{C}_5^*}, \mathrm{u}_2^*(\mathrm{t})=\frac{\left(\lambda_6^*-\lambda_3^*\right) \mathrm{F}+\left(\lambda_6^*-\lambda_5^*\right) \mathrm{N}}{\mathrm{C}_6}$

The property of the control space equations can be written as:

$\mathrm{u}_1^*(\mathrm{t})=\left\{\begin{array}{lll}0, & \text { if } & \frac{\left(\lambda_6^*-\lambda_1^*\right) \mathrm{C}+\left(\lambda_6^*-\lambda_2^*\right) \mathrm{N}}{\mathrm{C}_3^*} \leq 0 \\ \frac{\left(\lambda_6^*-\lambda_1^*\right) \mathrm{C}+\left(\lambda_4^*-\lambda_2^*\right) \mathrm{N}}{\mathrm{C}_3^*}, & \text { if } & 0<\frac{\left(\lambda_4^*-\lambda_1^*\right) \mathrm{C}+\left(\lambda_4^*-\lambda_2^*\right) \mathrm{F}}{\mathrm{C}_5^*}<1 \\ 1, & \text { if } & \frac{\left(\lambda_4^*-\lambda_1^*\right) \mathrm{C}+\left(\lambda_6^*-\lambda_2^*\right) \mathrm{N}}{\mathrm{C}_5^*} \geq 1\end{array}\right.$

$\mathrm{u}_2^*(\mathrm{t})= \begin{cases}0, & \text { if } \quad \frac{\left(\lambda_4^*-\lambda_3^*\right) \mathrm{D}+\left(\lambda_4^*-\lambda_3^*\right) \mathrm{F}}{\mathrm{C}_4^*} \leq 0 \\ \frac{\left(\lambda_4^*-\lambda_3^*\right) \mathrm{T}+\left(\lambda_4^*-\lambda_3^*\right) \mathrm{N}}{\mathrm{C}_4^*}, & \text { if } \quad 0<\frac{\left(\lambda_4^*-\lambda_3^*\right) \mathrm{F}+\left(\lambda_4^*-\lambda_5^*\right) \mathrm{N}}{\mathrm{C}_4^*}<1 \\ 1, & \text { if } \quad \frac{\left(\lambda_6^*-\lambda_3^*\right) \mathrm{F}+\left(\lambda_4^*-\lambda_5^*\right) \mathrm{N}}{\mathrm{C}_4^*} \geq 1 \end{cases}$

According to compact notation, $u_1^*(t)$ and $u_2^*(t)$ can be written as:

$u_1^*(t)=\min \left\{1, \max \left(\frac{\left(\lambda_4^*-\lambda_1^*\right) A+\left(\lambda_3^*-\lambda_2^*\right) N}{C_3^*}, 0\right)\right\}, u_2^*(\mathrm{t})=\min \left\{1, \max \left(\frac{\left(\lambda_4^*-\lambda_3^*\right) F+\left(\lambda_4^*-\lambda_3^*\right) N}{C_6}, 0\right)\right\}$

Figure 5 illustrates the dynamics of the proposed model under different control strategies. The red line represents the scenario with control factors applied, while the blue line shows the system’s behavior without control. Also show regulates a behavior of dynamic processes to achieve desired outcomes. This is where control strategies come into play. These strategies are designed to influence the system’s inputs in order to maintain stability, optimize performance, or achieve specific goals despite disturbances or uncertainties.

In this study, we developed and applied a numerical scheme to model the dynamics of a coupled environmental system using a system of ordinary differential equations (ODEs). By implementing the classical 4th-order Runge-Kutta (RK4) method, we simulated and analyzed the evolution of various state variables over time, including the concentration of $\mathrm{CO}_2$, the human population $N$, forest area $F$, and the number of planted trees $T$. For the numerical simulations, we used the following parameter values: $Q=1.78, \lambda=0.447, \alpha=0.003$, $\lambda_1=5.581 \times 10^{-7}, r=0.126, K=11, \theta=5.3765 \times 10^{-8}, u=0.005, M 5900, \phi=0.300371$, $\pi_1=0.1005, \gamma_1=3 \times 10^{-7}, K_1=200, \beta=3 \times 10^{-6}, v=0.2004$, and $v_0=0.1$.

This simulation defines the dynamic interaction among human population, forest area, tree plantations, and carbon dioxide (CO$_2$) emissions over time, based on the proposed CNFT model. The results are demonstrated in Figure 6 which contains four sub figures representing the evolution of key variables. In subgraph (a) of Figure 6 we observe the concentration of CO$_2$ increasing steadily with time $t$. This trend reflects the impact of rising human population and decreasing forest cover. As the human population $N(t)$ grows, it leads to an increased carbon emission rate $\lambda N$, while simultaneous deforestation $\varphi N F$ reduces the ability of forests to absorb CO$_2$, thus elevating atmospheric carbon levels. Subgraph (b) of Figure 6 displays the growth of the human population over time. It shows a logistic-type increase, influenced by the carrying capacity $K$ and environmental response. However, as CO$_2$ levels increase, the population growth rate may be affected destructively due to the term -$\theta C N$, which models carbon-related stress or mortality. In subgraph (c) of Figure 6 the forest area $F(t)$ is shown to decline initially, primarily due to human-induced deforestation. Over time, there is a slight recovery, influenced by natural forest growth and contact with tree plantations. The competition between deforestation and reforestation dynamics is clearly visible here. Subgraph (d) of Figure 6 illustrates the evolution of tree plantations $T(t)$. Initially, tree plantations grow due to efforts to mitigate forest loss and sequester carbon, governed by the term $v(M-F)$. However, carbon emissions from plantation activities $(v_0 T)$ and limited forest area slow down this growth over time. Overall, the simulation highlights that without effective forest management and active plantation efforts, CO$_2$ concentration will continue to rise due to increasing human activity and shrinking natural sinks. This underscores the need for strategic intervention in land use and environmental policies to ensure long-term sustainability and climate resilience.

This article provides a complete analysis of CO$_2$ concentration dynamics using both theoretical and computational techniques. The study presents a numerical approach to model the dynamics of a coupled environmental system through a system of ordinary differential equations. The simulation results demonstrate how factors such as deforestation, population growth, and reforestation efforts influence CO$_2$ concentration and overall environmental balance. The conclusions highlight the critical role of sustainable practices in maintaining environmental stability. This model aims to enhance understanding of how atmospheric carbon levels are affected by human populace growth, forest area changes, and tree plantation efforts. By capturing the dynamic interactions among these variables, the model highlights feedback loops between population growth, carbon emissions, forest management, and reforestation. It underlines the need to balance human development with environmental conservation, mostly in the context of combating climate change and promoting ecological sustainability.