Optimizing thinnings for timber production and carbon sequestration in planted teak (Tectona grandis L.f.) stands

Aim of study: We developed an optimization model for determining thinning schedules in planted teak (Tectona grandis L.f.) stands that maximize the financial output in terms of soil expectation value (SEV) and net present value (NPV) considering a) the simultaneous optimization of timber production and carbon (C) sequestration and b) only for C sequestration. Area of study: Planted teak forests in the western alluvial plains of Venezuela. Material and methods: We integrated a stand growth and yield model with a constrained optimization model based on genetic algorithms (GA) for determining optimal thinning schedules (number, age, and removal intensity) that maximize SEV when simultaneously managing for timber production and C sequestration. The data came from permanent plots established in planted teak stands with remeasurements from 2 to 32 yr.-old. Plots differ in site quality, initial spacing, and thinning schedules. We obtained optimal thinning schedules for several scenarios combining site quality, initial spacing, interest rates, harvest and transport costs, as well as timber and C prices. The stand growth and yield model estimates timber products and C flows (storage and emissions) until most stored C is reemitted to the atmosphere. Main results: When considering simultaneously both, timber production and C sequestration, the scenario with the maximum SEV consisted of initial stand densities = 1,111 trees ha-1, site quality (SQ) I, harvest age 20 years, and four thinnings (ages 6, 10, 14, 17 with removal intensities 26 %, 28 %, 39 %, and 25 % of stand basal area respectively). For maximizing C sequestration only, the best schedule consisted of 1,600 trees ha-1, SQ I, harvest age 25 years, with no-thinning. A sensitivity analysis showed that optimal schedules and SEV were highly sensitive to changes in interest rates, growth rates, and timber prices. Research highlights: • The management schedules favoring merchantable timber production are not the same that favor C sequestration. • For planted teak, the objectives of maximizing timber production and carbon sequestration are in conflict because the thinning schedules that maximize financial gains from C sequestration reduce economic gains from timber and vice versa. • With actual timber teak and market C prices, optimal NPVW is much larger than optimal NPVC. • For C prices under 40 $US MgC optimizing simultaneously for timber production and C sequestration is the best option, as additional although sub-optimal revenues can be obtained from C payments. • Lengthening the rotation, avoiding thinnings, or reducing their intensity increase carbon storage in planted teak, although, under the analyzed scenarios, after 120 yr. almost all carbon has been re-emitted to the atmosphere. Additional keywords: heuristics, genetic algorithms, operations research, forest management planning, stand level model, carbon stocks. Abbreviations used: C (Carbon); GA (genetic algorithm); NPVW, NPVC, NPVT (net present value from the cash flows of timber (wood), carbon, and total); SEV (Soil (land) expectation value); dbh (diameter at 1.3 m from the ground); G (stand basal area); Gp (potential site carrying capacity in terms of G); SQ (site quality); R (rotation, harvest age); A (age); I (thinning intensity); Vob, Vub (overbark, underbark volume); gr (basal area growth rate); r (interest rate); harvest and transport costs (Hc); Pc (C price). Authors ́ contributions: Conception, design, analysis, and interpretation of data (MAQM, MJR); design and programming of algorithms (MAQM); drafting manuscript, critical revision, statistical analysis (MJR, MAQM); provided the data (MJR). Both authors read and approved the final manuscript. Citation: Quintero-Méndez, M.A., Jerez-Rico, M. (2019). Optimizing thinnings for timber production and carbon sequestration in planted teak (Tectona grandis L.f.) stands. Forest Systems, Volume 28, Issue 3, e013. https://doi.org/10.5424/fs/2019283-14649 Received: 05 Feb 2019. Accepted: 20 Sep 2019. Copyright © 2019 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC-by 4.0) License. Competing interests: The authors have declared that no conflict of interests exist. Correspondence should be addressed to Mauricio Jerez-Rico: jerezorama@gmail.com Forest Systems 28 (3), e013, 14 pages (2019) eISSN: 2171-9845 https://doi.org/10.5424/fs/2019283-14649 Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA) Funding agencies/Institutions Project / Grant Consejo de Desarrollo Científico, Humanístico, Tecnológico y de las Artes (CDCHTA), Universidad de Los Andes, Mérida, Venezuela FO-748-17-01-B María-Alejandra Quintero-Méndez, Mauricio Jerez-Rico Forest Systems December 2019 • Volume 28 • Issue 3 • e013 2


Introduction
Given its high value as a precious tropical timber and its decline in natural forests, teak (Tectona grandis L.f.) is being planted at increasing rates in tropical regions of Asia, Africa, and America. So far, this species represents a small proportion of the total market of tropical timbers; however, it is the only precious wood tropical species that is becoming a commodity.
Commercial teak forest plantations must be managed intensively to obtain maximum profitability, however, an increasing awareness of society in environmental issues implies the need of reaching a trade-off between economic and environmental benefits. Today, many countries and private companies are investing in commercial teak plantations; however, there is a strong concern on global climate change, especially after the Paris 2015 Climate Change Conference, which raised the interest in planted forests as providers of environmental services, e.g., carbon sequestration. Planted forest for producing durable solid timber products can play an important role in carbon sequestration, as a large proportion of fixed C remains stored as solid, large size pieces of wood for long time after harvest. Thus, it should be appealing for planted teak forest managers to consider the additional potential benefits of producing timber and simultaneously providing environmental services such as C storage.
Within the frame of the intensive management of planted teak, managing the stand´s carrying capacity and optimizing the planting density, thinning, and pruning schedules is crucial for producing high quality timber in large logs with a high proportion of heartwood, as these have much higher value for the international markets.
Most knowledge on the effect of initial spacing and thinning schedules in teak growth and yield come from field trials based on permanent plots subjected to several combinations of spacing and thinning schedules, and from growers' experience and judgement. The latter is difficult to generalize; whereas, the former is limited by the number of testable combinations and circumscribed to specific sites. Mathematical models can integrate information from these experiences allowing generalization of growth responses and financial results to a large number of combinations including those for which no experimental studies exist. Moreover, models allow the assessment of the weight of intervening variables in the magnitude of observed changes in the biological, financial, and environmental outputs. Mathe matical models for analyzing responses to thinning schedules include stand density diagrams (e.g., Kumar et al., 1995;Jerez et al., 2003) and simulation models (e.g., Jayaraman & Rugmini, 2008;Tewari et al., 2014;Nӧlte et al., 2018). Less explored approaches for teak are optimization models, including classic optimization techniques (Mathematical Programming); meta-heuristics (e.g., genetic algorithms, simulated annealing), and multicriteria-decision-making techniques (e.g., goal programming) (Belavenutti et al., 2018, Pukkala & Kurttila, 2005. Thinning schedules can favor carbon sequestration in planted forests (Hoen & Solberg, 1994;Karjalainen, 1996;Pohjola & Valsta, 2007). Several works analyzed the carbon sequestration capacity of teak (Kraenzel et al., 2003;Gera et al., 2011;Takahashi et al., 2012;Sreejesh et al., 2013;Olayode et al., 2015;Nӧlte et al., 2018); however, for teak, optimization models incorporating a trade-off bet ween commercial (timber production) and environmental goals (carbon sequestration) considering biological and financial variables through meta-heuristics are rare. Quintero-Méndez & Jerez-Rico (2017) used a metaheuristics approach to develop a forest level (multiple stands) optimization model to determine alternative thinning schedules for a forest teak project. Optimize simultaneously for different objectives when managing planted forests can be a very complex task due to a large number of decision variables, constraints, and potential non-linear relationships among them. In this case, heuristic techniques can be applied as they can handle the model complexity more efficiently, although usually produce near optimal solutions rather than exact solutions.
Our objective was to develop and implement an optimization model integrating a growth and yield model with a heuristic optimization technique (genetic algorithms, GA) for determining thinning schedules that maximize the financial outputs (i.e., SEV, NPV) in teak plantations considering simultaneously merchantable timber production and C sequestration, and carbon sequestration only. In addition, the annual dynamics of C storage and emissions during and after the rotation is simulated for standing trees, solid products, and debris generated from thinning and harvest operations until the total re-emission of the stored carbon to the atmosphere. The model produces "optimal thinning" schedules under management scenarios varying in site quality (SQ), initial stand density at planting, variable rotation age (R), and financial variables: interest rates (r), harvest and transport costs (Hc), and timber and C prices. A comparative example shows technical and financial results for scenarios that combine rotation ages 20 and 25 yr., site index of 27 and 24 at base age 16 years, and initial spacings of 816, 1,111, and 1,600 trees ha -1 that are common for teak management in Venezuela, and other countries in Latin America and Africa.

Model description
The system represents planted, even aged teak stands managed with the primary purpose of producing solid merchantable timber and, in addition, storing carbon. Stands may differ in site quality, initial spacing, and rotation age (20 or 25 yr.). We developed a model that looks for combinations of thinning schedules (number, age, and intensity) that maximize the stand´s financial benefits. The model comprises three modules: 1) a stand growth and yield model (Jerez et al., 2015); 2) a module for computing carbon storage and emissions, and 3) an heuristic optimization module based on a Genetic Algorithm that looks for maximizing the financial benefits (NPV) for the goals of maximizing NPV W , NPV C , or both simultaneously (Fig. 1).

Growth and yield module
The growth and yield module consists of a wholestand model based on a system of differential equations describing the dynamics of the main state variables; e.g. top height (H, m), basal area (G, m 2 ha -1 ), and stand density (N, trees ha -1 ). The model projects growth in basal area, height, dbh, stand merchantable volume, total biomass, and stored carbon for various combinations of site quality, initial spacing, and thinning schedules. The basal area growth submodel follows the Pienaar & Turnbull (1973) approach to project the growth of thinned and unthinned plantations with a Richards´s growth equation (Jerez et al. 2015). The carrying capacity in terms of basal area (G p ) of the best sites for teak average 37.5 m 2 ha -1 (Zambrano et al., 1995, Bermejo et al., 2004. Dominant height growth was modeled as an anamorphic family of site index curves (base age = 16 yr.) that follows a Richards´s function fitted with the algebraic difference method (Clutter et al., 1983, Quintero et al., 2012. Site index classes were I = 24, II= 21, and III = 18 m. Stands with lower site index have corresponding lower G p . For SQs I, II, and III, G p is 37.5, 34, and 30 m 2 ha -1 , respectively. Quadratic diameter (dbh, cm) growth is computed from G and N. Stand average height is a function of dominant height and the ratio dbh of removed trees to dbh of remaining trees.
The mortality submodel has three components: a) density-independent mortality due to factors other than tree competition (e.g. drought, pests, and weeds); b) density-dependent mortality due to intraspecific competition; and c) removal of trees by thinning and final harvest. Density-independent mortality oc-

Carbon model
Calculates net C sequestration for management regimes from total C in biomass and C emissions Growth and yield model Computes stand growth and yield for a given management regime.
Heuristic optimizer Finds near optimal management regimes considering biological, financial, and technical variables.  curs only for trees between 0 and 3 yr.-old with decreasing probability. Density-dependent mortality is a decreasing exponential function of stand density. Estimates of equations´ coefficients came from permanent plots remeasured from 2 to 32 yr.-old established in the western plains of Venezuela. The model predicts the instantaneous change in stand quadratic diameter and the proportion of harvested trees with respect to removed basal area depending on a thinning selectivity coefficient th s , where th s = 1 for systematic thinning, th s < 1 for thinning from below, and th s > 1 for thinning from above (Jerez et al., 2015). Merchantable volume was calculated as overbark (V ob ) and underbark (V ub ) volumes in m 3 from stem base up to 5-cm diameter at tree top (Moret et al., 1998).

Carbon sequestration module
This module represents the yearly dynamics of stand carbon capture and emissions due to growth and death of above and belowground biomass. Carbon emissions come from the decomposition of branches, stumps, stems, and wood from dead or harvested trees throughout the stand rotation; plus decomposition of harvested forest products and wood losses from processing. Net carbon sequestration is the difference between carbon stored in biomass and carbon emitted to the atmosphere. Total above and belowground stored C was calculated from the total stand overbark volume using the conversion and expansion factors obtained by Kraenzel et al. (2003) for teak plantations. The model does not describe explicitly C stored or released from the soil, harvest operations, transport of products (i.e., carbon emitted by trucks and machinery), recycling, or effects of substituting fossil fuels by wood. Carbon emissions were estimated according to Hoen & Solberg (1994) who consider future C emissions as a function of the various decomposition and emission rates and organic lifetime (anthropogenic time) of the biomass components (categories). The carbon emission (E) in period i + j from category q removed in period i is: where RB i = removed biomass in period i, e k = the fraction of removed total biomass from category k, and Q k,i+j = emission rate of category k in i + j defined by: A Tk = anthropogenic time for category k, q k = annual decomposition fraction for category k: where D Tk = decomposition time of category k.
For each use category, emission and decomposition stop when a very tiny fraction (< 0.01 MgC ha -1 ) is emitted. Thus, a small fraction of biomass will remain without further decomposition, i.e., soil organic carbon (Hoen & Solberg, 1994). We assumed that carbon emission patterns are not appreciably affected by carbon mineralization processes.

Optimization module
This module determines the thinning schedule that maximizes NPV of cash flows related to stand timber production and carbon sequestration. Financial benefits of C sequestration come from the payment of environmental services made annually according to C prices following (Backéus et al., 2005, Díaz-Balteiro & Rodríguez, 2006, Baskent et al., 2008. The optimizer selects the best thinning schedule based on the outputs from the growth and yield model and from the C module ( Fig. 1).

Mathematical model
We developed a constrained optimization model that maximizes an objective function Z, where Z is the total net present value (NPV T ) of the cash flows through the rotation for a stand whose objectives are producing timber and sequestering carbon simultaneously: where NPV W = the net present value of the cash flows of timber (wood) from thinning and final harvest occurred through the stand life: [5] and NPV C = the net present value of the cash flows due to positive net C storage through the stand life: [6] where t = rotation age (yr.), r = interest rate, Cm i = establishment and maintenance costs in yr. considering the following decision variables: A j = number of yr. from planting (age = 0) to first thinning (j=1), or yr. between successive thinnings, where the subscript j is the j-th thinning, j = 2, 3, 4); I j = intensity of thinning j (j = 1, 2, 3, 4) as a percent of current ba sal area G; subject to the following constraints: [8] [9] [10] [11] [12] [13] Go i and Gf i are basal areas (m 2 ha -1 ) at the beginning and end of yr. i respectively subject to constraints (8-9-10); ΔG i+1,i = current annual increment in yr. i and Gthin i = removed basal area by thinning in yr. i. Constraint (11) specifies age = 3 yr. after planting as the minimum for carrying out the first thinning, and also the minimum interval between successive thinnings. Constraint (12) indicates that at least 25% G must be removed by a given thinning. The incomes by C sequestration (payment at market price per MgC) occur only the year in which sequestered C is larger than emitted C (Báckeus et al., 2005).

Optimization technique
We designed a heuristic procedure based on genetic algorithms (GA), a very robust optimization technique for solving efficiently many optimization problems (Dréo et al., 2006). The GA consists of searching throughout the possible solutions space by a process analogous to species evolution (Holland, 1975). The GA uses a binary codification to represent the possible problem solutions by generating a random initial solution and then applying a set of genetic operators: selection, crossing, and mutation.

Data
The data for growth and yield came from a network of permanent and temporal plots established on teak plantations in the western plains of Venezuela remeasured two or more times between 1.8 and 32 yr.old. Initial spacing varied from 2.0 × 2.0 to 4.0 × 4.0 m and thinning schedules comprised from 0 to 4 thinnings varying in intensity, age, and intervals of execution. Establishment, management, harvest, and operation costs are shown in Quintero-Méndez & Jerez-Rico (2017).

Model implementation
The model, implemented in Visual Basic 2015, comprises the growth & yield, carbon sequestration, and optimization modules. The program generates the best thinning schedules to optimize separately for timber production or carbon sequestration or for optimizing both objectives simultaneously. Inputs are site quality, initial spacing, rotation age, and desired number of thinnings, timber prices differentiated by diameter categories, carbon prices, establishment, main tenance and harvest operations costs, and interest rates. Outputs are optimal thinning schedules (age and intensity of thinnings), NPV T , NPV W , and NPV C . Soil expectation value (SEV) was calculated according to Bettinger et al. (2009): [14] where NPV = net present value, t = rotation age, and r = interest rate.
Optimal schedules are accompanied by the corresponding stand information on density, average tree diameter, basal area, dominant/average height, and merchantable volume on an annual basis. Outputs from carbon dynamics include storage, emissions, and annual C flows from stand initial conditions till 120 yr. after harvest considering standing trees, short, medium, and long term forest products, and wood debris from thinnings and final harvest.

Model runs
We determined thinning schedules for two optimization criteria: A) maximizing the NPV of timber production and carbon sequestration simultaneously (maximize NPV W + NPV C ); and B) maximizing the financial benefits associated with C sequestration only (maximize NPV C ). For each criteria 60 scenarios were defined combining SQ (I and II), initial planting density (816, 1,111, and 1,600 trees ha -1 ), thinning schedules (0 to 4 thinnings from below; i.e., the average dbh of removed trees was lower than that of the remaining trees for a given thinning (th s = 0.9), and harvest age (R = 20 and 25 yr.). Stand growth and yield was simulated by integrating the differential equations with the Runge-Kutta method, initial time was t 0 = 0 yr. at planting, initial G was calculated by multiplying the initial density times the root collar average diameter (1 cm). Initial stand height = 0.5 m.

Sensitivity analysis
We carried out a sensitivity analysis to determine the effects on the optimal solution (thinning schedule and maximum for the objective function) when assumed inputs (independent variables) were changed. Each input was changed within a given range and the other parameters kept fixed. The following parameters were varied: a) growth rate (g r ) at intervals of ±1%, (e.g., increases attributable to favorable climate conditions); b) thinning and harvest costs between ±10 and ±50% (base value =14.24 US$ m -3 ) at 10% intervals; c) r = 5, 8, 12, 14 %, base 10%); d) C prices (0,20,30,40,50,100,150, and 200 US$ Mg -1 C, base US$ 10); and e) ratios of timber prices per cubic meter among diameter classes. Teak timber prices depend largely on log size and age, as large logs with a high proportion of heartwood are preferred by the market (e.g. furniture, plywood). For young teak plantations, log dimension is a valid surrogate of wood quality. Four situations were considered: 1) all diameter classes have the same price in US$ m -3 (i.e., no premium for large diameter logs); 2) logs with diameter ≥ 25 cm are worth twice the price of logs 10 ≤ d ≤ 25 cm; 3) logs with diameter ≥ 25 cm are worth three times the value of 10 ≤ d ≤ 25 cm logs; and 4) logs with size d ≤ 10 cm have no value (Table 1).

Best thinning schedules
For the optimization criterion A (simultaneous maximization of NPV W and NPV C ), the largest SEV was reached for the scenario in SQ I, rotation (R) = 20 yr., and 1,111 trees ha -1 (US$ 14,542) and the NPV W is US$ 12,380, being the contribution of NPV C less than 3.5% ( Table 2). The SEV for scenario (SQ I, R= 20, 1,600 trees ha -1 ) was only slightly lower (US$ 14,318), but with a small increase in NPV C (4.2% of NPV T ). The scenario (SQ I, R = 20, 816 trees ha -1 ) had a considerably lower SEV (US$ 11,364). Thus, scenarios with 1,111 trees ha -1 always had the largest SEV as compared to the other stand densities. Also, stands in SQ I had always larger SEV than stands in SQ II; and scenarios with R = 20 had always larger SEV than R = 25. For optimization criteria B, i.e., only NPV C is optimized; the highest NPV C = US$ 763 is for scenario SQ I, 1600 trees ha -1 , and R =25; however, if we look for the largest SEV, then the best scenario is SQ I, 816 trees ha -1 , R = 20, despite that the NPV C is under the optimal. This is because the weight of the non-optimized NPV W , at base conditions, is very large when compared to the optimized NPV C (Table 2). Although longer rotations and higher densities increase the optimal NPV C , the SEV is strongly reduced, being as low as for scenario SQ II, R = 25, 1600 trees ha -1 . In this case, NPV C is even larger than NPV W .
The optimal number of thinnings in SQ II was always lower or equal than in SQ I scenarios given the same initial stand density and harvest age. On the other hand, the number of thinnings is always lower for R = 20 as compared with R = 25 yr. In the former case, 50% of scenarios showed a higher SEV when only two thinnings were executed. Conversely, for R =25, four thinnings  produce an optimal SEV 50 % of times. Furthermore, for the lower initial spacing (816 trees ha -1 ), in 50% of scenarios the best schedule is two thinnings. In contrast, with 1,111 trees ha -1 , in 50% of scenarios, the model indicated that the best schedule includes four thinnings. For 1,600 trees ha -1 , 75% of prescribed scenarios consis ted of three intensive thinnings (33, 31, 52% G) and R =20.
Overall, larger final diameters were reached for R = 25 years. The largest diameter was 42.1 cm for SQ I, 1,111 trees ha -1 , and four thinnings; whereas, the lowest diameter (30.3 cm) was reached for SQ =II, 816 trees ha -1 , harvest age = 20 yr. and two thinnings.
When only C storage was optimized, all scenarios showed no thinning schedules (Table 3). As occurred when maximizing simultaneously for C and timber, when optimizing only for C, scenarios with longer rotation age reached larger diameters, although these were comparatively low due to lack of thinnings. Thus, the highest dbh (27.6 cm) occurs for SQ I, 816 trees ha -1 , R = 25, and the lowest dbh (18.9 cm) occurs for SQ I, 1600 trees ha -1 , R = 20.
The curves for simulated G, stand density, dbh, and h for the scenario with the highest SEV (maximize NPV W + NPV C ) show a high contrast respect to the curves with the best SEV scenario that maximizes only C storage (Fig. 2).

Carbon sequestration and emissions
In addition to determining optimal SEV, NPV, and thinning schedules, the model generates information about annually stored and emitted C in standing trees and by type of product from initial planting to the time in which 99.99 % of all stored C has been released to the atmosphere for the scenarios that maximize NPV C and NPV W + NPV C (Fig. 3, A and B respectively).
Until harvest ages, most C remain stored in standing trees for both scenarios until harvest age. When maximizing only for NPV C , at age 25, just before harvest, C stored in standing trees peak at approximately 143 MgC ha -1 with only a small amount as debris from natural mortality (Fig. 3A). On the other hand for the scenario that maximizes NPV W + NPV C , the maximum stored C peaks around 13 years (94.1 MgC ha -1 ) with 73.4 MgC ha -1 in standing trees, and the rest in various products (Fig. 3B). By harvest age at year 20 only 103.5 MgC ha -1 remain stored, from which 55.6 is C in standing trees, and the rest in various products and debris. After peaking, in both cases C is emitted to the atmosphere till about 120 yr. when most C has been released (Fig. 3). In the period following thinnings or final harvest, the C stored in removed woody biomass is released according to the assumed emission rates. For this reason, when the objective is maximizing NPV W + NPV C (Fig. 3B), after the first and second thinning (ages 6 and 10) stored C in standing trees shows no or only a slight decrease in stored C because fast growing rates. After the third and fourth thinnings (ages 14 and 17); however, stored C showed relatively large reductions for standing trees due to the larger size of cut trees (Fig.  3B). For both optimization criteria, after the final cut, Table 3. Stand values for the best thinning schedules for A) Optimal combination of NPV W + NPV C and B) Optimal NPV C only (Interest rate =10%, timber prices (base) , carbon price = 10 US$ Mg -1 C). The best schedules are highlighted in grey. Site quality, R: Rotation age in yr., A: Thinning age (yr.), G: removed basal area (m 2 ha -1 ), I: thinning intensity (%G), dbh: quadratic diameter at breast height of harvested trees (cm), V ob : harvested over-bark volume (m 3 ha -1 ). C stored at a given moment begins to decrease due to C emissions from degradation of wastes remaining in the forest, wastes generated during the various stages of wood processing, and by the slower degradation rates of harvested products. By year 40 for both scenarios the amount of C is approximately equal to 40 MgCha -1 . This volume is comprised mainly by short and medium duration products (dbh = 21.3 cm) for scenario that maximizes C storage; whereas, for the other scenario, C remain stored in larger durability products (dbh = 35.0 cm). At age 60 from planting, most medium duration products have decomposed in both scenarios, thus, the scenario that maximizes C + wood, has a higher amount (20 MgCha -1 ) of C stored in long duration products. Afterwards, the remaining C is slowly released until 120 yr. when most of it has been reemitted in both scenarios, remaining only a small fraction that never decomposes or it is incorporated into the soil as organic C (Hoen & Solberg, 1994).

Sensitivity analysis
The sensitivity analysis for the model using as base the scenario with the maximum SEV (SQ I, 1,111 trees ha -1 , R= 20 yr., and four thinnings: ages 6, 10, 14, 17 and intensities 26, 28, 39, and 25% G) showed that SEV,   Figure 3. Carbon stored in standing trees, type of product and debris in a teak stand until 99.99 % has been released to the atmosphere. A) Carbon stored for the schedule that optimizes only the NPV C (SQ I, initial stand density = 1,600 trees ha -1 , R= 25 yr., no-thinning). B) Schedule that optimizes NPV W +NPV C (SQ I, initial stand density = 1,111 trees ha -1 , R= 20 yr., four thinnings ages 6, 10, 14,17, and thinning intensities 26 %, 28 %, 39 %, and 25 % of stand basal area).
NPV, and thinning schedules were mainly affected by changes in the interest rate and the growth rate (Table  4). Growth rate (g r ). The model was very sensible to variations in this parameter. In most scenarios, the optimal solution changes when g r varies ±1%.
Harvest costs (Hc). Changes within the chosen range for this variable did not affect the thinning schedules; however, they caused changes in the SEV of about 100US$ ha -1 for each change of ± 10% in harvest costs.
Interest rates (r). The optimal schedules were very sensitive to changes in r. For example, for r ≥ 12 %, the best schedule for SQ I, 1,111 trees ha -1 , R = 20 consisted of four thinnings. With r ≤ 8%, optimal solutions occur with only two thinnings, i.e., reduced interest rates favor lower number of thinnings. In addition the SEV increases sharply with lower r, e.g., from US$ 14,452 Table 4. Sensitivity analysis from the simultaneous optimization of NPV W + NPV C taking as basis the best scenario (SQ I, 1,111 trees ha -1 , R = 20 yr., and four thinnings: ages 6, 10, 14, 17 and intensities 26, 28, 39, and 25 % G (r =10%, C price = 10 US$ Mg -1 C, timber prices vary with diameter class). The largest SEV for each variable are highlighted in grey. when r = 10 % to US$ 41,755 for r = 5%. Conversely, increments in r, decrease the SEV, but in a lower magnitude (Table 4). Carbon prices (Cp). For C prices in the range 0 -40 US$ MgC -1 , the optimal thinning schedules did not change and the SEV remained almost constant. Therefore, the timber prices determine the best thinning schedules, as they have the largest weight on the objective function. Above 40 $MgC -1 the optimal solution changes by delaying the age and reducing the number or intensity of thinnings, and increasing the SEV.

Modified
Timber prices. Changing the relative prices among diameter categories changed the optimal solution by changing thinning number, ages, and intensities. When the difference in prices between the small and large diameter log was lower (Case 1 Table 1), the first thinning was done at earlier ages and greater intensity to produce monetary returns in the shortest time.

Optimal thinning schedules
As expected, plantations growing on SQ I, other conditions fixed, had the highest SEV; but also, more thinnings were needed that in lower site qualities. In this case, greater initial stand densities and longer ro tations also required greater thinning intensities. The sensitivity analysis showed that growth rate, interest rate, and timber price were the most influencing input variables to determine the best thinning schedules for each objective. Thus, focusing on a more precise estimation of these variables will increase the reliability of results. Also, for C prices 0-40 US$ MgC -1 , the best thinning schedule is the same, no matter how much C is stored, indicating that C sequestration is not essential for optimizing the thinning schedules, because its weight in the objective function is not significant as compared to the timber prices (53-400 US$ m -3 ). Thus, the optimal management schedules obtained for timber production and carbon sequestration (Table 2) can be compared with those reported in the literature for planted teak. Overall, our results agree with projected results of Pérez & Kanninen (2005) who propose three to four thinnings, each removing 25 to 50% of trees for 20-30 yr. rotations. Also, results agree with executing a first intensive thinning (40-60%) at ages 3-6 yr. (Chaves & Chinchilla, 1986, Jerez & Coutinho, 2017. The SEV values are in agreement with those reported for teak in other studies considering the range of interest rates (c.f., Restrepo & Orrego, 2015).
For teak plantations, the objectives of maximizing timber production and carbon sequestration are in conflict because the thinning schedules that maximize financial gains from C sequestration reduce economic gains from timber and vice-versa. Nepal et al. (2012) found a similar behavior when analyzing the financial relations between timber production and C sequestration in stands of Pinus taeda L. and Quercus pagoda Raf. Other authors (e.g., Báckeus et al., 2005, Keles & Baskent, 2007, Baskent et al., 2008, Raymer et al., 2009 found that when implementing management practices that increase C sequestration, the economic benefits from timber harvest decrease. In scenarios with economic incentives such as payments for environmental services, the best option is to choose schedules maximizing the joint NPV. In this case, NPV C will be lower than when maximizing only for C sequestration, but NPV W will be much higher, generating larger economic benefits, but keeping the environmental benefits of C fixation. Nӧlte et al. (2018) found a trade-off between C storage and economic return in planted teak in Costa Rica suggesting the need of creating economic incentives for increasing C sequestration that compensate for losses in timber production. For teak, with the current C assumed prices, carbon credits schemes such as those suggested by Báckeus et al. (2005) or Derwish et al. (2009 are insufficient because they represent only a small percentage (4-5% in our work) of the plantation total benefits.
The management schedule affects the C storage capacity of planted teak. Increasing rotation length, higher initial spacings, and fewer thinnings with lower intensity or no thinning, stored more C and generated larger benefits independently of site quality. These results agree with those from Nӧlte et al. (2018) for teak, and for other species (Liski et al., 2001;Pussinen et al., 2002;Kaipanen et al., 2004;Diaz-Balteiro & Rodríguez, 2006). For poor growing, non-profitable plantations, C storage could be attractive provided that it is paid as an environmental service. In this case, the best schedule would be no-thinning and delayed final harvest to reduce carbon emissions. Similar findings were reported by Lopera & Gutiérrez (2001) for Pinus patula, and Pohjola & Valsta (2007) in Pinus sylvestris.
In the analysis of C flows for a stand along the rotation under the model assumptions and scenarios, the average annual rate of C sequestration varied between 3.1 and 4.8 MgC ha -1 yr. -1 depending on the management schedule. These values agree with those reported for teak by the IPCC (1996), i.e., a mean annual increment of dry matter accumulation in planted forest equivalent to 8 MgC ha -1 yr. -1 representing a fixation rate of 4 MgC ha -1 yr. -1 . Also, they agree with those of Brown et al. (1986), who estimated a potential C fixation between 2.7 and 9.6 MgC ha -1 yr. -1 for tropical plantations.
For a 20-yr. rotation, the model estimated between 55.7 and 77.0 MgC ha -1 of C stored in standing trees depending on the thinning schedule and site quality. For a 25-yr. rotation this value fluctuated between 67.2 and 87.1 MgC ha -1 . Nӧlte et al. (2018) estimated between 76.9 MgC ha -1 and 89.5 MgC ha -1 for stands harvested at ages 20-25 respectively. Observed differences with our results are due mainly to the thinning schedule chosen by these authors (4, 8, 12, 18, and 24 yr. with remaining trees of 556, 333, 200, 150, and 120 trees ha −1 and initial stocking of 1111 trees ha −1 ) which differ from schedules generated by our model.
When the only objective was maximizing NPV C (SQ I, 1,600 trees ha -1 , R = 25 yr., and no-thinning), the mean annual carbon storage increased at a rate of 5.7 MgC ha -1 yr -1 , and the amount of C stored in standing trees at final cut was 143 MgC ha -1 . This result agrees with those of Kraenzel et al. (2003) for unthinned planted teak in Panamá (100-141 MgC ha -1 ). It is important to notice that for the scenario that optimized NPV W +NPV C a larger amount of C than for the scenario NPV C remained stored until 80 years because it was contained in long duration products.
The sensitivity analysis showed that optimal schedules and SEV are very sensible to interest rates (Table 4). Lower r (5-8 %) increased sharply the SEV and reduced the number of thinnings from four moderate thinnings to two more intensive thinnings. On the other hand, increases of r to 12 -14%, decreased SEV, but maintained the schedule of four thinnings, although the first one changed from age 6 to age 4. Increased growth rate increased SEV, although its effect in thinning schedules appears to be erratic, e.g., increasing g r by 2 % reduced the age of the first thinning to age 3; with an intensive thinning of 50 % intensity at age 16. On the other hand, reducing g r by 1% reduced the number of thinnings to three. Changing the relative timber prices also affected optimal schedules and SEV. Paying a high price for larger logs in relation to smaller logs changed from four moderate thinnings to two intensive thinnings. Harvest and transport costs did not affect the original thinning schedules. Finally, C prices only began to affect the schedules when the prices were above US$ 40.
According to the above arguments, our model generates reasonable and consistent results under the imposed assumptions and constraints, making it useful for analyzing the potential of planted teak forest for storing C in biomass during and after the harvest as long duration forest products, as well as analyzing the effect of thinning schedules and rotation age for timber and C sequestration.
Additional economic benefits to timber production can be obtained by accounting for C sequestration in teak plantations without substantial changes in the optimal schedules that maximize NPV W . In contrast, if the main objective was to maximize the benefits from C sequestration, moderate, infrequent, and late thinnings will be needed to reduce emissions. This leads to a much lower total SEV however, than that obtained by optimizing timber production. In scenarios with very low interest rates, long rotation (60-80 yr.) looking for very high quality veneer, and with an active market for carbon, this option could be attractive. If considering successive rotations on the same area for producing large duration solid timber products, teak plantations can become an efficient long-term system of C storage.