We developed a forest planning model integrating two operational scales (single-stand and forest levels) for the optimization of timber production and carbon sequestration in forest teak (
Teak (
Despite its long history as a planted species, reliable information regarding growth and productivity of teak is lacking, and gaps still exist regarding short rotations, growth dynamics and thinning regimes (
Different algorithms use different strategies for seeking the optimal solution through the solution space and various criteria allow to evaluate the quality of the solutions obtained. The maximization of profitability and monetary benefits is usually the optimal solution in industrial forest plantations. However, optimizing for carbon sequestration has lately become a major goal of many forest planning models, due to the significant role of forests in the reduction of atmospheric CO2 and climate change mitigation. The adoption of the Kyoto protocol and the creation of carbon markets stimulated the scientific community to focus on the assessment of the CO2 amount stored by forest ecosystems, thus bringing this issue to the attention of forest managers. Indeed, carbon uptake is currently considered a new goal of silvicultural management (
Many forest planning models that incorporate carbon sequestration as optimization criteria have been reported in the literature (
The aim of this work was to optimize management plans of multiple-stand teak plantations by developing a heuristic model which simultaneously maximizes financial benefits from timber production and carbon sequestration.
The model simulates a teak plantation project comprising many stands differing in age, site quality, and initial spacing. Time of thinning and final harvest for each stand depend on the trade-off between wood production and carbon storage to maximize financial benefits for the project over a planning period of 40 years. Two optimization models working at different planning levels were developed. The first one (sub-model 1) operates at the stand level, and provides a set of “best thinning schedules” for each stand, assuming a fixed final harvest after 30 years. The second model (sub-model 2) operates at the forest level, by choosing for each stand the “best schedules” that maximize financial gains for the whole project, taking into account both carbon and timber goals subject to various constraints. Although the rotation age initially assigned to all stand is 30 years, the forest-level sub-model can modify the optimal harvest age of each stand (
Sub-model 1 generates a set of thinning regimes for each stand with different initial spacing and site quality, and selects a subset of eight “best” thinning schedules, including the “no thinning” option. The objective function maximizes the financial gain, considering at the same time the merchantable timber production and the carbon sequestration in terms of Net Present Value (NPV, in US$). For each stand, several simulations are done taking into account the desired number of thinnings. The best schedules for each stand are the inputs for sub-model 2.
Sub-model 2 searches for the optimum harvest schedule (Maximum NPV) for the entire plantation, indicating which stands to thin or harvest each year, and considering both timber production and carbon sequestration. The objective function takes into account the incomes from timber and carbon capture, the cost of planting, maintenance, and transport of timber products. Constraints include a yearly timber production goal varying from year to year. In order to consider the additional costs due to spreading of harvest activities (moving personnel and machinery across the planted area), a spatial index seeks for minimizing the displacement of equipment within the plantation, and assumes no spatial adjacency constraints.
A whole-stand growth and yield model based on differential equations simulates the stand production by projecting stand variables such as basal area, stand density, average diameter, height, and volume for various combinations of site quality, initial spacing, and thinning regime. Basal area growth is calculated as (
where
Top height growth was derived from a site index (base age 16) family of anamorphic curves derived from the Schumacher-Hall model (
Model parameters were estimated based on a data set from trials and permanent/ temporal plots in plantations ranging from 2 to 32-yrs-old and located in the Western Plains of Venezuela (Tab. S1 and Tab. S2 in Supplementary material). Trials and plots were established on adequate climatic conditions for teak,
where
For each year of the planning horizon, the net carbon balance (
To estimate the stem biomass, the outputs of the growth and yield model were taken as inputs for the sub-model 2. Woody biomass dry weight was assessed from stem timber volume based on the basic density of teak wood (0.54 Mg m-3 -
The estimation of carbon emissions from the decomposition of harvest residues (such as roots, branches, stumps, and dead wood) was included in sub-model 2. The rates of carbon emissions generated by the forest products after harvesting or thinning were projected according to product classes: (i) long (i.e. large solid wood products), (ii) medium, and (iii) short duration. According to
where
where
Eqn. 6 implies that the biomass in each category decomposes as a constant proportion of the remaining biomass. For each use category, emission and decomposition will stop when a very tiny fraction (<0.01 MgC ha-1) is emitted. Thus, a small fraction of biomass will not further decompose,
Seven use categories were considered in the model: (1) roots; (2) deadwood; (3) branches and stumps; (4) bark and debris; (5) short duration wood products (small poles for fencing); (6) midterm duration products (struts and large poles); and (7) long term duration products (sawn timber). Decomposition time (DT) and anthropogenic time (AT) were taken from
We developed a restricted optimization model using integer variables to represent ages of thinning, and continuous variables for thinning intensities (percentages of the stand basal area).
The objective function was (
where
The model constraints were (
where
Constraint in
A binary variable constrained optimization was modeled to achieve the optimum harvest plan. The objective function was (
where
The net present value of merchantable timber was calculated as follows (
where
The net present value of carbon was calculated as follows (
where
The spatial index (
with i ≠ w, where
where
The decision variables (
Management regime was defined based on initial spacing, thinning regime, and the final harvest age. Thinning ages and intensities were determined by sub-model 1 depending on the planting year of stands and the age of last thinning, which must be done at least three years before the final harvest. The constraints were (
with
Constraint in
The existence of nonlinear relationships among integer and continuous decision variables, and the model structure make extremely difficult to solve the above equations with classic mathematical techniques; therefore, heuristic techniques were used.
For solving sub-model 1, we used Genetic Algorithms (GA), a robust and efficient technique successfully applied to find optimal or close-to-optimal solutions when exact methods are difficult to implement (
An optimization algorithm based on Simulated Annealing (SA) was incorporated in sub-model 2 in order to search for the best management plans at the forest level. The SA technique has been previously applied to forest harvest planning (
We used the data set obtained from teak plantations in Venezuela (described above) to simulate a teak plantation covering 10.000 ha and comprising 200 stands (50 ha each). Simulated stands differed in initial spacing, with 100 stands having an initial density (Di) of 1600 trees ha-1 (2.5×2.5 m) and the remaining 100 stands having Di = 1111 trees ha-1 (3.0×3.0 m). Moreover, the simulated stands differed in site quality: for each type of initial spacing, 50 stands were assigned to site quality I (SQ I), with a site index of 27 m (base age 16 years) and maximum basal area (
The planning project spanned 40 years, from initial planting of the first stand until the final harvest of the last stand. Stand planting took place in the first ten years of the project (0 to 10). The timber goal was to obtain merchantable wood from logs ≥ 20 cm in diameter. At the end of the planning period (40 years), all stands should have been harvested. We assumed that stands were not replanted after harvest, thus allowing for secondary vegetation to recover and some teak trees to re-sprout; however, the additional carbon sequestration of these abandoned stands was not considered in the project.
A different associated transport cost (in US$ m-3) was assigned to each stand, depending on its distance from the sawmill. Roundwood prices were assigned according to diameter classes (between 53 and 400 US$ m-3 - see Tab. S6 in Supplementary material), and base carbon price was 10 US$ Mg-1 of C, as commonly used in financial analyses (
Three alternative management scenarios were set with the objective of maximizing: (A) total NPV for merchantable timber production only; (B) total NPV for both timber and carbon sequestration; (C) total NPV for carbon sequestration only.
We carried out a three-stage sensitivity analysis for scenario B only (Timber production + C sequestration). The stage 1 consisted of a sensitivity analysis for individual stands (sub-model 1). The magnitude of changes in the optimal solution was analyzed considering changes in harvest costs, interest rates, growth rates, timber prices, and carbon prices. Sensitivity was tested for: (i) growth rate (
At stage 2 we carried out a sensitivity analysis for changes at forest level considering the sensitivities obtained in stage 1 (stand level), as changes in the optimal solution at stand level propagate to the forest level changing the optimum harvest plan. If no changes occurred at the stand level for a given variable, then the sensitivity analysis was carried out at the forest level for that variable. Finally, in stage 3 we carried out an additional sensitivity analysis at forest level, but considering only forest-level input variables such as transport costs from the stands to external roads and production quotas. Transport costs were changed from the base costs up to ±50%, at 10% intervals for each of the four categories. We tested the model sensitivity to modifications of the annual production quota by raising or lowering the quota in certain periods (Tab. S7 in Supplementary material).
The above validation approach was carried out at both stand and forest levels for the scenario of simultaneously optimizing timber production and carbon sequestration (scenario B). Fifty runs at each planning level were done and the optimum NPV and its basic statistics were recorded.
The developed model simulated the best thinning regimes for each stand according to initial spacing and site quality under the optimization scenarios. For those scenarios that included timber production, either alone (scenario A) or combined with carbon sequestration (scenarios B), the best thinning regimes were the same, but they differed only in the magnitude of NPVtotal, which was larger when both objectives were considered. On the other hand, when NPVcarbon was the only optimization/maximization goal (scenario C), thinnings were prescribed at later ages and lower intensities. For example, for scenario A (Di = 1600 trees ha-1, SQ = I, harvest age = 30), the best thinning regime for scenarios A and B consisted of three thinnings at ages 5, 10, and 19, removing 46.5%, 47.7%, and 29% of the total basal area, respectively (
We ran the model to obtain the optimal harvest sequence for 200 stands during the planning period under the proposed scenarios. For each year, the model indicated which stands should be thinned or harvested to maximize the expected benefits. Again, the simulated NPVtotal value obtained for the scenario B (both optimization criteria) was higher than those obtained for either scenario A (maximizing timber production) or scenario C (maximizing carbon sequestration) alone. In addition, there were no changes in the optimal harvest sequence between scenarios A and B.
The average rotation age obtained by optimizing for NPVwood was shorter than that obtained by maximizing for NPVcarbon only (
For the base scenarios (
The developed model provided estimates of the carbon stored in standing trees, died trees, forest products, harvest residues, and woody debris for each year of the planning period (
In our model, at the year 40 all stands not harvested during the project life were cut (76.879 m3 in the A-B scenarios, and 273.600 m3 in scenario C), thus the largest part of the carbon previously sequestered was kept as merchantable products and as debris from harvest and wood processing. In scenario C (max NPVcarbon), the amount of carbon stored remained larger until year 69, because the final harvest was much larger than in scenario A, but as short to medium duration products. On the other hand, in scenario A (max NPVwood), the total amount of C stored as harvested products was lower, but as long duration products. In both scenarios, the remaining fraction of carbon returned to the atmosphere, according to decomposition and anthropogenic times for products and debris; hence, the curve for scenario C decreased gradually (
The model was very sensitive to changes in growth rates of basal area and in timber prices, as very small changes generated different optimum harvest plans. In addition, the model was sensitive to changes in interest rates by changing the optimal solutions in terms of thinning regimes and harvest ages (
Changes in allowable production quotas for the whole project modified the optimal harvest plan. When annual quotas were raised from five to 10.000 m3 between the years 7 and 20 (scenarios 1, 2, 3 - Tab. S7 in Supplementary material), the optimum harvest plan changed because thinning regimes and final harvest for each stand changed to meet the new quotas. When the quota for years 21-40 rose from 60 to 65.000 m3, an optimal solution could not be found. For scenarios where production quotas were reduced (scenarios 5, 6, 7, 8 - Tab. S7 in Supplementary material), the optimal solutions were the same as the base scenario; yet, volume yield per year surpassed the assigned annual quotas. Changes in harvest and transport costs did not change the optimal solution, thus the proposed model is not sensitive within the range of changes assumed for these variables.
To validate the model, descriptive statistics for a sample of 50 solutions for the thinning regime sub-model (stand level) were produced considering one to four thinnings (
At forest level, model validation showed that the difference between the highest and lowest NPV was 4.68 % (
The model proposed in this study simulates a full teak plantation project at the stand and forest levels and provides a optimized harvest plan by identifying which stands to thin or harvest each year in order to both maximize the benefits and meet the operating constraints. The simulation carried out in this study (200 stands, 10.000 ha) was solved by the model in a few seconds. The results obtained were satisfactory, considering the model complexity, assumptions and constraints. For example, for maximizing NPVwood, the model increased the number of thinnings and their intensities to favor larger diameter growth and logs with higher value as sawn timber. Moreover, the simulated dynamics of the stand variables (diameter, basal area, stand density) were within the expected limits for acceptable quality teak stands in terms of both growth and yield (
The proposed model can be usefully applied to optimize harvest plans in teak plantations by maximizing wood production and/or carbon sequestration as separate objectives or simultaneously. According to the model results, for a price of 10 US$ Mg-1 of C and the current timber prices (53-400 US$ m-3, depending on log diameter) there were no differences in cutting patterns (best thinning regimes and optimal harvest plan) between the scenario optimizing timber production (scenario A) and that maximizing both wood production and C sequestration (scenario B). This is because timber prices were much higher than the market carbon prices. The lack of differences between the above scenarios can be explained because the model is not sensitive to changes in C prices between 0 and 40 US$ Mg-1 (
When maximization of NPVcarbon was set as the only management goal in the model, the generated harvest plans had lesser, later, and lower intensity thinnings, as well as longer rotations.
The simulations carried out revealed that maximizing for NPVcarbon C stock in the plantation was larger until year 69 as compared to scenarios maximizing NPVwood; however, for the latter scenarios a larger fraction of C remained stored in long duration products until the year 116.
Our results showed that implementing management plans to favor carbon sequestration reduce the economic benefits from timber production, and
Our results suggest that teak has a good potential for carbon storage due to its fast growth and wood durability. Additional benefits could be obtained from C sequestration only by optimizing the NPVwood, as long as C prices are below 40 US$ Mg-1. Managing teak projects as infinite rotation plantations would transform them in permanent carbon sinks, also increasing their capacity for soil carbon storage.
Most teak plantations, especially in Latin America, are small to medium size private-owned projects (
Although our model is calibrated on teak trees growing in the Venezuelan Western Plains (alluvial soils), growth and yield performances of teak stands are comparable with those observed in other countries, like Brazil, Mexico, Colombia, Ecuador, and other states of Central America, in terms of growth rates, mean annual increments (8-20 m3 ha-1 yr-1), yield, and financial results for a given site index. Moreover, volume and yield tables from Venezuela have been successfully extrapolated to other countries. However, we caution to extrapolate our results to Asian or African Countries, or to clonal plantations under intensive management. Furthermore, as the growth and yield model used is based on the basal area, its calibration to more site specific conditions and management regimes can be done with no risk of obtaining unrealistic results. Finally, the modular structure of the proposed model allows for changes in the growth and yield model without affecting the other model components.
So far, the high prices of high quality teak timber has reduced the need of harvest optimization techniques for increasing efficiency in timber production and planting operations. As many teak plantations are reaching commercial ages in Latin America, timber supply could exceed the demand, thus lowering timber prices. Further, certification of sustainability is increasingly required on the market for timber products from plantation projects, which implies many environmental requirements. In addition, many projects receive subsidies in their initial phases and are subject to enter the incipient carbon market. Thus, optimizing carbon sequestration and timber production in teak plantations is going to be crucial in the upcoming years.
In this study, a model was developed for optimizing both thinning regimes in individual teak stands and the harvest plan for the whole plantation, aimed at simultaneously maximizing timber production and carbon sequestration. The adopted modeling approach allowed the integration of both stand and forest planning levels. At the stand level, we approached the stocking-regime optimization problem using a heuristic model based on Genetic Algorithms. Production curves were generated according to stand characteristics by a forest stand growth and yield simulator. At forest level, we used a simulated annealing heuristic technique to determine optimal harvest ages and assign an appropriate thinning regime to each stand matching the annual production quotas.
The proposed harvest planning model may help managing small to large scale teak plantations, by providing information on plantation growth and yield, and considering the impact of stocking and harvest regimes. The possibility of quickly running and examining several potential harvest plans under various scenarios makes the model useful as a starting point for designing management plans, fieldwork task assignments, road construction, and equipment displacement.
This work was supported by the
Basic structure of the forest planning model.
Dynamics of net carbon storage (from total accumulated biomass). Average life cycle of carbon stored as wood in standing trees, forest products, and harvest residues for a 10.000 ha teak plantation under the base optimization scenarios (maximum NPV criteria). Carbon price 10 US$ Mg-1of C and timber prices between 53 and 400 US$ m-3 according to log size categories (Tab. S6 in Supplementary material). (Scenario 1): maximizing NPV for carbon sequestration (green dashed line); (scenario 2): maximizing NPV of wood production only (orange dashed line); (scenario 3): Maximizing NPV of wood and carbon simultaneously (brown line).
Dynamics of net carbon storage under scenarios of wood production only (dark green dashed line), wood + carbon sequestration (40 US$ Mg-1 - brown solid line); wood + carbon sequestration (100 US$ Mg-1 - orange dotted line), wood + carbon sequestration (150 US$ Mg-1 - green dashed line). Timber prices are reported in Tab. S6 (Supplementary material).
Net present value (NPV) of the best thinning regimes for teak stands (initial density: 1600 trees ha-1; Site quality: I; Rotation age: 30) for scenarios: A (wood production only); B (wood production and carbon sequestration); C (carbon sequestration only). (NT): number of thinnings for each regime. For each thinning, the age (years) and the percentage of removed basal area (in parentheses) are reported.
Scenario | NT | Thinnings | NPV(US$ ha-1) | |||
---|---|---|---|---|---|---|
First | Second | Third | Fourth | |||
A/B | 0 | - | - | - | - | A: 7100; B: 1390 |
A/B | 1 | 9 (79.9) | - | - | - | A: 6290; B: 6780 |
A/B | 2 | 5 (55.8) | 14 (53.7) | - | - | A: 8980; B: 9470 |
A/B | 3 | 5 (46.5) | 10 (47.7) | 19 (29.0) | - | A: 9090; B: 9580 |
A/B | 4 | 8 (44.6) | 16 (25.6) | 20 (47.7) | 25 (26.7) | A: 7900; B: 8470 |
C | 0 | - | - | - | - | 680 |
C | 1 | 27 (27.9) | - | - | - | 672 |
C | 2 | 23 (25.2) | 27 (25.7) | - | - | 670 |
C | 3 | 20 (26.2) | 24 (27.4) | 27 (48.1) | - | 666 |
C | 4 | 17 (26.0) | 21 (25.4) | 24 (36.3) | 27 (39.0) | 651 |
Rotation age for stands under optimal planning schedules for each management scenario: (A) Wood production; (B) Wood production + carbon sequestration; (C): Carbon sequestration.
Scenario | Clearcut age (years) | Clearcut stands(age < 25 years) | Clearcut stands(age ≥ 25 years) | ||||
---|---|---|---|---|---|---|---|
Average | Min | Max | N | % | N | % | |
A and B | 27.2 | 20 | 38 | 61 | 30.5 | 139 | 69.5 |
C | 30.5 | 20 | 39 | 25 | 12.5 | 175 | 87.5 |
Number of stands with a given number of prescribed thinnings under the optimization scenarios: (A) wood production only; (B) wood production + carbon sequestration; (C): carbon sequestration only.
Scenario | Number of thinnings | Total | ||||
---|---|---|---|---|---|---|
0 | 1 | 2 | 3 | 4 | ||
A and B | 29 | 39 | 55 | 55 | 22 | 200 |
C | 62 | 39 | 38 | 31 | 30 | 200 |
Net present value for wood (NPVwood) and carbon sequestration (NPVcarbon) according to optimization objectives.
NPV | Optimization objectives | Difference(US$ million) | |
---|---|---|---|
Woodproduction | Carbonsequestration | ||
NPVwood (US$ million) | 31.03 | 8.08 | -22.95 |
NPVcarbon (US$ million) | 8.47 | 11.94 | +3.47 |
Total NPV | 39.50 | 20.04 | - |
Sensitivity analysis for selected variables. Data indicate the ranges in which the optimal solution (harvest patterns) is modified with respect to the base solution.
Variable / Parameter | Model sensitive tochanges (range) | Model not sensitiveto changes (range) |
---|---|---|
Basal area growth rate | [0.10, 0.25] | None |
Harvest cost | None | [7.12, 21.36] US$ m-3 |
Interest rate | > 10% | < 10% |
Wood prices | All scenarios | No scenarios |
Carbon price | ≥ 40 US$ Mg-1 of C | [0, 40] US$ Mg-1 C |
Transport cost | None | All scenarios |
Production quotas | Scenarios 1-4 | Scenarios 5-8 |
Model validation. Net present value (NPV) statistics from 50 runs for the stand-level model and for the forest-level model. Study case: 10.000 ha teak plantation project. Statistics for NPV are: average, minimum (Min), maximum (Max), standard deviation (SD) in thousand US$; and coefficient of variation in percentage (CV); (D1): percentage difference between the maximum value and the minimum value; (D2): percentage difference between the average and the maximum value; (NT): number of thinnings.
Statistics | Stand level | Forestlevel | |||
---|---|---|---|---|---|
NT = 1 | NT = 2 | NT = 3 | NT = 4 | ||
Average | 6.55 | 8.83 | 8.70 | 7.12 | 38185.53 |
Min | 6.32 | 8.10 | 7.29 | 6.13 | 39500.00 |
Max | 6.78 | 9.47 | 9.58 | 8.47 | 37649.18 |
SD | 0.20 | 0.37 | 0.47 | 0.52 | 384.99 |
CV | 3.12 | 4.18 | 5.40 | 7.23 | 1.01 |
D1 | 6.87 | 14.48 | 23.95 | 27.67 | 4.68 |
D2 | 3.53 | 6.76 | 9.14 | 15.90 | 3.33 |
Tab. S1 - Growth and yield data from
Tab. S2 - Parameter estimates for the growth and yield model.
Tab. S3 - Parameter estimates for calculating carbon storage in biomass.
Tab. S4 - Decomposition and anthropogen times for the woody biomass of teak trees (adapted from
Tab. S5 - Proportion of product type by diameter categories of teak trees.
Tab. S6 - Roundwood prices of teak timber from Latin American plantations (2012).
Tab. S7 - Scenarios for determining the model’s sensitivity to changes in required annual production quotas for the whole project.