Simulating Harvest Schedule for Timber Management and Multipurpose Management in Teak Plantations

Sustainable management of teak plantations in Java requires an improvement of the existing yield regulation method to optimize multiple benefits of the plantations at risk of stand destruction. This study was therefore aimed to formulate an alternative harvest scheduling model that integrates risk of stand destruction for supporting multipurpose management of teak plantations. The proposed model used a state-space planning model to simulate the dynamic of plantations due to timber harvesting and stand destruction, and then sought optimal solutions for 2 management scenarios, i.e. timber management that optimized total harvest volume and multipurpose management that optimized net present value (NPV) while increasing carbon stocks. Using a case study on a typical teak plantation, this study confirmed that increasing destruction rates reduced harvest volumes, NPV, carbon stocks, and resulted in imbalanced ending age-class structures. Reducing cutting-age limit increased harvest volumes and NPV, but it also reduced carbon stocks of the plantations. Although the multipurpose management generated lower financial benefit, it maintained carbon stocks and produced better ending age-class structures compared to timber management. The proposed harvest scheduling model provides a useful planning tool for managing teak plantations.


Introduction
In recent years, there has been an increasing public concern on sustainability of teak plantations in Java, Indonesia.This concern emerged mainly because most of the plantations are dominated by young stands (≤ 30 years old) and always face the risk of stand destruction (Tiryana et al. 2011a;Rohman 2014).In the presence of stand et al. destruction, it is unclear whether the plantations could achieve a sustained timber yield (Smartwood 2000) and perpetually produce social and environmental benefits.Thus, forest managers ( Perum Perhutani ) are facing a i.e.
, PP greater challenge to sustainably manage teak plantations.
Although the paradigm of sustainable forest manag m nt e e (SFM) has been recently adopted by P P , the erum erhutani sustained timber yield is still becoming a main concern on managing teak plantations.PP uses a neoclassical yield regulation method called Burns method (Perum Perhutani 1974), which is a combination of area and volume controls, to determine annual allowable cut (AAC) for a 10-years planning period.In this method (Perum Perhutani 1974;Parthama 1995), AAC area is simply calculated by dividing total standing forest areas with a rotation age, whereas AAC volume (AAC-v) is determined using steps from Burns 2 (1951).In the first step, an initial AAC-v is computed by dividing total harvest volumes at a mean-cutting-age with a rotation age.The mean-cutting-age is defined as the average weighted age (by age-class areas) of growing stocks plus one-half of the rotation age.In the second step, an iterative cutting-time test (also called tabular check) is carried out to determine the cutting-time of each age-class.The initial AAC-v is accepted if it results in a small deviation (≤ 2 years) between the sum of cutting-times of each age-class and the rotation age.Otherwise, the initial AAC-v is adjusted by repeating the cutting-time test several times, resulting in a tedious and cumbersome work.
The Burns method is relatively simple, but it has some drawbacks.First, it ignores the risk of stand destruction during a rotation period.Like other classical harvest scheduling methods, this method assumes no forest disturbances, so that a normal forest condition (balanced ageclass structure) might be expected after a certain rotation period.This assumption, however, is not valid for the current condition of teak plantations because stand destruction always occurs from time to time.The current age-class structures of teak plantations (similar to ) provide an Figure 1 empirical evident that the plantations have never reached a normal forest condition.Second, like other volume control methods, the Burns method only includes the objective of achieving harvest volumes (Davis 2001); hence, it has et al. no ability to include other management objectives, such as optimizing financial benefits or environmental benefits.In the context of SFM, forest managers have a greater challenge to be able to manage their forests for achieving sustainability of multiple benefits.Thus, the Burns method is less appropriate so that an alternative method is required to support a better management of teak plantations.
Few studies have attempted to develop alternative harvest scheduling models for teak plantations in Java.Parthama (1995) proposed a harvest scheduling model based on linear and chance-constrained programming methods to include risks due to timber yield variability of teak plantations.The proposed model, however, ignored the risk of stand destruction and only concerned with optimizing timber benefits.Kuncahyo (2006) proposed a simulation model using system dynamic approach by considering stand destruction and focusing on social benefits of teak plantations.He argued that the simulation model was more appropriate because the Burns method produced static and overestimated harvest volumes for disturbed plantations.Such model, however, requires many data on forest resources and socio-economic conditions, and formulating dynamic equations to describe relationships among model components would be not an easy task for forest managers.Although a simulation model provides a greater flexibility to develop management scenarios, such model may not produce optimal solutions because "there are no general solution algorithms to identify an optimal solution" (Buongiorno & Gilless 2003).Thus, an alternative harvest scheduling model is still needed to assist forest managers in achieving multiple benefits of teak plantations at risk of stand destruction.
The objective of this study was to formulate an alternative harvest scheduling model for simulating and optimizing multiple benefits (i.e.timber volume, financial benefit, and carbon sequestration benefit) of teak plantations at risk of stand destruction.Specifically, by using the proposed harvest scheduling model for a typical teak plantation, this study were aimed: 1 to analyze the effect of stand destruction and cutting-age limit on optimal harvest volumes, financial benefits, carbon stocks, and ending age-class structures, and 2 o evaluate trade-offs between timber management that t optimize harvest volume and multipurpose management that optimize financial benefit while maintaining carbon sequestration benefit of teak plantations.This study extended the basic model of Reed and Errico (1986) by including harvest volumes from thinning, financial benefits, environmental benefits from carbon sequestration, and various cutting-age limits, which were not considered in some previous harvest scheduling studies (Reed & Errico 1986;Armstrong 2004;Leduc 2014) et al. .

Methods
This study developed a harvest scheduling model using a case study in teak plantations of Kebonharjo Forest Management Unit (FMU).There are main phases in 3 developing the model: 1 imulating the dynamic of age-class structures from s period to period 2 ormulating optimization models to develop f management scenarios, and 3 sing the optimization models to simulate harvest u schedules for timber management and multipurpose management of the plantations.

Case study area
The Kebonharjo FMU is one of the most productive FMUs for producing high quality teak timber in Central Java.The total area of teak plantations in this FMU is 12,678.8ha with an imbalance age-class structure (Figure 1), indicating that young stands are more dominant (81%) than mature stands (Perum Perhutani 2007).Similar to other FMUs in Java, this FMU also uses the Burns method for determining AAC of teak plantations.For the period 2007-2016, the FMU uses a rotation of 60 years to obtain AAC volume of 11,168.6 m y r (with AAC area of 3 ea -1 69.1 ha y r ) from clear-cuttings and thinning volume of ea -1 1,921.2m y (Perum Perhutani 2007).

-1 ear
Model formulation Simulating dynamic of age-class : structures The dynamic of teak plantations from one period to another was simulated using a state-space planning model (Garcia 1990), which is also called Model III (Boychuk & Martell 1996;Davis 2001).In this planning model, the et al. state of plantations at period is described by plantation areas t in each age-class.This study used a 5 year age-class and period (with 15 age-classes) to minimize variation in stand yields due to area grouping and to account for 5-year thinning treatments.
During period , forest managers may conduct clear-t cutting of some productive stands, salvage-cutting of damaged stands, and thinning of standing forests.While harvested stands are then regenerated as age-class 1 (under an assumption of no delay in regeneration), remaining stands move to the next age-classes at period +1 (Garcia 1990) The area of age-class 2, 3, …, -1 is composed of the k surviving areas that were not harvested in period : 3 k k This study assumed that destruction rates at age-class ( ) j q j of a certain destruction level in matrix and were constant R S for all periods (Reed & Errico 1986)  In addition, a zero destruction rate was also used to analyze optimum harvest levels at no risk of stand destruction.

Developing optimization models
To find optimal solutions for harvest levels at risk of stand destruction, this study used a linear programming (LP) model, which was proposed by Reed and Errico (1986) .maximizing total harvest i.e volume and maximizing total net present value (NPV) of timber harvest while increasing carbon stocks.The objective of maximizing total harvest volume or NPV was explicitly formulated in the objective function of LP models, whereas the objective of increasing carbon stocks was formulated as a model constraint (Baskent 2008;McCarney et al. et al. 2008).

Total harvest volume
To optimize total harvest volume (THV) from clear-cutting, salvage-cutting, and thinning, the following objective function Eq was used in LP ( uation ) [10] models (modified from Reed & Errico 1986): , , and are volume-at-age vectors of clear-cutting, salvage-cutting, and thinning: The elements of these vectors are average timber volumes (m ha ) in age-class from clear-cutting ( ), salvage-

F NPV
NPV of salvage-cutting, NPV of thinning, and total fixed cost (for administration, production activities, social services, and environmental services) at period , t respectively, which were calculated using Eq uation [12a] Eq Eq Eq : uation , uation , and uation  , , and harvesting, planting, and fixed cost, respectively; and is the l length of period.
This study assumed an interest rate of 8% per year, timber prices, and activity costs that were constant over the planning horizon.Timber prices according to stand age were approximated using log prices based on log's diameter classes applied in the FMU.The log diameters at a certain stand age were estimated from the existing yield table at site class 3.5 (for logs of clear-cutting and thinning) and site class 1 (for logs of salvage-cutting).In addition, harvesting cost, planting cost, and total annual fixed cost were calculated based on financial data from the FMU.

Formulation of model constraints
The LP models used 4 constraints: area (Eq ), harvest volume flow, uation [ ] 9 minimum cutting-age limit, and carbon stocks flow.The use of these constraints depended on the management scenarios used in the LP models (Table 1).

Harvest volume flow
The sequential flow of clear-cutting volume (SFV) was used to control the variation of harvest volumes by allowing clear-cutting volumes ( ) to decrease V c or increase within a specified percentage from one period to the next (Reed & Errico 1986).The SFV constraint (Eq uation [ ] 13 ) reflected the objective of forest managers in maintaining continuous flows of timber volume from scheduled harvestings.

[ ] 13
This study used a 10% maximum proportional volume decrease ( ) and volume increase ( ) to facilitate a greater   flexibility for the LP models in finding optimal solutions.To accommodate the SFV constraint, the lower parts of matrix A and vector were modified as shown in Appendix 1c (Reed & b Errico 1986).

Minimum cutting-age limit
This study used minimum cutting-age limits (MCA, an analogue of rotation age) of 71, 61, and 51 years to analyze the effect of cutting-ages to harvest levels.To implement this constraint, the elements of vector at age-classes below a specified MCA limit (i.e.≤ j h c 14 for MCA 71 years, ≤ 12 for MCA 61 years, and ≤ 10 for j j MCA 51 years) were set equal to zero, meaning no clearcuttings at these age-classes.

Carbon stocks flow
The use of carbon stocks flow (CSF) constraints was aimed to reflect the objective of forest managers to maintain environmental benefits from carbon sequestration.The flow of total carbon stocks ( ) of standing C forests from one period to another was formulated using Eq : uation [ ] 14 [14] where and denote the maximum proportions of carbon γ θ stocks decrease and increase, respectively.This study used γ = 0 and = 0.1, meaning that the carbon stocks stored in θ standing forests were increase by 10% from period to period to represent the increasing role of forest ecosystems in maintaining environmental balances.The carbon stocks at period ( ) were calculated using Equation : where the elements of vector and were calculated using To accommodate the CSF constraint, the lower parts of matrix and vector were modified as shown in Appendix A b 1c.The average carbon stocks (tons ha ) at age-class that assumed four destruction rates that might occur during the planning horizon: zero (no destruction), low (4.4-6.3% per period, Equation 7a), medium (9.4-15.9% per period, Equation 7b), and high (18.5-49.6% per period, Equation 7c).Thus, to find optimal solutions for the possible combinations of scenarios, destruction rates, and cutting-age limits, this study analyzed 24 LP models as specified in Table 1.
Each LP model was run for a 140 year planning horizon (28 periods), which was 2-2.8 times the considered rotation lengths (50, 60, and 70 years); hence, it was a sufficiently long time to assess the sustainability of plantations (Clutter et al. 1983;Vanclay 2015).To avoid LP models liquidated forest stands at the last period (Garcia 1990), additional 25 years (5 periods) were added to the planning horizon but their solutions were then ignored from the models (Reed & Errico 1986;McCarney 2008).The LP models were et al. programmed in the R software (R Core Team 2015).

Results and Discussion
The LP-based harvest scheduling models generated optimal solutions for each scenario in relation to harvest volumes, financial benefits, carbon stocks, and ending ageclass structures of teak plantations.The optimal solutions of each scenario were not only affected by the rate of stand destruction but also by the cutting-age limit of teak plantations.
T he effect of stand destruction All scenarios confirmed that total harvest volumes (THV) during the planning periods generally decrease with increasing destruction rates.The decrease in THV corresponded to clear-cutting volumes (CCV) that decreased from period to period with increasing destruction rates (Figure 2  (Figure 3).This finding is reasonable, because when the plantations were damaged by various destruction agents (Tiryana 2011a), some portions of productive stands et al. were no longer available so that the areas and volumes of clear cuttings become considerably less.Meanwhile, the areas and volumes of salvage-cuttings increase with increasing destruction rates (Figure 2 and ) because Figure 3 the LP models emulated the PP's management policy (Perum Perhutani 1974) that permit the salvage-cutting of damaged stands at age-classes below an MCA limit.Reed and Errico (1986) also reported that annual harvest volumes and areas of clear-cuttings decreased when the probability of fire destruction increased.
Interestingly, except for high destruction, periodic harvest volumes (especially under TM scenario) tended to gradually increase from period to period until certain peaks at about rotation ages and then decrease towards the end of planning horizon.This behavior is reasonable because the plantations, which are initially dominated by young stands (Figure 1), continuously grow that eventually increase the growing stock of harvestable mature stands from period to period.Because of a lack of mature stands at initial periods, the LP models allocated less harvest volumes for some initial periods and then allocated more harvest volumes for the next periods when mature stands increased.Armstrong (2004) also confirmed that forests with less mature stands produced allowable harvest volumes that gradually increased from period to period in the presence of wildfire disturbances.This finding, however, contrasts to Reed and Errico (1986) who observed that annual harvest volumes at a certain destruction rate were initially decrease and then either slightly increase or relatively stable towards the end of planning horizon.A possible reason for this discrepancy is that the initial ageclass structure used by Reed and Errico (1986) had dominant mature and over-mature stands, which contradicts to the initial age-class structure of this study.For a forest with surplus mature stands, Armstrong ( 2004) also observed a gradual decline of harvest levels.These findings confirm that optimal harvest levels may also depend on the initial ageclass structures (Uusivuori & Kuuluvainen 2005).The uneven flow of timber volumes during a planning horizon could also reflect fluctuation in market demand that would not impact on the sustainability of plantations (Vanclay 2015).
The dynamic of harvest volumes due to stand destruction has affected total NPV, which decreased with increasing destruction rates (Table 2).Depending on the MCA limits, total NPV of TM scenario decreased up to 16.6-18.4%at low destruction rates, 26.3-33.4% at medium destruction rates, and 68.9-69.7% at high destruction rates.Similarly, the decrease in NPV of MM scenario reached up to 9.1-9.4% at low destruction rates, 38.5-45.8% at medium destruction rates, and 69.1-70.3% at high destruction rates.These findings confirmed that financial loss is inevitable when the plantations experienced stand destruction.Thus, to ensure a stable or increased income in the future, forest managers must be able to eliminate or minimize the risk of stand destruction.
Stand destruction greatly affected the carbon stocks of standing forests (Figure 4).At zero, low, and medium destructions, carbon stocks in each scenario generally increased up to period 10-15 and then decreased towards the end of planning horizon.At high destruction rates, however, all scenarios showed that standing forests stored similar amounts of carbon, which decreased from about 205,000 tons period (at the first period) to about 145,000 tons period -1 -1 (at the end of planning horizon).In each period, however, MM scenario generated more carbon stocks than TM scenario.The increase in destruction rates lead to the increase in salvage-cuttings, which then remove some amount of carbon stocks in standing forests.Similarly, thinning and clear-cuttings contribute to the removal of carbon stocks within a period.These results (Figure 4) confirmed that at high destruction rates, it seems impossible to store more carbon stocks in standing forests because the LP models allocated greater harvest areas from salvage-cutting that increased over time.These findings suggest that an effective strategy to save carbon in the plantations is to eliminate or reduce stand damage (Boscolo 2001;Martin et al. et al. 2015).
The effect of stand destruction is also obvious to age-class structures at the end of planning period (Figure 5, after regrouping the 15 age-classes for simplifying the results).The increase in destruction rates generally resulted in imbalanced ending age-class structures (EAS).In TM scenario, EAS were fairly good (where young and mature stands were exist in the plantations) when the destruction rates were zero (Figure 5a) or low (Figure 5b).But, the EAS of TM scenario at the other destruction rates showed negative exponential distributions (Figure 5c, Figure 5d), indicating that young stands were more dominant than mature stands.The MM scenario generally generated better EAS than those of TM scenario, in which mature stands (≥71 years old) were dominant when zero destruction rates (Figure 5e).At high destruction rates, however, all scenarios consistently showed similar EAS in which the plantations were mostly dominated (95-96%) by young stands (Figure 5d and Figure 5h).These findings are similar to those of Reed and Errico (1986), who observed that at no risk of stand destruction the EAS was close to that of a normal forest condition; but at the presence of destruction, the EAS was dominated by young stands.Parthama (1995) also observed imbalanced EAS for teak plantations in Cepu FMU.Thus, it seems to be impossible for teak plantations to reach a normal forest condition when stand destruction always occurs from time to time.Forest managers are then face a great challenge if they insist to achieve a normal age-class structure and a fully stocked forest with even flow timber volumes.However, forest managers may obtain fairly good EAS if they are able to eliminate stand destruction and willing to reduce harvest levels to some extents in order to increase standing forests (Boychuk & Martell 1996;Leduc 2014) and carbon et al. stocks from one period to another.
The effect of cutting-age limit Another factor affecting harvest volumes is a cutting-age limit, in which reducing the MCA limits generally increase total harvest volumes (Figure 2 and Figure 3) and total NPV (Table 2).When a lower MCA limit ( .51 years) is used in the LP models, more harvest e.g areas are available, greater harvest volumes are obtained, and higher NPV are then generated for each period.Parthama (1995) also reported that total NPV increased up to 43.9-96.8%when the rotation age of teak plantations was reduced from 80 to 60 years.Similarly, Boscolo .( 2001) et al showed that reducing diameter limits (analogue of MCA limits) of a natural forest increased NPV.
Although the reduction of MCA limit is financially profitable, this management option negatively affects the carbon stocks of standing forests (Figure 4), because the LP models allocated more harvest volumes for each period that resulted in greater removals of carbon stocks.Liski . et al (2001) also reported that reducing rotation lengths decreased carbon stocks of European forests.Therefore, increasing an MCA limit (lengthening a rotation period) would be an effective way to increase carbon stocks (Boscolo 2001;et al. Liski 2001;Raymer 2011), even though it reduces et al.
et al. the financial benefit of timber productions.
All scenarios also showed that reducing MCA greatly affected the pattern of EAS at zero and low destruction rates; but, it slightly affected the patterns of EAS at medium and high destruction rates (Figure 5).In TM scenario, reducing MCA from 71 years to 61 and 51 years removed the availability of mature stands over 70 and 60 years old, respectively, because the LP models allocated more harvest volumes and areas for mature stands in order to maximize the total harvest volumes.This is contradict to MM scenario, which still kept some mature stands over 80 years old (except at high destruction rates), because the optimization of timber volumes for generating NPV was constrained by the objective of increasing carbon stocks.As the result, timber harvests at each planning period only removed some portions of mature stands that accumulated to the end of planning horizon.
Trade-offs between timber management and multipurpose management Regardless of the effect of destruction rates and cutting-age limits, this study also reveals some trade-offs related to the management scenarios.In MM scenario, the financial benefit of timber productions  are generally lower than those of TM scenario; but, this management scenario provides greater environmental benefits from increasing carbon stocks of standing forests from period to period.Furthermore, the MM scenario produces better ending age-class structures (Figure 5e-h) compared with the TM scenario (Figure 5a-d).Many previous studies have also reported that some trade-offs occur when forests are managed for multiple-uses.For example, Baskent .( 2008) reported that a multipurpose et al management reduced the NPV of timber by 42% compared with a timber management.The lower financial benefit of MM scenario in this study is due to the lost of timber revenue that was not compensated by the financial benefit of increasing carbon stocks.In a recent climate initiative called r educing emission from deforestation and forest degradation (REDD+), forest managers could get financial benefits from conserving and enhanching carbon stocks.For production forests, however, forest managers need to determine appropriate opportunity costs and allocate most suitable areas for REDD+ projects to make a balance between timber production and conservation objectives (Aziz 2015).et al.Indeed, tradeoffs between timber and non-timber benefits are inherent features of SFM (Luckert & Williamson 2005).

Applicability of the proposed harvest scheduling model
The proposed harvest scheduling model provides some advantages compared to the Burns method.The model takes into account the risk of stand destruction, which is urgently required to determine a reliable AAC (Smartwood 2000) when the plantations always suffer from destruction.The proposed model also offers some flexibilities for forest managers to develop appropriate management scenarios by simulating several destruction rates, cutting-age limits, management objectives, and constraints.The flexibility of developing management scenarios is not an inherent feature of the Burns method, which is only capable to calculate annual harvest volumes.
This study only provides a case study of multipurpose management for optimizing financial and carbon sequestration benefits, partly because of a lack of basic models for quantifying other forest benefits.Further studies may extend the proposed harvest scheduling model by integrating other forest benefits.For example, the model can be extended to include social benefits of timber-sharing mechanism, which is currently practiced by PP through a community-based forest management program.Another possible extension of the model is integrating environmental concerns on maintaining protected areas, such as riparian zones, water springs, and wildlife habitats.

Conclusion
The proposed harvest scheduling model confirmed that optimal harvest volumes generally decreased with increasing destruction rates, resulting in financial losses.The increase in destruction rates resulted in imbalanced ending age-class structures.In the presence of stand destruction, it seemed impossible for teak plantations to reach a normal forest condition.Although reducing cutting-age limits increased the financial benefit of timber production, this management option resulted in greater removals of carbon stocks.The financial benefits of timber production were decrease under the multipurpose management.Nevertheless, the multipurpose management maintained carbon stocks and produced better ending age-class structures compared to the timber management.The proposed harvest scheduling model can be used to assist forest managers in determining optimal harvest volumes and developing relevant management scenarios for teak or other plantations at risk of stand destruction.Thus, it provides a useful forest planning tool to support SFM of plantation forests in Java or other regions with similar problems.
) For LP models that use CSF constraints coupled with area and SFV constraints, the elements of matrix and vector were

(
Figure 2).Similarly, in MM scenario CCV at MCA 61 years decreased from 41,330-113,210 m period at no destruction

Figure 2 Figure 3
Figure 2 Optimal harvest volumes of the timber management (TM) scenario from clear-cutting ( ), salvage-cutting ( ), and thinning ( ).The left (a-d), middle (e-h), and right (i-l) figures show the harvest volumes of four destruction rates (zero, low, medium, and high) at minimum cutting ages (MCA) of 71, 61, and 51 years, respectively.

Figure 4 Figure 5
Figure 4 Carbon stocks of standing forests at zero ( ), low ( ), medium ( ), and high ( ) destruction rates of the timber management (TM) scenario (a-c) and multipurpose management (MM) scenario (d-f).The left, middle, and right figures of each scenario show the carbon stocks at minimum cutting age (MCA) of 71, 61, and 51 years, respectively.
In the constraint of LP model, matrix (with × kN kN A elements) is composed of identity matrix (with × k k I elements), matrix and matrix while vector (with R, S, b kN × 1) is composed of matrix and vector which 1 R x , are concatenated as follows (see also Reed & Errico 1986): The upper parts (with k rows) of matrix and vector correspond to area constraints, the middle parts (with 2( -1) rows) correspond to SFV N N A b constraints, and the lower parts (with 2( -1) rows) correspond to CSF constraints.
greater than or equal to zero.This study extended such basic LP model by formulating objective 2 functions and imposing some constraints related to harvest volumes, minimum cutting-age limits, and carbon stocks.

Table 1
Specification of the LP models for each management scenario