2. Materials and Methods
2.1. The Mathematical Model
Referring to the models of [33,34], we created a mathematical model that studies optimal rotation strategies that maximize carbon sequestration in forest plantations. This model
is based on a system of three ordinary differential equations (ODEs) that are governed
by three state variables, B(t), r(t), and I(t), that denote the amount of living biomass
(considers aboveground and belowground biomass), intrinsic biomass growth, and burned
area, respectively. From the state variables, four control strategies are associated: reforestation R(t), felling F(t), fire prevention S(t), and thinning T(t).
In the following, the
notation “dot” represents the derivative of a variable with respect to time t. The dynamics
of the live biomass has a logistic growth with a carrying capacity K. The biomass also
increases proportionally and indirectly with respect to reforestation activities and decreases
proportionally due to the immediate effects of fire, felling, and thinning. The contribution
of relative humidity is not considered in the biomass dynamics, since it exceeds the carrying
capacity [33]. The differential equation governing biomass is as follows.
.
B(t) = r(t)B(t)
�
1 −
B(t)
K
�
+ [βR(t)]B(t) − [µ1 I(t) + σF(t) + τ T(t)]B(t). (1)
From Equation (1), β is the rate at which biomass increases with respect to reforestation,
µ1 is the rate at which biomass decreases due to fire effects, σ is the rate at which biomass
decreases due to felling effects, and τ is the rate at which biomass decreases due to thinning
(there is an instantaneous decrease in biomass).
Intrinsic growth was originally studied in [34]. Here, we consider that thinning has
an indirect contribution with respect to individual growth in a linear manner, since, in the
long term, it is related to biomass to control density and timber quality [35], such that
.
r(t) = r0 − ρr(t) + νT(t). (2)
In Equation (2),
r0 represents the maximum individual growth rate under ideal
conditions [36] and ρ is the effect of the natural mortality rate on individual growth. In the
case of thinning, in the long term, this silvicultural operation enhances the growth of the
remaining trees, and ν is the parameter by which the intrinsic growth of the living biomass
is increased by the effects of thinning.
Finally, the burned area increases as biomass increases, but this increasing behavior
cannot be unlimited. Therefore, a relationship between burned area and biomass is considered that limits the growth of burned area without inhibiting fire. On the other hand,
Forests 2023, 14, 82 4 of 17
burned area decreases due to fire prevention, thinning, and relative humidity. Thus,
the
dynamics of the burned area are
.
I(t) = µ2 I(t)
�
B(t)
1 + B(t)
�
− θS(t) − ηT(t) − hI(t). (3)
From Equation (3), µ2 is the fire rate, θ is the fire prevention rate, η is the thinning
rate, and h is the relative humidity threshold at which fire occurs.
Considering Equations (1)–(3), the mathematical model is presented by means of a
system of nonlinear ODEs, as follows
.
B(t) = r(t)B(t)
�
1 −
B(t)
K
�
+ [βR(t)]B(t) − [µ1 I(t) + σF(t) + τ T(t)]B(t)
.
r(t) = r0 − ρr(t) + νT(t)
.
I(t) = µ2 I(t)
�
B(t)
1+B(t)
�
− θS(t) − ηT(t) − hI(t),
(4)
satisfying the initial conditions B(t0) = B0, r(t0) = r0, and I(t0) = I0. It is assumed that
the amount of carbon stored in the biomass is proportional to the amount of biomass [37],
therefore C(t) = αB(t) where α is the rate of carbon capture. The following set of assumptions allows building a model that simulates control strategies that maximize carbon
capture in forest plantations:
• The model is formulated for fast-growing managed forest plantations;
• There are many plantations with different ages, thus biomass never goes to zero;
• The ambient humidity is considered constant for simplicity;
• No soil fertilization in each cycle of the forest regeneration is considered;
• The area burned per year is considered, but human intentionality is not taken into account;
• There are no incentives for reforestation or carbon capture;
• The harvesting method corresponds to clear-cutting;
• In the thinning, the thinner, lower quality and less commercially valuable trees will be
removed. Two types of thinning effects are considered, which are explained below;
• The presence of artificial irrigation is neglected in the model;
• The budget for fire prevention is limited;
• Intensive management of forestry is not included in our model;
• Trees burned by fire are replaced by new plants and natural regeneration is not used;
• The mortality rate of extreme events is neglected in our model.
Where the variables B(t) ≥ 0,r(t) ≥ 0, I(t) ≥ 0, R(t) ≥ 0, F(t) ≥ 0, S(t) ≥ 0 and
T(t) ≥ 0, the parameters (K, σ, h, µ1, µ2, τ, β, r0, ρ, η, θ, α, ν) ∈ R
13
+ .
The relationships between the variables and parameters used in the controlled mathematical model (4) are schematically represented in Figure 1.
Their notations, definitions, and units for variables and parameters are described in
Tables 1 and 2, respectively.