Optimal Management Strategies to Maximize Carbon Capture in Forest Plantations: A Case Study with Pinus radiata D. Don

Abstract

Plantations with fast-growing species play a crucial role in reducing global warming and have great carbon capture potential. Therefore, determining optimal management strategies is a challenge in the management of forest plantations to achieve the maximum carbon capture rate. The objective of this work is to determine optimal rotation strategies that maximize carbon capture in forest plantations. By evaluating an ecological optimal control problem, this work presents a method that manages forest plantations by planning activities such as reforestation, felling, thinning, and fire prevention. The mathematical model is governed by three ordinary differential equations: live biomass, intrinsic growth, and burned area. The characterization of the optimal control problem using Pontryagin’s maximum principle is analyzed. The model solutions are approximated numerically by the fourth-order Runge–Kutta method. To verify the efficiency of the model, parameters for three scenarios were considered: a realistic one that represents current forestry activities based on previous studies for the exotic species Pinus radiata D. Don, another pessimistic, which considers significant losses in forest productivity; and a more optimistic scenario which assumes the creation of new forest areas that contribute with carbon capture to prevent the increase in global temperature. The model predicts a higher volume of biomass for the optimistic scenario, with the consequent higher carbon capture than in the other two scenarios. The optimal solution for the felling strategy suggests that, to increase carbon capture, the rotation age should be prolonged and the felling rate decreased. The model also confirms that reforestation should be carried out immediately after felling, applying maximum reforestation effort in the optimistic and pessimistic scenarios. On the other hand, the model indicates that the maximum prevention effort should be applied during the life cycle of the plantation, which should be proportional to the biomass volume. Finally, the optimal solution for the thinning strategy indicates that in all three scenarios, the maximum thinning effort should be applied until the time when the fire prevention strategy begins.

1. Introduction

Carbon dioxide (CO₂) is one of the main greenhouse gases (GHG) in the atmosphere. Multiple human activities in most industrialized countries have contributed to the increase in this gas and have exacerbated the negative effects of climate change. According to the latest report of the Intergovernmental Panel on Climate Change (IPCC), climate change is devastating today, in particular, because of the changes in the patterns of humidity, temperature, winds, snow, and ice, especially in coastal zones. These changes in climate conditions could have negative impacts on human health, agriculture, and the economy [1,2,3]. Under this worldwide situation, governments are making cooperative efforts agreements (e.g., the Paris Agreement and the Kyoto Protocol) to create new forest areas to help prevent the global average temperature rising more than 2 °C during the 21st century [4,5,6]. Forest ecosystems cover approximately 4100 billion hectares of the Earth’s surface and have a huge potential for carbon capture [7]. Of this total area, approximately 45% are exotic plantations whereas the other 55% corresponds to native forests [8]. Because forest ecosystems can store the largest amounts of carbon [9], it has been suggested that expanding forest areas and prolonging the rotation age (i.e., the growth period required to derive maximum value from a stand of timber), especially in exotic forest plantations [10], are key strategies to maximize carbon capture and mitigate the negative effects of global climate change [11]. There is a large body of literature where carbon capture is estimated [12,13]. In a temperate forest in Southern Europe, the aboveground carbon capture in the species Eucalyptus nitens (Deane and Maiden), Eucalyptus globulus Labill, and P. radiata, with rotation ages ranging from 10 to 35 years, was estimated to be from 443 to 634 Tn C ha⁻¹ [14]. The carbon sequestration with the same species established in Chile was 212 Tn C ha⁻¹ for P. radiata, 180 Tn C ha⁻¹ for E. nitens, and 117 Tn C ha⁻¹ for E. globulus (age of 20–24 years for P. radiata and 10–14 years for Eucalyptus) [15]. On the other hand, in Panama the carbon stored in Tectona grandis E.L (Teca) plantations during 1 and 10 years was estimated to be 2.9 Tn C ha⁻¹ and 40.7 Tn C ha⁻¹, respectively [16].

On the other hand, there are studies on the oil palm (Elaeis guineensis Jacq) which, due to its high biomass production and expansion dynamics, plays an important role in carbon capture [17]. By means of mathematical modelling the dynamics of both oil production and carbon capture have been studied [18]. In [19], they formulated an optimal control problem based on a system of ordinary differential equations that relate the dynamics of young and mature trees and considers felling as a control variable. The authors concluded that palm oil production and carbon capture increases with a controlled felling rate.

Notwithstanding, to increase CO₂ capture the trees must remain for longer periods in the field, which delays the rotation age [20,21]. However, in some situations, it is risky to prolong the rotation age in order to increase carbon capture, since it increases the probability of forest fires when there is more fuel in the field. More frequent forest fires will increase CO₂ levels in the atmosphere, causing extreme climate events and decreasing relative humidity in many regions of the world [22]. To model the probability of forest fire occurrence some authors have used the Faustmann model generalized to the stochastic Poisson process [23], whereas others have studied this phenomenon by using the Bellman equation to determine the optimal rotation age in a forest stand that produces timber and carbon benefits under fire risk [24]. The authors showed that higher fire risk will reduce the optimal rotation age due to a lack of fire prevention and low carbon prices, while a higher carbon price will increase the rotation age, thus obtaining a higher ecological benefit. It is known that fires contribute to the increase of CO₂ in the atmosphere. In [25] they developed a meteorological fire index to predict the risk of fire occurrence and help forest managers take appropriate preventive measures. The authors determined that relative humidity is a simple and feasible parameter to describe the occurrence of fires. Several mathematical models have been developed to describe the dynamics of CO₂ capture in reforestation projects [26,27,28]. The atmospheric CO₂ concentration decreases as the rate of reforestation increases. Also in [29], they presented a study to model the greenhouse effect caused by CO₂ emissions through the optimal control theory. In the model, the authors addressed the optimization of investments in reforestation and clean technologies associated with state variables such as CO₂ emissions, planted area, and Gross Domestic Product (GDP). They concluded that it is more efficient to invest in reforestation than in clean technologies.

Because forested areas can contribute to climate change mitigation, it is necessary to find optimal management strategies that maximize carbon capture. Strategies such as large-scale reforestations are efficient in capturing huge amounts of carbon [30], whereas the optimization of thinning, fire prevention, and harvesting strategies can also reduce CO₂ emissions in forest plantation management [31]. In [32] they applied a thinning strategy in Korean pine (Pinus koraiensis Sieb. et Zucc.) forest plantations and determined that the optimal rotation age that maximizes wood production and carbon capture was at the age of 86 years. In another study on oil palm [18], they applied the optimal control theory to model the dynamics of biomass growth and intrinsic biomass growth as state variables and considered felling as a control variable. The authors showed that the maximum oil production and carbon capture was reached at the age of 20 years. However, to our knowledge, no mathematical models have simultaneously modeled the relationship between the living biomass, the intrinsic biomass growth, the burned area, the reforestation, the felling and thinning, the fire prevention, and the relative humidity. Recently, [33] modeled the effects of the dynamics of living biomass, intrinsic growth, and burned area on carbon capture in forest plantations. The authors showed that biomass decreases in each cycle of regeneration because of forest fires, and suggested a strategy based on fire prevention in order to obtain maximum carbon capture. In this context, the objective of the present work is to determine optimal rotation strategies that maximize carbon sequestration in forest plantations. Based on the optimal control theory, a mathematical model is proposed to describe the dynamic relationship of carbon capture in forest plantations with control strategies such as reforestation, felling, fire prevention, and thinning, which are associated with state variables such as living biomass, intrinsic growth, and burned area. To verify the efficiency of the model, three scenarios are considered: realistic, pessimistic, and optimistic, using numerical methods to approximate its solution. In the case of the realistic scenario we tested with data of the species P. radiata.

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\left(t\right), r\left(t\right)$, and $I\left(t\right),$ 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\left(t\right)$, felling $F\left(t\right)$, fire prevention $S\left(t\right)$, and thinning $T\left(t\right)$. 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.

$$\dot{B}\left(t\right)=r\left(t\right)B\left(t\right)\left(1-\frac{B\left(t\right)}{K}\right)+\left[\beta R\left(t\right)\right]B\left(t\right)-\left[{\mu }_{1}I\left(t\right)+\sigma F\left(t\right)+\tau T\left(t\right)\right]B\left(t\right).$$

(1)

From Equation (1), $\beta$ is the rate at which biomass increases with respect to reforestation, ${\mu }_1$ is the rate at which biomass decreases due to fire effects, $\sigma$ is the rate at which biomass decreases due to felling effects, and $\tau$ 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

$$\dot{r}\left(t\right)=r_0-\rho r\left(t\right)+\nu T\left(t\right).$$

(2)

In Equation (2), $r_0$ represents the maximum individual growth rate under ideal conditions [36] and $\rho$ 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 $\nu$ 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, burned area decreases due to fire prevention, thinning, and relative humidity. Thus, the dynamics of the burned area are

$$\dot{I}\left(t\right)={\mu }_{2}I\left(t\right)\left(\frac{B\left(t\right)}{1+B\left(t\right)}\right)-\theta S\left(t\right)-\eta T\left(t\right)-hI\left(t\right).$$

(3)

From Equation (3), ${\mu }_2$ is the fire rate, $\theta$ is the fire prevention rate, $\eta$ 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

$$\{\begin{array}{c}\dot{B}\left(t\right)=r\left(t\right)B\left(t\right)\left(1-\frac{B\left(t\right)}{K}\right)+\left[\beta R\left(t\right)\right]B\left(t\right)-\left[{\mu }_{1}I\left(t\right)+\sigma F\left(t\right)+\tau T\left(t\right)\right]B\left(t\right) \hfill \\ \dot{r}\left(t\right)=r_0-\rho r\left(t\right)+\nu T\left(t\right) \hfill \\ \dot{I}\left(t\right)={\mu }_{2}I\left(t\right)\left(\frac{B\left(t\right)}{1+B\left(t\right)}\right)-\theta S\left(t\right)-\eta T\left(t\right)-hI\left(t\right), \hfill \end{array}$$

(4)

satisfying the initial conditions $B\left(t_0\right)=B_0, r\left(t_0\right)=r_0$, and $I\left(t_0\right)=I_0$. It is assumed that the amount of carbon stored in the biomass is proportional to the amount of biomass [37], therefore $C\left(t\right)=\alpha B\left(t\right)$ where $\alpha$ 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\left(t\right)\ge 0,r\left(t\right)\ge 0,I\left(t\right)\ge 0,R\left(t\right)\ge 0,F\left(t\right)\ge 0,S\left(t\right)\ge 0$ and $T\left(t\right)\ge 0$, the parameters $\left(\mathrm{K},\mathsf{\sigma },\mathrm{h},{\mu }_1,{\mu }_2,\tau ,\beta ,r_0,\rho ,\eta ,\theta ,\alpha ,\nu \right)\in 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

Table 1. Notation, definition, and units of each variable.

Notation	Definition	Unit
$B\left(t\right)$	Volume of living biomass	m3 ha−1
$r\left(t\right)$	Intrinsic growth of biomass	year−1
$I\left(t\right)$	Burned area per year	m2 year−1
$C\left(t\right)$	Carbon capture	Tn C ha−1year−1
$R\left(t\right)$	Forest reforestation	ha year−1
$F\left(t\right)$	Forest felling	ha year−1
$S\left(t\right)$	Fire prevention	US$m2 ha−1year−1
$T\left(t\right)$	Forest thinning	ha year−1

and Table 2, respectively.

The parameters, their units, definition, and notation are as follows.

Table 2. Notation, definition, and units for each parameter.

Notation	Definition	Unit
$\beta$	Rate of increase in biomass due to the effects of reforestation	ha−1
$h$	Relative humidity threshold to reduce fire	year−1
${\mu }_1$	Rate at which biomass decreases due to fire effects	m−2
${\mu }_2$	Fire parameter	year−1
$\sigma$	Rate at which biomass decreases due to felling effects	ha−1
$\tau$	Rate at which biomass decreases due to thinning effects	ha−1
$r_0$	Maximum growth rate	year−2
$\rho$	Natural mortality rate	year−1
$\nu$	Rate of increase in thinning over individual growth	ha−1
$\theta$	Fire prevention rate	ha US$−1year−1
$\eta$	Thinning rate	m2 ha−1 year−1
$\alpha$	Carbon capture rate	Tn C m−3 year−1

Table 2. Notation, definition, and units for each parameter.

Notation	Definition	Unit
$\beta$	Rate of increase in biomass due to the effects of reforestation	ha−1
$h$	Relative humidity threshold to reduce fire	year−1
${\mu }_1$	Rate at which biomass decreases due to fire effects	m−2
${\mu }_2$	Fire parameter	year−1
$\sigma$	Rate at which biomass decreases due to felling effects	ha−1
$\tau$	Rate at which biomass decreases due to thinning effects	ha−1
$r_0$	Maximum growth rate	year−2
$\rho$	Natural mortality rate	year−1
$\nu$	Rate of increase in thinning over individual growth	ha−1
$\theta$	Fire prevention rate	ha US$−1year−1
$\eta$	Thinning rate	m2 ha−1 year−1
$\alpha$	Carbon capture rate	Tn C m−3 year−1

2.2. The Optimal Control Problem

An optimal ecological control problem is presented with the objective to maximize the objective function $J$ representing the carbon capture in a fixed period $\left[0,T_f\right]$, such that

$$J\left(S,R,F,T\right)={{\displaystyle \int }}_0^{T_f}\left[{\omega }_{1}C\left(t\right)+{\omega }_2{S}^2\left(t\right)+{\omega }_3{R}^2\left(t\right)-{\omega }_4{F}^2\left(t\right)-{\omega }_5{T}^2\left(t\right)\right]dt,$$

(5)

subject to the state variables of the system (4). The state variables $B$, $r$, and $I$ are assumed free at the final time. The initial conditions are the real and adjusted values for the species P. radiata from studies carried out in central Chile:

$$\begin{array}{c}B\left(0\right)=B_0 \left(\mathrm{real}\right), r\left(0\right)=r_0\left(\mathrm{real}\right), I\left(0\right)=I_0\left(\mathrm{adjusted}\right)\\ B\left(T_f\right)=\mathrm{free}, r\left(T_f\right)=\mathrm{free}, I\left(I_f\right)=\mathrm{free},\end{array}$$

(6)

where the quantities ${\omega }_i, i=1, \cdots ,5$ are the weight parameters that balance the units of the terms of the objective function. It is considered $S^2\left(t\right)$, assuming that the contribution to carbon capture will be much higher if high prevention is applied and much lower if low prevention is applied. It is considered $R^2\left(t\right)$, assuming that the greater the reforestation, the greater the carbon capture. It is also considered $T^2\left(t\right)$, which implies that reducing the availability of fuel in the forest will contribute to decreasing the fire risk. Finally, it is considered $F^2\left(t\right)$, which increases the harvest rate to not compromise the owner’s timber production. On the other hand, quadratic terms are introduced in the Lagrangian framework to avoid the problem being linear, thus usual techniques can be applied [38].

To prove the existence of the optimal control problem we will follow the results of [39,40], which is demonstrated in detail in Appendix A. So, now we are in the condition to characterize the optimal control problem by means of Pontryagin’s maximum principle.

To solve the optimal control problem, Pontryagin’s maximum principle is used to characterize the optimal controls $S^{*}, R^{*}, F^{*},T^{*}$ [41]. To do this, it is necessary to determine the expressions for the adjoint variables using the Hamiltonian function [42].

$$\begin{array}{c}H\left(B,r,I,R,F,T,S,{\lambda }_1, {\lambda }_2,{\lambda }_3,t\right)=\left[{\omega }_{1}C+{\omega }_2{S}^2+{\omega }_3{R}^2-{\omega }_4{F}^2-{\omega }_5{T}^2\right]\\ +{\lambda }_1\left[rB\left(1-\frac{B}{K}\right)+\left[\beta R\right]B-\left[{\mu }_{1}I+\sigma F+\tau T\right]B\right] \\ +{\lambda }_2\left[r_0-\rho r+\nu T\right]+{\lambda }_3\left[{\mu }_{2}I\left(\frac{B}{1+B}\right)-\theta S-\eta T-hI\right].\end{array}$$

(7)

From Equation (7), the adjoint variables ${\lambda }_1,{\lambda }_2,$ and ${\lambda }_3$ satisfy the following adjoint ODE.

$$\begin{array}{c}\dot{{\lambda }_1}\left(t\right)=-\frac{\partial H}{\partial B}=-\alpha {\omega }_1-{\lambda }_1\left[r\left[1-\frac{2B}{K}\right]+\beta R-\left[{\mu }_{1}I+\sigma F+\tau T\right]\right]-{\lambda }_3\frac{{\mu }_{2}I}{{\left(1+B\right)}^2}\hfill \\ \dot{{\lambda }_2}\left(t\right)=-\frac{\partial H}{\partial r}=-{\lambda }_{1}B\left(1-\frac{B}{K}\right)+{\lambda }_2\rho \hfill \\ \dot{ {\lambda }_3}\left(t\right)=-\frac{\partial H}{\partial I}={\lambda }_1{\mu }_{1}B-{\lambda }_3\left[\frac{{\mu }_{2}B}{1+B}-h\right].\hfill \end{array}$$

(8)

Since the initial conditions for the state variables given in Equation (6) are free, the transversality conditions for the adjoint variables given in Equation (8) are:

$${\lambda }_1\left(T_f\right)=0, {\lambda }_2\left(T_f\right)=0, {\lambda }_3\left(T_f\right)=0.$$

(9)

Using Pontryagin’s maximum principle, the optimality of the control variables $S, R, F$ and $T$ establishes the characterization of the optimal controls $S^{*}, R^{*}, F^{*}$, and $T^{*}$ which satisfy the necessary first-order conditions

$$\begin{array}{c}{S}^{*}=min\left(max\left(0,\frac{{\lambda }_3\theta }{2{\omega }_2}\right),1\right)\\ \begin{array}{c} R^{*}=min\left(max\left(0,-\frac{{\lambda }_1\beta B}{2{\omega }_3}\right),1\right)\\ F^{*}=min\left(max\left(0,-\frac{{\lambda }_1\sigma B}{2{\omega }_4}\right),1\right)\\ T^{*}=min\left(max\left(0,-\frac{\left({\lambda }_1\tau B-{\lambda }_2\nu +{\lambda }_3\eta \right)}{2{\omega }_5}\right),1\right).\end{array}\end{array}$$

(10)

For solving optimal control problems three approaches are feasible: direct methods, indirect methods, and dynamic programming [43,44]. We used the indirect methods because they are more robust than the base on the classical theory of Pontryagin’s maximum principle that reduces the optimal control problem to the solution of a boundary value problem. We used the fourth-order Runge–Kutta scheme using Octave/MATLAB software [44,45]. In the following section, we performed numerical simulations to approximate the solutions to the optimal control problem. We simulate using the fourth-order forward Runge–Kutta method [38]. This iterative method consists of solving the controlled system (4) using a fourth-order forward Runge–Kutta scheme and the traversal or terminal conditions given in Equation (9) in a time interval $\left[0,T_f\right]$. Then, the adjoint system (8) is solved by a fourth-order backward Runge–Kutta scheme using the solution of the current iteration of the controlled system (4). The characterization of the controls of the Equation (10) is updated through a convex combination of the above controls. The procedure stalls if the values of the variables of the previous iteration are very close to the present iteration [44].

3. Results and Discussion

The main objective of this work is to determine the optimal strategies of reforestation $R\left(t\right)$, felling $F\left(t\right)$, fire prevention $S\left(t\right)$, and thinning $T\left(t\right)$ that maximize carbon capture in fast-growing forest plantations. To start with, the final time is set equal to $T_f=200$ years. The initial conditions assumed for the biomass of P. radiata plantations are $B\left(0\right)=7.5$ m³ ha⁻¹ and, for the intrinsic growth, a maximum growth $r\left(0\right)=0.0725$ year⁻¹, given in [46], is considered. For the burned area, an average $I\left(0\right)=0.07$ m² year⁻¹ is adjusted. Then, the Runge-Kutta solver of order four is run for the following three scenarios:

Realistic Scenario. Real parameters from

Table 3. Real parameters for P. radiata.

$\mathit{K}$	${\mathit{\mu }}_1$	$\mathit{\beta }$	$\mathit{\sigma }$	$\mathit{\tau }$	$\mathit{\rho }$	${\mathit{r}}_0$	$\mathit{\nu }$	${\mathit{\mu }}_2$	$\mathit{\theta }$	$\mathit{\eta }$	$\mathit{h}$	$\mathit{\alpha }$
400	0.055	2 × 10−5	2 × 10−5	1 × 10−5	0.06	0.07275	0.159	0.4097	0.15	9.3 × 10−7	0.135	0.5

(see Appendix B) were considered for the species P. radiata from studies carried out in central Chile.

Pessimistic Scenario. According to the IPCC, in recent decades global warming has increased GHG emissions and their accumulation is causing irreversible climatic changes. As a result, the Earth is experiencing an increase in temperature, which will lead to a sustained decrease in ambient relative humidity and a prolongation of the duration and frequency of drought events [47] and fires [48]. In this context, this scenario considers a decrease in ambient relative humidity, which will have a negative effect on biomass growth, and will also favor conditions for an increase in the frequency and intensity of fires which, in turn, will increase the burned area [49]. This scenario also considers a decrease in the budget dedicated to fire prevention, due to the effects that the global health and war situation will have on local and global economies [50,51]. For this scenario, the fixed values of the parameters shown in

Table 4. Pessimistic parameters for forest plantations.

$\mathit{K}$	${\mathit{\mu }}_1$	$\mathit{\beta }$	$\mathit{\sigma }$	$\mathit{\tau }$	$\mathit{\rho }$	${\mathit{r}}_0$	$\mathit{\nu }$	${\mathit{\mu }}_2$	$\mathit{\theta }$	$\mathit{\eta }$	$\mathit{h}$	$\mathit{\alpha }$
400	0.02	0.1	0.09	1 × 10−5	0.06	0.061	0.03	0.47	0.052	9.3 × 10−7	0.11	0.5

, which correspond to an artificial data set, are considered.

Optimistic Scenario. This scenario assumes that there will be a substantial increase in the annual rate of afforestation and reforestation, due to the growing concern to prevent the global average temperature from rising over the next century and the growth of new commodities based on environmental values, such as carbon capture [5,52,53]. The Paris Agreement and the Kyoto Protocol have designed economic instruments that provide financial incentives for governments to protect the environment using the private sector for strict pollution standards [29]. In this scenario, most of the parameters of the realistic scenario in Table 3 were considered, assuming an afforestation and reforestation rate higher than that of the realistic scenario $\left(\beta =0.04\right)$. By increasing the rate of afforestation and reforestation, carbon capture will increase, decreasing the greenhouse effect and regulating environmental temperature and humidity patterns. According to [54,55], they showed that trees have a positive effect on relative humidity increase and temperature reduction. Thus, this scenario assumes an increase in ambient relative humidity in the area where P. radiata plantations are concentrated ($h=0.14$). On the other hand, the presence of financial incentives would increase the fire prevention budget by 120% ($\theta =0.18$) in activities such as thinning, which will indirectly contribute positively to individual growth in the long term ($\nu =0.42$) and decrease the continuity of fuel susceptible to fire.

The optimal trajectories for live biomass are similar for the three scenarios (Figure 2a), however, as expected, in the optimistic scenario there is greater biomass accumulation (the total area under the curve is 12 300) compared to the other scenarios (10 739 in the realistic and 5 125.4 in the pessimistic scenarios, respectively). In general, the minimum biomass volume coincides with the maximum burned area (Figure 2c), although in the pessimistic scenario this occurs at earlier ages as a result of lower ambient humidity and higher fire propensity [56], which shortens the rotation age. It is also observed that in the following plantation rotation cycles the biomass is lower than in the first cycle for the three scenarios. It could be due to the presence of fires in the three scenarios which causes the volatilization of the main soil nutrients [57]. Volatilization of some major soil nutrients, such as nitrogen and phosphorus, can affect tree growth and limit terrestrial carbon sequestration [58,59]. This model does not consider artificial fertilization of the soil and the plantation starts growing when the minimum burned area occurs. In [10], they argue that a higher reforestation rate, together with prolonging the rotation age, are key strategies to maximize carbon capture and mitigate the negative effects of global climate change. The model corroborates the above since the optimistic scenario considers a higher rate of afforestation and reforestation than the realistic scenario, which increases the volume of forest biomass and there is a prolongation of the rotation age to 29 years, as opposed to the earlier rotation ages determined by the pessimistic (23 years) and realistic (27 years) scenarios. The rotation age of the realistic scenario is within the rotation age range reported by [60] in operational plantations of P. radiata in central Chile. In the realistic scenario the model maintains the same rotation age as in [60], while for the pessimistic scenario the rotation age was reduced since the forest was affected by the fire, which forces early felling to avoid damage by new fires, with the consequent emission of CO₂ to the atmosphere. In the case of the optimistic scenario, the model prolongs the rotation age due to the positive impact of market incentives for environmental protection.

Carbon capture and burned area follow the same trend in all three scenarios (Figure 2b,c). As expected, the maximum burned area occurs years after the maximum biomass volume is produced. The realistic scenario shows a small burned area under the curve of 6258.8, which suggests that the model realistically reflects the current fire prevention and firefighting situation in Chile, which is more efficient [61,62]. It is also observed in the optimistic scenario that the burned area under the curve is 8313.2, which is greater than the realistic scenario. That scenario happens because in the optimistic scenario there is a greater volume of biomass (Figure 2a). However, in the pessimistic scenario, despite a low volume of biomass, there is a greater burned area under the curve of 10,974 compared to the other two scenarios. This situation is because the relative humidity threshold is the lowest of the three scenarios and the forest is more prone to burning, which decreases the biomass.

Finally, Figure 2d shows that the intrinsic growth variable is higher in the optimistic scenario, while the pessimistic scenario is lower and with a flat trend due to the negative effect of the higher number of fires, which affects individual plantation growth.

Figure 3 shows that the carbon accumulated in 200 years is 6050, 5240, and 2560 Tn C for the optimistic, realistic, and pessimistic scenarios, respectively. Overall, the optimistic scenario is 57.69% higher than the pessimistic scenario and 13.4% higher than the realistic scenario.

As abovementioned, the optimal control theory was applied in [18,29] by independently assessing incentives that promote reforestation, the use of clean technologies, or other factors that modify the rotation age as control strategies to maximize carbon capture. However, in the present work, optimal control theory is applied with four strategies accordingly in three simulation scenarios (Figure 4). The results indicate that in the realistic and optimistic scenarios the felling strategy considers an optimal rotation age of 29 and 27 years, and no maximum felling effort is applied, while in the pessimistic scenario the optimal rotation age is shorter and maximum felling effort should be applied throughout the plantation life cycle, and the maximum effort has a limit value of one (Equation (10)). The previous facts are due to the presence of fire risk, which decreases the optimal rotation age [20]. Then, as fires increase, short rotations could be used and maximum felling effort applied as in the pessimistic scenario, so as not to increase the release of CO₂ to the atmosphere and not compromise the forest owner business. As previously mentioned, in the case of the optimistic and realistic scenarios, no maximum felling effort is applied. This coincides with the Chilean national reality since, despite the increase in demand for timber in the country in the last 15 years, the felling rate has dropped from 55,000 ha year⁻¹ to 25,000 ha year⁻¹ [63,64]. While this is positive from an environmental point of view, as more planted hectares are being maintained and contribute to mitigating climate change, it could lead to a shortage of wood for processing plants, which would increase the pressure to shorten the rotation age. However, this hypothesis needs further analysis. According to the IPCC, a significant percentage of GHGs (30%–40%) can be reduced by avoiding felling, forest degradation, and the recovery of forest areas.

On the other hand, in all three scenarios, the reforestation strategy is applied when felling begins. The model shows that maximum reforestation effort should be applied for the pessimistic and optimistic scenarios, while in the realistic scenario no maximum reforestation effort is applied. Thus, the model reports that today 100% of the reforestation rate is not reached, despite the Paris Agreement and the Kyoto Protocol that offer financial incentives to governments to increase forest area. According to Chilean Decree Law 701, between 1998 and 2015 36% of the area was forested, 39% of the area corresponds to reforestation, and more than 20,000 hectares are subsidized for forest management, which meant an investment of 388 million USD for the State in that period. According to [29] they suggest that it is more efficient to invest in reforestation than other practices such as clean technologies, since forested areas can contribute to climate change mitigation in a more efficient and environment friendly way. However, even though maximum reforestation effort is applied in the pessimistic scenario, the biomass volume is lower than in the other two scenarios (Figure 2a), which corroborates the negative effect of the larger burned area on biomass growth.

In the fire prevention strategy, the model suggests that in the three scenarios presented, the maximum prevention effort should be applied years after reaching the maximum biomass volume since the plantation in its adult state forms a continuous mass of fuel that is more difficult to control in case of fire. Therefore, the model based on this strategy, considering all plantations, suggests investing in fire prevention on average four years before and four years after felling (Figure 2c). Prevention applied before helps to minimize the fire risk in plantations ready to be felled, whereas prevention applied after helps to manage the forest residues in the soil (i.e., fuel) by felling operations. On the other hand, in the pessimistic scenario, it is observed that there is a decrease in the prevention strategy of the following rotation cycles because the biomass decreases (Figure 2a). Therefore, the model communicates that the owner should invest according to the amount of existing biomass.

Finally, for the thinning strategy in the three scenarios, it is observed that the maximum thinning effort should be applied from the initial moment to the moment of fire prevention. It is also observed that in the case of the pessimistic scenario the thinning periods are longer because the forest is more prone to burning.

In general, the optimal rotation age for the realistic scenario in the following rotation cycles turned out to be 29 years, whereas for the optimistic scenario in the following rotation cycles this is shortened to 24 years; this is because the presence of fires shortens the rotation age. Even though the rotation age is shortened there is a greater biomass and therefore greater carbon capture than in the other two scenarios (Figure 2a,b). However, in the pessimistic scenario in the following cycles, the rotation age is prolonged to ca. 45 years, because forests need a long time to recover from the fire. Therefore, the rotation age should be applied at the maximum time that the burned area reaches.

4. Conclusions

In this study, we determined optimal management strategies that maximize carbon capture in fast-growing forest plantations using the optimal control theory. The model effectively simulates the optimal dynamics of live biomass, intrinsic growth, and burned area to consider four strategies such as reforestation, felling, fire prevention, and thinning. To evaluate the effectiveness of the model, three scenarios have been considered: realistic, pessimistic, and optimistic. The parameters for modeling the realistic scenario were based on real data for the exotic species P. radiata.

The model predicts a higher biomass volume for the optimistic scenario, with a consequent higher carbon capture than in the other two scenarios. The optimal solution for the felling strategy suggests that, in order to increase carbon capture, the rotation age should be extended and the felling rate decreased. Likewise, the model corroborates that reforestation should be carried out immediately after felling, applying maximum reforestation effort in the optimistic and pessimistic scenarios. It is suggested that in order to increase CO₂ capture, large forestry companies and government agencies should increase investment in afforestation and reforestation. On the other hand, the model indicates that, although maximum prevention effort should be applied during the life cycle of the plantation, it should be proportional to the volume of biomass. Finally, the optimal solution for the thinning strategy indicates that in all three scenarios maximum thinning effort should be applied until the time when the fire prevention strategy is applied. Here, the optimistic scenario is considered as a possible alternative in forestry activities to maximize carbon capture. While this scenario yields the largest carbon capture, its implementation requires joint efforts between forest companies and the government to prolong the rotation age and economically incentivize reforestation, respectively.