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.