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 whole-stand 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 ² 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 ² 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 SQ s I, II, and III, G _p is 37.5, 34, and 30 m ² 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­curs only for trees between 0 and 3 yr.-old with decreasing probability. Density-dependent mortality is a decreasing exponential function of stand den­sity. 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 ³ 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:

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:

and NPV ^C = the net present value of the cash flows due to positive net C storage through the stand life:

where t = rotation age (yr.), r = interest rate, Cm _i = establishment and maintenance costs in yr. i , B _i = benefits in yr. i from timber harvested, P _c = carbon price (US$ Mg ^-1 C), F _i = carbon fixed in yr. i (MgC), E _i = carbon emitted in yr. i (MgC), Td = time (yr.) in which 90% of C has decomposed, I n _i = indicator variable (I n _i = 1 if F _i > E _i , I n _i = 0 otherwise). There are three maximization options for Z : maximizing NPV ^W only, maximizing NPV ^C only, or maximizing both values simultaneously:

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 possi­ble problem solutions by generating a random initial solution and then applying a set of genetic operators: selection, crossing, and mutation.

The data for growth and yield came from a net­work 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) .

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) :

where NPV = net present value, t = rotation age, and r = interest rate.

Optimal schedules are accompanied by the corres­ponding 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.

We determined thinning schedules for two opti­mization 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 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.

The model was set to execute 50 runs for r = 10 %, harvest and transport costs ( Hc ) = 14.24 US$ m ^-3 and C price = 10 US$ Mg ^-1 C (equivalent to

2.72 US$ MgCO ₂ where 1 MgC = 3.67 MgCO ₂ ), commonly used in financial analysis including C sequestration ( Álvarez, 2009 ). Merchantable timber prices according to diameter category (Table 1) correspond to 2013 international market prices ( De Camino & Morales, 2013 ).

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).

For the optimization criterion A (simultaneous maxi­mization 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 thinning schedule with the best SEV , i.e., maximizing simultaneously NPV ^W and NPV ^C for SQ I, 1,111 trees ha ^-1 , R = 20, included four thinnings at ages 6, 10, 14, 17 with intensities 26, 28, 39, and 25% G (Table 3). However, for SQ II, the best schedule was for 1,111 trees ha ^-1 , R = 20, but only two thinnings (5 and 11 yr.-old) and removal intensities of 30.5 and

46.2% G ).

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 in­dicated 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).

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 stan­ding 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, 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 ).

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, 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 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 .

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.

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 econo­mic 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 ⁻¹ and initial stocking of 1111 trees ha ⁻¹ ) 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.