IntroductionTop

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 com­mercial 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 genera­lization 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 va­riables in the magnitude of observed changes in the biological, financial, and environmental outputs. Ma­the­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 tech­niques (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.

Material and MethodsTop

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 NPVW, NPVC, or both simultaneously (Fig. 1).

Figure 1. Model components and their relationships.

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, m2 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 (Gp) of the best sites for teak average 37.5 m2 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 Gp. For SQs I, II, and III, Gp is 37.5, 34, and 30 m2 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 ths, where ths = 1

for systematic thinning, ths < 1 for thinning from below, and ths > 1 for thinning from above (Jerez et al., 2015). Merchantable volume was calculated as overbark (Vob) and underbark (Vub) volumes in m3 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 RBi = removed biomass in period i, ek = the fraction of removed total biomass from category k, and Qk,i+j = emission rate of category k in i + j defined by:

ATk = anthropogenic time for category k, qk = annual decomposition fraction for category k:

where DTk = 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 emi­ssion patterns are not appreciably affected by carbon mineralization processes.

Seven carbon compartments were considered: 1) roots, 2) deadwood, 3) branches and stumps, 4) bark and debris, 5) short-term duration products (small poles for fencing), 6) midterm duration products (struts and large poles), and 7) long term duration products (sawn wood). Information on decomposition and anthropogenic times are from Hoen & Solberg (1994). Fractions of wood products per category came from Quintero-Méndez & Jerez-Rico (2017).

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 (NPVT) of the cash flows through the rotation for a stand whose objectives are producing timber and sequestering carbon simultaneously:

where NPVW = the net present value of the cash flows of timber (wood) from thinning and final harvest occurred through the stand life:

and NPVC = 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, Cmi = establishment and maintenance costs in yr. i, Bi = benefits in yr. i from timber harvested, Pc = carbon price (US$ Mg-1C), Fi = carbon fixed in yr. i (MgC), Ei = carbon emitted in yr. i (MgC), Td = time (yr.) in which 90% of C has decomposed, Ini = indicator variable (Ini = 1 if Fi > Ei, Ini = 0 otherwise). There are three maximization options for Z: maximizing NPVW only, maximizing NPVC only, or maximizing both values simultaneously:

considering the following decision variables: Aj = 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); Ij = intensity of thinning j (j = 1, 2, 3, 4) as a percent of current ba­sal area G; subject to the following constraints:

Goi and Gfi are basal areas (m2 ha-1) at the beginning and end of yr. i respectively subject to constraints (8-9-10); ΔGi+1,i = current annual increment in yr. i and Gthini = 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 possi­ble 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 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).

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), NPVT, NPVW, and NPVC. 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.

Model runs

We determined thinning schedules for two opti­mization criteria: A) maximizing the NPV of timber production and carbon sequestration simultaneously (maximize NPVW + NPVC); and B) maximizing the financial benefits associated with C sequestration only (maximize NPVC). 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 (ths = 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 t0 = 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-1C (equivalent to

2.72 US$ MgCO2 where 1 MgC = 3.67 MgCO2), 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).

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 (gr) 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-1C, 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).