A mathematical model to study the dynamics of carbon capture in forest plantations

Fast-growing forest plantations play an important role in reducing global warming and have great potential for carbon capture. In this study, we aimed to model the dynamics of carbon capture in fast-growing plantations. A mathematical model is proposed consisting of a tridimensional nonlinear system. The variables involved are the amount of living biomass, the intrinsic growth of biomass, and the burned area by forestry fire. The environmental humidity is also considered, assumed as a parameter by simplicity. The solutions of the model are approximated numerically by the Runge-Kutta fourth-order method. Once the equilibria of the model have been obtained and its local stability determined, the analysis of the model reveals that the living biomass, as well as the stored carbon, decreases in each harvest cycle as a consequence of the negative effects of fire on soil properties. Furthermore, the model shows that the maximum area burned is attained always after the maximum volume of biomass is obtained. Numerical simulations show that the model solutions are reasonable for the growth dynamics of a plantation, from a theoretical perspective. The mathematical results suggest that a suitable optimal management strategy to avoid biomass losses in the successive regeneration cycles of the plantation is the prevention of fires together with soil fertilization, applied to fast-growing plantations.


Introduction
Carbon dioxide (CO 2 ) is one of the main greenhouse gases (GHG) in the atmosphere; in the period of 1750-2019, the CO 2 concentration has increased from 280 to 380 parts per million (ppm) and it is estimated that this increase could raise the global mean temperature [1,2], generating adverse effects such as an increase in the sea level, increased frequency of fires, decreased productivity in global agriculture, and will even cause severe negative impacts on the economy and human survival [3]. Due to the above, governments, through cooperation agreements (e.g., the Paris agreement and the Kyoto Protocol), are creating new forest areas to prevent the global average temperature from increasing over the next century above 2 • C [4].
Worldwide, the total area of plantations currently reaches 294 million hectares; of this surface, approximately 45% are forest plantations, that is, intensively managed forests, and the other 55% are native species [5]. Because forested areas store the highest amount of atmospheric carbon [6], it has been suggested that expansion of these areas and the increasing of the harvest age; especially in fast-growing forest plantations [7], are important strategies to maximize the timber yield and also to increase the carbon capture and the mitigation of global climate warming [8].
In some climatic and site conditions it is risky to maintain the forest for longer periods as a strategy to increase carbon capture. This is because the probability of forest fires occurrence increases as there is more fuel present in the field.
The problem of forest fires has been incorporated into the analyzes of various forest management problems using different models [9]. To model the probability of occurrence of forest fires, some authors used the Faustmann model generalized to Poisson's stochastic process [10,11]. By analyzing a mathematical model of non-linear differential equations utilizing the Hopf bifurcation [12] showed that the CO 2 concentration decreases as the reforestation rate increases. In an oil palm study [13], optimal control theory is applied to maximize oil production and carbon absorption. The authors found that longer harvest ages are optimal to maintain biomass and maximize carbon uptake.
However, despite that there are mathematical models originally presented in [13,14], which incorporates the dynamics of living biomass and the intrinsic growth of biomass and consider the felling rate as a control variable, they do not analyze variables such as the burned area and the ambient humidity in joint form. In this context, the novelty of this research will be to jointly analyze the dynamics of living biomass, the intrinsic growth of biomass, and the burned area, under different environmental humidity conditions. The main goal of this work is to model the dynamics of carbon capture in fast-growing forest plantations by using a mathematical model of non-linear ordinary differential equations.

Mathematical model
A non-linear mathematical model is presented to model the dynamics of carbon capture in fastgrowing forest plantations. Given at any time t and in any region, it is considered that B(t), r(t) and I(t) denote the amount of living biomass (it considers the above and below-ground biomass of the tree), the intrinsic growth of biomass and the burned area per year, respectively. In addition, we make the following assumptions for our model.
(i) The biomass follows a logistic growth model with a carrying capacity K.
(ii) There is a variable intrinsic growth rate. (iii) The biomass increases proportionally by the environmental humidity (assumed constant for simplicity) and decreases proportionally by the presence of fire. (iv) Soil fertilization and reforestation are not considered.
Under the above assumptions, the dynamics of the mathematical model is governed by the following system of non-linear ordinary differential equations (ODE), as expressed in the system of Equations (1).
Satisfying initial conditions B(t 0 ) = B 0 , r(t 0 ) = r 0 y I(t 0 ) = I 0 . The amount of carbon stored in biomass is assumed to be proportional to the amount of biomass, therefore C(t) = αB(t) [13]. The model applies to a forest with few silvicultural operations. In this model h 1 is the environmental humidity parameter and µ 1 is the rate at which the biomass decreases due to the effects of fire. In the case of intrinsic growth r(t), it was originally proposed in [14], where r 0 represents the maximum rate of individual growth under ideal conditions and ρ is the effect of the natural mortality rate for individual growth. Finally, the burned area increases as biomass increases, but this increasing behavior cannot be unbounded. For this reason, a relationship between the burned area and the biomass is considered that limits the growth of the burned area without spread the fire, µ 2 is the fire rate, and h 2 is the threshold of ambient humidity that reduces the incidence of fires, all parameters are positive. The variables are all positive and their notation, definition, and units are shown in Table 1. Table 1. Notation, definition and units of each variable.

Notation
Definition Unit

B(t)
Volume of living biomass Intrinsic growth of biomass m 3 ha −1 yr −1

I(t)
Burned area per year m 2 yr −1 The region of attraction [15] for Equation (1), which shows all solutions of the Equation (1) are bounded, is given by the following result.
The set Ω is the region of attraction for the dynamic system of the Equation (1), is positively invariant and attracts all the solutions of the Equation (1) with initial conditions in R 3 + .
Proof. Applying the upper limit to the first and second equations of the system of the Equation (1), it has to be bounded [15], we have B m = K + h 1 r 0 , r m = r 0 ρ . Now, it is determined that the function Z is bounded. Replacing the derivatives of the Equations (1), we obtain Z ≤ Bm(rm+h 1 ) µ 1 (1+Bm) + r m ρ; therefore, the attraction set Ω is bounded.

Analysis of the model
This section deals with the dynamical behaviour of the Equation (1). The existence of the equilibrium points of the system and the stability properties of the equilibria are determined.

Equilibrium points at the boundary
Equation (1) has the following two points of feasible frontier equilibrium, these are:

Equilibrium points at the interior
The interior equilibrium P * (B * , r * , I * ), is obtained by setting of the Equation (1) equal to zero, such that is a sufficient condition for the existence of a positive equilibrium.

Stability of equilibrium points
In this section, the local asymptotic stability of the equilibrium points mentioned above is determined. The stability of an equilibrium point of the Equation (1) is related to the signs of eigenvalues of its corresponding Jacobian matrix; therefore, to analyze the stability at the points, we will linearize the system around each equilibrium point. Also, we will apply the Routh-Hurwitz criterion [16,17] for obtain the sufficient conditions for local stability of these points. The Jacobian matrix of the Equation (1) is given in Equation (2).
The linearization technique is used to study the behavior of the system in a neighborhood of their equilibrium points to obtain a qualitative understanding of the solutions, using the previous Jacobian matrix in the Equation (2). Proof. First, the stability of point P 1 is determined; with the Jacobian matrix, J evaluated in P 1 the eigenvalues are calculated. We have the Equation (3).
Therefore of the Equation (3), P 1 is a saddle point, locally unstable in the I direction and locally stable in the B − r plane. Now, the stability of point P 2 will be determined. The eigenvalues of matrix J evaluated in P 2 are given in Equation (4).
Therefore of the Equation (4), P 2 is locally stable in the B − r plane and locally stable in the I direction as long as µ 2 < h 2 + r 0 h 2 K(h 1 ρ+r 0 ) . It is locally unstable in the I direction if the opposite is true µ 2 > h 2 + r 0 h 2 K(h 1 ρ+r 0 ) . Now we study the stability at point P * , Jacobian matrix evaluated at point J(P * ) the matrix entries are: j 21 = j 32 = j 33 = j 23 = 0, j 22 = −ρ, the characteristic equation is determined by det(J(P * ) − ψI) = 0, then Equation (5) is obtained.
the coefficients of the characteristic equation are in Equation (6) and Equation (7).
In the following theorem, we will demonstrate the stability of the interior point for Equation (1 Proof. If µ 2 > h 2 , then it satisfy that D 1 > 0. On the other hand, if µ 2 > h 2 and K > h 2 µ 2 −h 2 , then it satisfy that D 3 > 0. Furthermore, it is known that in the Equation (8).

Numerical simulation
This section presents the numerical solution of the model system, Equation (1). One of the most used methods to numerically solve ODE problems with initial conditions is the fourthorder Runge-Kutta method, which behaves well in its approximation [18].
In the Figure ( Figure 1(a), shows that the biomass reaches a maximum in year 15 approximately, then decreases due to fires; while in Figure 1(b), it shows that the burned area occurs years after the maximum biomass volume, then begins to decrease in an oscillatory manner as the biomass decreases, but reaches a point where it remains constant. In Figure 1(c), the intrinsic growth has an asymptotic growth. Finally, in Figure 1(d), the behavior of carbon sequestration is observed, which is half of the total biomass following the same behavior of the biomass. Subject to the previous conditions, considering the fixed parameters K = 400, µ 1 = 1.   The model indicates that the volume of biomass has a maximum of 185 m 3 ha −1 by the age of 10 years, and this biomass decreases over time as a consequence of fires Figure 2(a) and Figure 2(b). Although the model indicates that the plantation regenerates in each growth cycle, it only reaches a maximum volume of 52 m 3 by the age of 45 years, observing a decrease in biomass in each subsequent cycle Figure 2(a). This might be because the soil fertility decreases in each cycle; with the consequent negative effect on the growth and accumulation of biomass in the plantation, because the constant presence of fires Figure 2(b) causes a decrease in the main soil minerals [19]. Likewise, since the model does not include the practice of reforestation, it is assumed that in each regeneration cycle there will be a high competition between the new regenerated plants, which will decrease their growth potential.
A direct relationship is observed between the accumulated burned area and biomass per year, that is, the smaller the area burned, the higher the biomass; the behavior of the burned area is quite regular Figure 2(b), and follows a similar trend to that of biomass. In the following cycles of the plantation, the maximum amount of biomass is reached at age 45 to 50 years, while the maximum of burned area is reached at age 40 to 60 years. In other words, for a burned area to exist, the biomass must first exist. The intrinsic growth rate of biomass obtained shows an asymptotic behavior, this means that when the biomass reaches its maximum, the intrinsic growth is slow Figure 2(c). Finally, it is observed that the dynamics of carbon capture follows the same decreasing trend as that shown by biomass Figure 2(d). This is because the model assumes that the carbon stored by the tree is a fraction of all living biomass. On average, it is accepted that the stored carbon content of a tree corresponds to c.a. 50% of its biomass [13].

Discussion
In this work, we studied the local stability of a mathematical model that analyzes the dynamics of carbon stored in fast-growing exotic species and provide new insights into the modeling of carbon capture dynamics in forest plantations. We start from a mathematical model proposed in [14], in which the authors demonstrate how the sustainability of the lethal and non-lethal harvest depends on the demographic effect of each type of harvest on growth. Furthermore, in [20] by applying the optimal control theory and assuming the mathematical model presented in [14], the authors showed that slow-growing species have optimal collection rates and lower profits than fast-growing species. In [13] they studied the dynamics of maintaining biomass at a certain level and maximizing oil production and carbon absorption at longer harvest ages. However, they do not include the probability of natural or anthropogenic disasters, such a foresty fires. We extended the model presented in [14] and included the burned area as a new variable.
The model confirms that the influence of fires causes a decrease in the volume of biomass, which in turn causes a reduction in carbon captured by trees and might imply an increase in the atmospheric CO 2 concentration. The inclusion of fire in the model could help forest owners' decision making as they can decide when to harvest their forests in order to minimize financial losses from fires. Our model corroborates that the presence of fires and the lack of artificial fertilization to the soil negatively affect the growth of biomass during each regeneration cycle.
Optimal strategies to achieve maximum biomass yield and carbon capture consistently over time should be considered for this purpose, could be the application of silvicultural management operations such as fertilization and fire prevention. However, this assertion needs further research; due to the contribution of physics to our work (CO 2 flux) it was possible to understand the importance of trees as carbon sinks, which can be useful for civil society to value the contribution of fast-growing forest plantations on preventing global warming.

Conclusion
The model studied confirms that the dynamics of stored carbon is similar to the growth of biomass and that both reach a maximum in the first 10 years of the plantation life cycle. Therefore, the consideration of the fixed parameters can help to predict how the influence of fires in forest plantations affects the dynamics of carbon capture.