 Research Paper
 Published:
Optimal management of larch (Larix olgensis A. Henry) plantations in Northeast China when timber production and carbon stock are considered
Annals of Forest Science volume 75, Article number: 63 (2018)
Abstract
• Key message
The optimal management of larch (Larix olgensis) plantations in Northeast China consisted of 2 or 3 thinnings and a rotation length of 55–61 years when economic profitability, wood production, and carbon sequestration were simultaneously maximized. Wood production ranged from 5.4 to 11.7 m^{3} ha^{−1} a^{−1}, depending on site quality.
• Context
L. olgensis is an important tree species in the northeast forest region of China, playing a significant role in the establishment of fastgrowing and highyielding plantation forests in China. However, the management of these plantations has not been optimized in previous studies.
• Aims
The objective of the study was to find the optimal combinations of thinning times, thinning types, and rotation length for L. olgensis stands when both timber production and carbon stock are considered.
• Methods
First, a growth and yield model was developed to simulate the dynamics of larch plantations. Then, the models were linked with the Hooke and Jeeves optimization algorithm to optimize forest management for two commonly used planting densities and three site qualities.
• Results
Two thinnings were found to be suitable for larch plantations when the stand density at 10 years was 2125 trees/ha (corresponding to a planting density of 2500 trees/ha) whereas three thinnings were recommended when the density at 10 years was 2800 trees/ha (planting density of 3300 trees/ha). When the stand density was 2800 trees/ha, the optimal rotation length was 61, 58, and 55 years for site indices (SI) 12, 16, and 20 m (dominant height at 30 years), respectively. The mean annual wood production was 5.4 m^{3} ha^{−1} for SI 12, 8.2 m^{3} ha^{−1} for SI 16, and 11.7 m^{3} ha^{−1} for SI 20. The results were nearly the same for the lower initial stand density. The better the site quality of the stand, the earlier the thinnings were conducted.
• Conclusion
In multifunctional forestry, optimal rotation lengths of larch plantations were 10–20 years longer than advised in the current silvicultural recommendations for Northeast China.
1 Introduction
With the current fast economic growth, the demand for wood in China is increasing constantly, and the gap between supply and demand is becoming wider. One way to narrow this gap is the establishment of fastgrowing and highyielding plantation forests. As the main forest region in China, the northeast part of the country is responsible for both timber production and the protection of forest environments. According to the Eighth National Forest Resources Inventory (State Forestry Administration 2014), the area of larch plantations in Northeast China was 2.06 million ha, with a growing stock volume of 133.6 million m^{3}, which accounts for 36% of the total area and volume of plantation forests in the Heilongjiang, Liaoning, and Jilin provinces. However, many larch plantations have poor quality, low productivity, and low carbon stocks. Therefore, finding ways to improve the quality of plantations through forest management is an important scientific issue.
Larix olgensis A. Henry is a precious native tree species in Northeast China with strong adaptability, fast growth at young age, resistance against insects and diseases, and good technical properties of wood (Yu 2008). It is an important timber tree species in the northeast, northwest, north, and south subalpine areas of China (Chen 2010). It plays an important role in the establishment of fastgrowing and highyielding plantation forests in China (Zhu et al. 2013). Larch is widely used for housing, furniture, plywood, flooring, and decorative purposes (Chen 2010), and it is one of four coniferous pulp species in China.
Climate change has become an important environmental problem. Forests, as the largest terrestrial ecosystems on earth, not only provide wood for humankind, but also play an important role in sequestering CO2 from the atmosphere and mitigating the adverse effects of climate change. Since the multiple benefits of forests often correlate negatively (Jin et al. 2017), it is necessary to balance between the various functions of the forest to realize sustainable forest management. Previous studies on larch plantations have mostly concentrated on growth and yield modeling (e.g., Yasaka et al. 2008; Liu and Li 2010; Lei et al. 2016; Zang et al. 2016), and biomass and carbon storage (Ju 2010). There has been little research on management optimization (time, intensity and type of thinning, and rotation length) for multiobjective forestry.
The most representative study on plantation management in China is the one that was conducted by Li et al. (1997) who developed a simulation model for supporting the management of larch plantations assigned to the production of construction wood. Several studies on other species have been conducted elsewhere on multiobjective optimization of forest management, mainly focusing on the joint production of wood and nonwood benefits such as timber and bilberries (Miina et al. 2010), timber and scenic value (Pukkala and Miina 1997), timber and carbon sequestration (Olschewski and Benítez 2010), and cone and timber production (PasalodosTato et al. 2016).
To reach particular management objectives, plantation management is optimized by finding the ideal combination of thinning times, thinning types, and rotation length. Algorithms available for stand management optimization include dynamic programming, nonlinear programming such as the method of Hooke and Jeeves (1961), and populationbased methods (Pukkala 2009). Optimization of the management of a tree stand requires models for simulating stand development under alternative management schedules. Avery and Burkhart (1994) divide the growth models into standlevel models, diameter distribution models, and individual tree models. Individualtree models may be divided into distancedependent and distanceindependent models, depending on how the competition among trees is described.
The set of required models depends on the optimization approach so that dynamic programming is commonly used with standlevel growth models whereas nonlinear programming methods can use also other types of growth models. For example, Brodie et al. (1978) used dynamic programming with standlevel growth models, de Miguel et al. (2014) employed nonlinear programming with distanceindependent individualtree growth models, and Muchiri et al. (2002) used nonlinear programming with distancedependent individualtree growth models. Since individualtree models allow the most detailed simulation of stand development, they were used in this study, together with nonlinear programming.
The objective of this study was to develop optimizationbased management instructions for larch plantations in Northeast China, taking into account that management should be profitable, plantations should produce large amounts of wood and, at the same time, they should store significant amounts of carbon. As there are no suitable models for simulating the dynamics of larch plantations, the first step of our study consisted of fitting the necessary models required for the simulations. The maximized objective function included three management objectives that were maximized simultaneously, namely net present value, mean annual wood production, and average carbon stock of the plantation during the rotation.
2 Materials and methods
2.1 Simulation of stand development
The models needed to simulate the development of larch plantations included a site index model, individualtree diameter increment, survival and height models, a taper equation, and biomass models for different biomass components of a tree. Simulation of stand development began from a young initial stand. First, the site index model was used to calculate site index from the age and dominant height of the initial stand. Then, the individualtree height model was used to calculate the heights of all trees of the plot. The taper model was used to calculate the assortment volumes of each tree (saw log and pulpwood volume), and the biomass models were used to calculate the biomasses of trees, separately for stem, branches, foliage, and roots. The carbon content of each biomass component was used to convert biomass into carbon stock.
After calculating the growing stock characteristics for the initial stand, the development of the stand was simulated to the end of the rotation in 5year steps using the diameter growth model to increment diameters, and the survival model to calculate the proportion of surviving trees. The height, taper, and biomass models were used to calculate updated values for tree height, volume, and biomass. Thinning treatments were simulated during the rotation by removing a part of the trees.
2.2 Data for modeling stand development
Most L. olgensis plantations of China are located in the Heilongjiang and Jilin provinces of Northeast China, encompassing the mountains of Wanda (132° 48′ E–133° 56′ E, 46° 07′ N–47° 01′ N), Zhangguangcai (126° 31′ E–131° 44′ E, 43° 35′ N–46° 09′ N), and Changbai (127° 40′ E–128° 16′ E, 41° 35′ N–47° 57′ N). The climate is continental monsoon climate, with a mean annual temperature ranging from − 7 to 4.2 °C, mean annual rainfall from 475 to 800 mm, and elevation from 300 to 1800 m above sea level. The soils are mostly Haplumbrepts and Eutroboralfs (dark brown forest soil in the Chinese taxonomic system).
The data for modeling stand dynamics were collected from two sets of permanent sample plots established in the forest resources inventory of Heilongjiang province during 1990–2017. A set of 205 plots was used for site index and dominant height modeling and another set of 148 plots was used to model diameter increment, tree survival, and tree height. Each plot was located by GPS, and the topographic information (elevation, slope, aspect) was recorded. The plot size was always 20 m by 30 m (600 m^{2}).
The data for site index and dominant height modeling (Table 1) consisted of age and dominant height measurements in 205 permanent plots located in Changbai, Wanda, and Zhangguangcai mountains. Dominant height was calculated as the mean height of the 100 tallest per hectare, i.e., six tallest trees of the plot.
The data for diameter increment and survival models were obtained from the measurements of 148 plots, located in the larch plantations of Wanda and Zhangguangcai mountains (Fig. 1). The plots were measured at 5year intervals. Only one 5year interval per plots was used for modeling. The number of individualtree observations for modeling the 5year diameter increment was 7114, and 7838 observations were available for modeling the 5year survival (Table 2).
The data for modeling tree height consisted of 4288 tree height and diameter measurements in 59 plots (Table 2). These plots are a subset of the 148 plots that were used in diameter increment and survival modeling.
The data for the taper equation were measured in Zhangguangcai Mountain during 2008 and 2015 from 75 trees located in the proximity of 25 permanent plots differing in stand density, age, and site productivity. The diameter distribution of each plot was divided into five classes, each representing the same basal area. Then, one tree from each DBH (diameter at breast height) class was selected from the surroundings of the plot, cut down, and measured for taper data.
Diameter measurements from different relative heights (0 h, 0.02 h, 0.04 h, 0.06 h, 0.08 h, 0.10 h, 0.20 h, 0.25 h, 0.30 h, 0.40 h, 0.50 h, 0.60 h, 0.70 h, 0.75 h, 0.80 h, 0.90 h) resulted in 1200 diameter observations for developing the taper model (Table 3).
Of the 75 trees felled for the taper modeling, 68 trees were also measured for biomass (Table 3). The stems were cut into 1m sections and each section was weighed in the field. A 2–3cmthick disc was cut from the end of each section, weighed, and taken to the laboratory for determining the dry matter content. All live branches were weighed. Then, 3–5 branches were selected from the lower, middle, and upper part of the crown and weighed separately for branch and foliage mass. Samples of 50–100 g were taken from each selected branch, weighed, and taken to the laboratory for the determination of dry matter content.
Tree roots were excavated using both a lifting machine and manual digging. The zone of excavating roots was limited to a circle of 3 m in radius, and fine roots (< 5 mm) were not included. The roots of the sampled trees were divided into three categories: large roots (diameter ≥ 5 cm), medium roots (diameter 2–5 cm), and small roots (diameter < 2 cm). Samples of 100–200 g were taken from each category, weighed, and taken to a laboratory for the determination of dry matter content. All samples were ovendried at 80 °C and weighed. The dry biomass of a biomass component was calculated by multiplying the fresh weight of the component by its dry mass/fresh mass ratio (dry matter content).
The carbon content of each sample was measured with a Multi C/N 2100S carbon nitrogen analyzer. Using these measurements, the average carbon content was calculated separately for stem, root, branch, and foliage samples. These carbon contents were used in simulations, to convert the dry biomass of growing stock into carbon stock.
2.3 Modeling stand dynamics
The site index model was constructed using the guide curve method (Clutter et al. 1983). The guide curve shows the average development of dominant height along stand age. Several equations proposed in literature were tested, and the following variant of the Mitscherlich function (Mitscherlich 1919; Elfving and Kiviste 1997) was used to model the relationship between stand age (T) and the mean dominant height of the plots (HD_{guide}):
Diameter increment was modeled in the same way as described in Jin et al. (2017), for example. Diameter increment was modeled as a function of tree size, site productivity, and competition (Vanclay 1994) using multiple linear regression analysis. Tree size was described by DBH (d), site productivity by site index (SI), and competition by stand basal area (G) and basal area in larger trees (BAL). BAL is the basal area of those trees that are larger in DBH than the tree for which growth is predicted (Vanclay 1994). Logarithmic transformation was applied to the response variable to guarantee that growth predictions are never negative. To describe diameter increment, the logarithm of the difference of squared diameters at the end (d_{5}) and beginning (d_{0}) of a 5year period was used:
Binary logistic regression analysis was used to model the probability that a tree will survive for the coming 5year period (s). This model form guarantees that the predicted survival probability is always between 0 and 1. The predictors were the same as used in diameter increment modeling:
The Näslund function (Näslund 1937) was selected as the heightdiameter curve. The equation has shown a good fit for different species (Mehtätalo et al. 2015) and has been already used in Heilongjiang for Korean pine (Jin et al. 2017). The model is as follows:
where h is tree height (m), d is DBH (cm), and a and b are the model parameters. To describe the change of the height curve during stand development, parameters a and b were made dependent on the dominant height of the stand (HD):
The Kozak model was used as the taper model as it has performed well in previous studies (Kozak 2004; Heiðarsson and Pukkala 2011; de Miguel et al. 2012):
where d_{i} is the diameter (cm) at any given height h_{i} (m), q is h_{i}/h, t = 1.3/h, and b_{0}–b_{8} are parameters.
The modeling methods for biomass were similar to those in Dong et al. (2016). The following system of additive equations was simultaneously fitted for different biomass components:
where w is the ovendry biomass of the biomass component of a tree; subscripts r, s, b, f, and t represent root, stem, branch, and foliage, respectively; d is diameter at breast height and β_{i}_{0} and β_{i}_{1} are parameters estimated in regression analysis; and ε_{i} is the model error.
2.4 Initial stands
Simulations were started from young (10 years old) simulated initial stands. Three different site indices (dominant heights 12, 16, and 20 m at 30 years) were analyzed in this study. In each site index, the planting density was 3300 or 2500 ha^{−1}, which are the most common planting densities in China (Table 4). Assuming an 85% survival rate after 10 years since planting, the initial stand densities were 2800 and 2125 trees/ha.
2.5 Economic data
The economic data required for calculating the net present value (NPV) were the discount rate, timber prices, and the costs of site preparation, planting, and tending. Timber prices and treatment costs were based on interviewing the managers at the Mengjiagang Forest District in Heilongjiang Province. The site preparation cost was 2000 CNY/ha in year 0, planting cost was 5000 CNY/ha in year 0, and the tending cost was 1500 CNY/ha in year 10. CNY stands for Chinese Yuan Renminbi (sometimes abbreviated as RMB). In January 2018, one CNY was worth 0.13 EUR. The specifications and prices of timber assortments are shown in Table 5.
2.6 Objective function
Both singleobjective and multiobjective optimizations were used in this study. The following utility functions were used in the multiobjective optimization:
where NPV, WP, and C are net present value, mean annual wood production (total volume harvested during rotation divided by rotation length), and average carbon stock during a rotation, respectively. NPV_{max} is the NPV when it is the only objective (CNY ha^{−1}). Similarly, WP_{max} (m^{3} ha^{−1} a^{−1}) and C_{max} (t ha^{−1}) were obtained by maximizing WP and C as the only objective.
Therefore, all variables were divided by their maximum value to eliminate the influence of units. MOF_{2} was the baseline objective function based on the opinions of the authors and local experts on the importance of the three management objectives. The other functions (MOF_{1} and MOF_{3}) were used to analyze the sensitivity of management to the importance of management objectives.
2.7 Optimization
Optimizing stand management is equal to finding the optimal values for the decision variables.
Decision variables (DV) refer to a set of controllable variables that define the stand management schedule. The decision variables of this study were:

For each thinning (three decision variables):

Years since previous thinning, or if it is the first thinning, stand age at thinning

Parameter a_{1} of the thinning intensity curve (see below)

Parameter a_{2} of the thinning intensity curve

For final felling:

Years since the last thinning
Thinning intensity and thinning type were defined with the following logistic function, which gives the proportion of trees removed from different diameter classes (Pukkala et al. 2014b; Pukkala 2015; Jin et al. 2017):
where TI(d) is the proportion of harvested trees when DBH is d cm; and a_{1} and a_{2} are parameters to be optimized. Parameter a_{2} gives the diameter at which thinning intensity is 0.5, and a_{1} defines the type of thinning. When a_{1} is negative, the thinning is low thinning (thinning from below), a_{1} > 0 corresponds to high thinning (thinning from above), and a_{1} = 0 leads to uniform thinning (the same thinning intensity in all diameter classes). The number of decision variables (NDV) depends on the number of thinnings (N_{Thin}) as follows: NDV = 3 × N_{Thin} +1.
The optimization algorithm used in this study was the direct search method of Hooke and Jeeves (Hooke and Jeeves 1961), which has been widely used in previous optimization studies (e.g., Roise 1984; Valsta 1990; Pukkala and Miina 1997; Hyytiäinen 2003; Palahí et al. 2009; PasalodosTato 2010; Jin et al. 2017). A detailed description of the Hooke and Jeeves method is given in Appendix.
Data availability
The datasets generated and/or analyzed during the current study are not publicly available at the time of publication due to ongoing research project. However, partial data are included in the supplementary information files.
3 Results
3.1 Models for stand development
The model fitted for average dominant height development was as follows (Fig. 2):
The root mean square error (RMSE) of the model was 2.07 m. The base age was selected as 30 years. The formula for the site index of a sample plot is therefore:
where HD(T) is the measured dominant height at T years and HD_{guide}(T) is the prediction of the guide curve model at T years. The model was converted to the following form to simulate the development of dominant height (Jin et al. 2017):
After fitting, we obtained the following model for diameter increment:
where I_{d} is the 5year diameter increment (cm). The RMSE of the model was 0.41 and its R^{2} was 0.49. The model showed that diameter increment first increased, then decreased with increasing DBH (Fig. 3). DBH grew fastest when the diameter was 5–15 cm. As one would expect, there was a negative correlation between diameter growth and competition factors (G and BAL). The rate of diameter increment correlated positively with site index.
The model for the probability that a tree survives for 5 years (s) was as follows:
The AUC (area under receiver operating characteristic curve, see Heagerty and Zheng 2005) of the model was 0.847, which indicates good performance. The probability of correct predictions was 93.8% in the modeling data. As seen in Fig. 4, the survival probability increased with increasing DBH and site index, but stronger competition decreased the probability of survival.
The following heightdiameter model was fitted to calculate the individual tree height:
The RMSE of the model was 1.49 m. It can be seen in Fig. 5 that tree height increased with the DBH under the same dominant height, but only slowly when the DBH was greater than 30 cm.
The DBHheight curve “climbed up” when the stand developed and its dominant height increased.
The Kozak model (Kozak 2004) fitted for stem taper was as follows:
where d_{i} is the diameter (cm) at height h_{i} (m), d is the DBH of the tree (cm), h is the total tree height (m), q is h_{i}/h, and t = 1.3/h (RMSE = 0.96 cm).
The parameters of the biomass model (Eq. 7) for different biomass components of the tree are shown in Table 6, which also shows the carbon content of different parts of a tree.
The biomass models give the dry weights of the biomass components in kilograms. The unit of the predictor (DBH of the tree) is centimeter.
3.2 Optimal number of thinnings
Optimizations with the multiobjective utility function (MOF_{2}) with 0, 1, 2, 3, and 4 thinnings under different discount rates (1%, 2%, 3%) for different initial stand densities (2800 and 2125 trees/ha) showed that utility increased with the number of thinnings, but the increase was slow after 1–2 thinnings. Taking site index 16 m as an example, the effects of using 0, 1, 2, 3, or 4 thinnings on the objective function value in the two initial stands are shown in Fig. 6. Based on these results, we suggest three thinnings for an initial stand density of 2800 trees/ha, and two thinnings for 2125 trees/ha. These numbers of thinnings were used in all subsequent analyses in this study.
3.3 Effect of management objective on optimal management
The results of the singleobjective (NPV or wood production) and multiobjective optimization for different site indices and initial stand densities are compared in Tables 7 (for 2800 trees/ha) and 8 (for 2125 trees/ha). In all cases, NPV was calculated with a 2% discount rate. The results showed that the timber production of the maximum WP schedule was only slightly higher than in the schedule for maximal NPV. Rotation length was extended by 1–5 years, and the average carbon stock increased by 3–10 t ha^{−1}when WP was maximized instead of NPV. Maximizing the baseline multiobjective utility function (MOF_{2}) further increased rotation length by 4–13 years compared to the max NPV schedule. The average carbon stock also increased, but wood production was not much affected. Wood production, net income, NPV, and carbon stock all increased with increasing site index under the same objective function. Increasing the weight of carbon stock in the multiobjective utility function extended the rotation length.
The tradeoffs between the different outputs of larch plantations are visually depicted in Fig. 7 for site index 16 m. It can be seen that wood production was not sensitive to the objective function used in optimization. The mean carbon stock of the forest and economic profitability (NPV) were most in tradeoff. Maximizing NPV led to a 30–35% decrease in the average carbon stock of the forest, when compared to maximizing MOF_{3}, where the weight of carbon stock was 0.33. The baseline objective function, MOF_{2}, was good in all other outputs except carbon stock, where it was 5–10% worse than in the management schedules that maximized MOF_{3}.
When MOF_{2} was maximized in site index 16 m, the three thinnings took place at stand ages of 27, 37, and 48 years, respectively, when the initial stand density was 2800 trees/ha (Fig. 8). All thinnings were from above, and the DBH for which thinning intensity was 50% was 15.3 cm (first thinning), 17.9 cm (second thinning), and 19.8 cm (third thinning). The rotation length of the multiobjective schedule was a few years longer than in the schedules that maximized NPV.
The management schedules that were evaluated as the most recommendable are shown in Table 9 and Fig. 9. In general, thinnings began earlier in sites of better quality, which is logical since trees grow faster with good site quality. For the same site, intense competition commenced earlier with larger initial density. Therefore, thinnings should be conducted earlier in the stands with higher initial density.
3.4 Comparison with current silvicultural recommendations
The maximum thinning intensity is usually not more than 45% in Chinese tree plantations, and the thinnings are conducted as thinning from below. To evaluate the optimality of the current silvicultural recommendations, optimizations for initial stand density of 2800 trees/ha with 2% discount rate and three thinnings were repeated using constraints for maximal thinning intensity (max 45%) and thinning type (only thinning from below was accepted) for all site indices using MOF_{2} as the objective function. The constraints reduced economic profitability (by 1.3–11.8%), timber production (0.4–5.5%), and utility (1.4–15.4%), but increased the average carbon stock of the rotation by 1.8–15.5%. The two constraints shortened the rotation length by 1–10 years. This is logical given that thinning from below increases the average financial maturity of trees (Fig. 10).
4 Discussion
This study developed management instructions for L. olgensis plantations in Northeast China for multifunctional forestry. Three management objectives were simultaneously maximized: net present value (NPV), wood production (WP), and carbon stock. Whereas NPV describes the profitability of forestry from the viewpoint of forest landowners or managing organizations, WP indicates how much wood is supplied by forestry for the various needs of society. Forest landowners are not necessarily worried about wood production, but NPV is a very important criterion for the landowner as it describes the economic efficiency of forest management.
Maximizing NPV as the only objective may lead to low volume production. For example, if stand establishment costs are high and the discount rate is also high, maximal NPV is reached with minimal investments in stand establishment and young stand management, which may lead to low wood production. Since NPV measures the benefits to the landowner and WP is related to societal benefits, and the two variables do not necessarily correlate strongly, it was justified to include both NPV and WP in the objective function.
Carbon storing is another benefit that forestry generates to society. Using the average carbon stock of the rotation as the third objective variable leads to favoring management schedules where the growing stock volume is high. As seen in Tables 7, 8, and 9, the carbon goal also leads to longer rotations and higher net incomes, most probably due to a greater share of saw logs and largesized timber assortments. Largesized timber and a high share of saw logs make it possible to use much of the harvested wood for construction, and the constructed buildings serve as longterm carbon stores. A gradual shift towards longer rotations, higher growing densities, and saw logoriented management would increase the carbon stocks for many decades, both in forests and in different woodbased structures. Therefore, the management instructions developed in this study serve as a good compromise between economic efficiency, wood supply, and carbon sequestration.
Our results are in agreement with many earlier studies, which show that thinning from above is more profitable than thinning from below (e.g., Haight et al. 1985; Roise 1986; Haight and Monserud 1990; Solberg and Haight 1991; Valsta 1992; Vettenranta and Miina 1999; Hyytiäinen et al. 2005; Miina et al. 2010; Pukkala et al. 2014a, b; Pukkala 2015, 2017; Vauhkonen and Pukkala 2016; Jin et al. 2017). The effects of discount rate and site productivity on optimal management were also in line with previous studies: optimal rotation lengths are shorter for good sites and high discount rates (e.g., Jin et al. 2017; de Miguel et al. 2014; Palahí and Pukkala 2003; PasalodosTato and Pukkala 2007; PasalodosTato et al. 2009). Previous studies have shown that increasing importance of carbon sequestration leads to longer optimal rotations, higher growing stock volume, increased harvesting of saw logs, and decreased harvesting of pulpwood (Liski et al. 2001; Pohjola and Valsta 2007; Pingoud et al. 2010; Pukkala 2011; Pukkala et al. 2011; Niinimäki et al. 2013; Pihlainen et al. 2014). These earlier findings correspond to our results on the effect on carbon sequestration on optimal stand management.
We are aware that there were some limitations in our data, which affected the simulation of the dynamics of larch stands. Since the modeling data did not have mortality among large trees, the 5year survival rate was still greater than 0.95 when tree diameter was larger than 25 cm, and survival rate never started to decline when the trees became larger and older. This made the optimal rotation for maximal carbon stock very long. As this result was partly a consequence of inadequate models, the maximum weight of carbon stock was restricted to be 0.33. With this constraint, all management schedules reported in this study were reasonably well within the ranges of variation in the modeling datasets. Despite this, all the models developed and used in this study should be regarded as preliminary. New models should be developed when older plantations become available, and the number of permanent sample plot increases.
Even though the models used in optimization may have limitations and were based on small datasets, it was important to fit the models and conduct the optimizations reported in this study. This is because sciencebased instructions for multifunctional management of larch plantations are urgently needed right now. Therefore, the results presented in this study are important for forestry practice in Northeast China, despite some potential limitations in the models.
The present forest harvesting and regeneration regulations of China recommend final felling for larch plantations at 41–50 years. However, these clearcutting ages should not be interpreted as fixed rules, but as flexible management guidelines, which may be modified depending on planting density, number, and type of thinnings, site quality, and several other factors. The clearcutting age should also be adjusted in accordance with timber prices and the discount rate. Li et al. (1994) pointed out that larch plantations should rarely aim at the production of large logs, but mainly mediumsized and small logs. From the technical and economic maturity viewpoints, rotation lengths as short as 17–21 years have been proposed for L. olgensis plantations in the production of small logs while 31–35 years are recommended for the production of mediumsized logs. Chen (2010) concluded that site index 16 m is suitable for the production of mediumsized logs, and three thinnings are optimal for planting densities of 3300–4400 trees/ha. Stands of site index 20 m may be used to produce largediameter timber, and the stand should be thinned four or five times during a rotation of 45–54 years.
In our study, the optimal rotation length for maximal wood production was 50–56 years. The use of the multiobjective function increased the rotation length by 1–16 years when compared to the management schedule that maximized timber production. Rotation length increased with increasing weight of carbon stock, which is a logical result.
5 Conclusion
This study presented a set of growth and yield models for larch plantations of Northeast China and used the models in management optimization. The study determined the optimal combination of thinning times, thinning types, and rotation length when both timber production and carbon stock are considered. Two thinnings and three thinnings were recommended for planting densities of 2500 and 3300 trees/ha, respectively. Wood production, net income, carbon stock, and NPV increased with the increase of site quality under the same objective function. In multiobjective optimization, increasing the weight of carbon stock led to longer rotation length. Wood production was not sensitive to the objective function while NPV and carbon stock were the most sensitive. Competition between trees starts earlier with higher density and better site quality of the initial stand. Therefore, the thinnings are to be conducted earlier in stands with better site and higher initial density. The two constraints of silvicultural recommendations, namely a maximum thinning intensity of 45% and forced low thinning, reduce economic profitability, timber production and utility, but increase the average carbon stock of larch plantations.
References
Avery TE, Burkhart HE (1994) Forest measurements, 4th edn. McGrawHill, Inc, New York
Bazaraa MS, Sherali HD, Shetty CM (1993) Nonlinear programming: theory and algorithms, 2nd edn. Wiley, Hoboken 639 p. ISBN 0471557935
Brodie JD, Adams DM, Kao C (1978) Analysis of economic impacts on thinning and rotation for Douglasfir, using dynamic programming. For Sci 24:513–522
Chen DS (2010) The optimal management mode of large and medium diameter timber in larch plantation. Ph.D thesis, Northeast Forestry University. in Chinese
Clutter JL, Forston JC, Piennar LV, Brister GH, Bailey RL (1983) Timber management—a quantitative approach. Wiley, New York, p 333
de Miguel S, Mehtätalo L, Shater Z, Kraid B, Pukkala T (2012) Evaluating marginal and conditional predictions of taper models in the absence of calibration data? Can J For Res 42:1383–1394
deMiguel S, Pukkala T, Yesil A (2014) Integrating pine honeydew honey into forest management optimization. Eur J For Res 133(423):432
Dong LH, Zhang LJ, Li FR (2016) Developing two additive biomass equations for three coniferous plantation species in Northeast China. Forests 7:136
Elfving B, Kiviste A (1997) Construction of site index equations for Pinus sulvestris L. using permanent plots data in Sweden. For Ecol Manag 98:125–134
Haight RG, Monserud RA (1990) Optimizing anyaged management of mixedspecies stands. II. Effects of decision criteria. For Sci 36:125–144
Haight RG, Brodie JD, Dahms WG (1985) A dynamic programming algorithm for optimization of lodgepole pine management. For Sci 31:321–330
Heagerty PJ, Zheng Y (2005) Survival model predictive accuracy and ROC curves. Biometrics 61:92–105
Heiðarsson L, Pukkala T (2011) Taper functions for lodgepole pine (Pinus contorta) and Siberian larch (Larix sibirica) in Iceland. Icel Agric Sci 24:3–11
Hooke R, Jeeves TA (1961) “Direct search” solution of numerical and statistical problems. J Assoc Comput Mach 8:212–229
Hyytiäinen K (2003) Integrating economics and ecology in standlevel timber production. Finnish Forest Research Institute, Research Papers 908: 42p.+ pp
Hyytiäinen K, Tahvonen O, Valsta L (2005) Optimum juvenile density, harvesting and stand structure in evenaged Scots pine stands. For Sci 51:120–133
Jin XJ, Pukkala T, Li FR, Dong LH (2017) Optimal management of Korean pine plantations in multifunctional forestry. J For Res 28:1027–1037
Ju WZ (2010) Age effects on stand biomass and carbon storage of Larix olgensis plantation: a case study in DongzheLengHe forest station of Yichun City. Master thesis, Beijing Forestry University. in Chinese
Kozak A (2004) My last words on taper equations. For Chron 80:507–515
Lei XD, Yu L, Hong LX (2016) Climatesensitive integrated stand growth model (CSISGM) of changbai larch (Larix olgensis ) plantations. For Ecol Manag 376:265–275
Li M, Li CS, Li CC, Zhang YC (1994) Study on the rotation length of larch plantation. Sci For logging 4:1–6 (in Chinese)
Li M, Li Y, Man DB, Deng BZ, Pan JZ, Yao JQ, Zhu YJ, Zou DD (1997) A computer system of the stand silviculture models of building wood for Larix olgensis plantation. J Northeast Forestry Univ 25:26–29 (in Chinese)
Liski J, Pussinen A, Pingoud K, Mäkipää R, Karjalainen T (2001) Which rotation length is favourable to carbon sequestration? Can J For Res 31:2004–2013
Liu W, Li FR (2010) Distanceindependent individualtree growth models of Larix olgensis plantation. J Northeast Forestry Univ 38:24–27 (in Chinese)
Mehtätalo L, deMiguel S, Gregoire TG (2015) Modelling heightdiameter curves for prediction. Can J For Res 45:826–837
Miina J, Pukkala T, Hotanen JP, Salo K (2010) Optimizing the joint production of timber and bilberries. For Ecol Manag 259:2065–2071
Mitscherlich EA (1919) Vorschriften zur Anstellung von Feldversuchen in der landwirtschaftlichen Praxis. P. Parey, Berlin 32 pp
Muchiri MN, Pukkala T, Miina J (2002) Optimising the management of maize—Grevillea robusta fields in Kenya. Agrofor Syst 56:13–25
Näslund M (1937) Skogsförsöksanstaltens gallringsförsök i tallskog (Forest research intitute’s thinning experiments in Scots pine forests). Meddelanden frstatens skogsförsöksanstalt Häfte 29 (in Swedish)
Niinimäki S, Tahvonen O, Mäkelä A, Linkosalo T (2013) On the economics of Norway spruce stands and carbon storage. Can J For Res 43:637–648
Olschewski R, Benítez PC (2010) Optimizing joint production of timber and carbon sequestration of afforestation projects. J For Econ 16:1–10
Palahí M, Pukkala T (2003) Optimising the management of Scots pine (Pinus sylvestris L.) stands in Spain based on individualtree models. Ann For Sci 60:95–107
Palahí M, Pukkala T, Bonet JA, Calinas C, Fischer CR, Martinez de Aragon JR (2009) Effect of the inclusion of mushroom values on the optimal management of evenaged pine stands of Catalonia. For Sci 55:503–511
PasalodosTato M (2010) Optimising forest stand management in Galicia, northwestern Spain. Dissertationes Forestales 102:1–52
PasalodosTato M, Pukkala T (2007) Optimising the management of evenaged Pinus sylvestris L. stands in Galicia, northwestern Spain. Ann For Sci 64:787–798
PasalodosTato M, Pukkala T, RigueiroRodríguez A, Fernández Núñez E, Mosquera Losada MR (2009) Optimal management of Pinus radiata silvopastoral systems established on abandoned agricultural land in Galicia (northwestern Spain). Silva Fenn 43:831–844
PasalodosTato M, Pukkala T, Calama R, Cañellas I, SánchesGonzález M (2016) Optimal management of Pinus pinea stands when cone and timber production are considered. Eur J For Res 135:607–619
Pihlainen S, Tahvonen O, Niinimäki S (2014) The economics of timber and bioenergy production and carbon storage in Scots pine stands. Can J For Res 44:1091–1102
Pingoud K, Pohjola J, Valsta L (2010) Assessing the integrated climatic impacts of forestry and wood products. Silva Fennica 44:155–175
Pohjola J, Valsta L (2007) Carbon credits and management of Scots pine and Norway spruce stands in Finland. Forest Policy Econ 9:789–798
Pukkala T (2009) Populationbased methods in the optimization of stand management. Silva Fennica 43:261–274
Pukkala T (2011) Optimising forest management in Finland with carbon subsidies and taxes. Forest Policy Econ 13:425–434
Pukkala T (2017) Optimal crosscutting: any effect on optimal stand management? Eur J For Res 136:583–595
Pukkala T, Miina J (1997) A method for stochastic multiobjective optimization of stand management. For Ecol Manag 98:189–203
Pukkala T, Lähde E, Laiho O, Salo K, Hotanen JP (2011) A multifunctional comparison of evenaged and unevenaged forest management in a boreal region. Can J For Res 41:851–862
Pukkala T, Lähde E, Laiho O (2014a) Stand management optimization—the role of simplifications. Forest Ecosystems 1:1–11
Pukkala T, Lähde E, Laiho O (2014b) Optimizing anyaged management of mixed boreal under residual basal area constraints. J For Res 25:627–636
Pukkala T, Lähde E, Laiho O (2015) Which trees should be removed in thinning treatments? Forest Ecosystems 2:2–12
Roise JP (1984) A nonlinear programming approach to stand optimization. For Sci 32(32):735–748
Roise JP (1986) An approach for optimizing residual diameter class distribution when thinning evenaged stands. For Sci 32:871–881
Solberg B, Haight RG (1991) Analysis of optimal economic management regimes for Picea abies stands using a stagestructured optimalcontrol model. Scand J For Res 6:559–572
State Forestry Administration (2014) The eighth forest resource survey report. http://211.167.243.162:8085/8/index.html
Valsta LT (1990) A comparison of numerical methods for optimizing even aged stand management. Can J For Res 20(7):961–969
Valsta L (1992) An optimization model for Norway spruce management based on individualtree growth models. Acta For Fenn 232:20 p
Vanclay JK (1994) Modelling forest growth and yield: applications to mixed tropical forests. CABI, Walingford 312 pp. ISBN 0851989136
Vauhkonen J, Pukkala T (2016) Selecting trees to be harvested based on the relative value growth of the remaining trees. Eur J For Res 135:581–592
Vettenranta J, Miina J (1999) Optimizing thinnings and rotation of Scots pine and Norway spruce mixtures. Silva Fenn 33:73–84
Yasaka M, Oono Y, Nakagawa M, Fukuchi M, Akashi N, Seiwa K (2008) Prediction of the diameter growth of individual trees in the larch plantation. The Bor For Soci 56:55–57
Yu B (2008) Artificial cultivation techniques of Larix olgensis. Heilongjiang Sci Tec Info 35:193 (in Chinese)
Zang H, Lei X, Zeng W (2016) Heightdiameter equations for larch plantations in northern and northeastern China: a comparison of the mixedeffects, quantile regression and generalized additive models. Forestry 89:434–445
Zhu Y, Dumroese RK, Pinto JR, Li GL, Liu Y (2013) Fall fertilization enhanced nitrogen storage and translocation in Larix olgensis seedlings. New For 44:849–861
Acknowledgments
The authors warmly thank Dr. Lihu Dong for his help in developing the individual tree growth model. The authors also thank the teachers and students of the Department of Forest Management, Northeast Forestry University (NEFU), China, who provided and collected the data for this study.
Funding
The research was financially supported by the National Key R&D Program of China (No. 2017YFD0600402), National Natural Science Foundation of China (31600511), and the Fundamental Research Funds for the Central Universities of the People’s Republic of China (2572017CA04).
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Handling Editor: Andreas Bolte
Contribution of the coauthors
Wei Peng developed the growth and yield model and wrote most of the paper. Timo Pukkala provided help in simulation and optimization and modified the paper. Xingji Jin helped to deal with the data.
Fengri Li supervised and coordinated the research, designed and installed the experiment, and contributed to the writing of the paper.
Electronic supplementary material
ESM 1
(RAR 29 kb)
Appendix. Optimization with the Hooke and Jeeves method
Appendix. Optimization with the Hooke and Jeeves method
When the method of Hooke and Jeeves (HJ) is used to optimize stand management, and the objective function is a utility function, the problem is formulated as follows:
where U is utility, calculated from the weights and values of the objective variables (net present value, wood production, carbon stock), and x is a vector of n real numbers which define the management schedule of the stand (thinning years, thinning intensity, rotation length, etc.). The elements of x are called decision variables. Different sets of values of decision variables are passed from the HJ optimization algorithm to a simulator, which uses these values to simulate stand development for one rotation and calculates the value of the utility function (Fig. 11). The utility value of the simulated management schedule is returned to the HJ algorithm. Based on the new and previous utility values, HJ alters the values of decision variables and passes the new values again to the simulation software. This process is repeated for many iterations.
The HJ method starts from an initial vector of decision variables (x_{1}). The search begins with exploratory search in which one decision variable is altered at a time (each decision variable corresponds to one coordinate direction of the search space). After completing one round of exploratory searches, the algorithm goes to pattern search if the exploratory search is successful (able to improve the solution). The pattern search makes simultaneous changes to more than one decision variable. If exploratory search cannot improve the solution, the step size used in exploratory search is divided by two and the search is repeated. The search is terminated once the step size becomes smaller than a predefined stopping criterion.
Below is an exact description of the HJ algorithm (the description is based on Bazaraa et al. (1993) and Pukkala (2009)). Let d_{1},…,d_{m} be the coordinate directions and let y_{1} = x_{1} and k = j = 1.

Exploratory search

1.
If f(y_{j} + Δd_{j}) > f(y_{j}) (success), let y_{j + 1} = y_{j} + Δd_{j} and go to Step 2. Otherwise if f(y_{j}Δd_{j}) > f(y_{j}) (success), let y_{j + 1} = y_{j}Δd_{j} and go to Step 2. Otherwise, let y_{j + 1} = y_{j} and go to Step 2 (no improvement found in direction j).

2.
If j < m, replace j by j + 1 and repeat Step 1 (go to next decision variable). Otherwise, go to Pattern search if f(y_{m + 1}) > f(x_{k}) (at least one successful change detected in the directions of coordinate axes). If f(y_{m + 1}) ≤ f(x_{k}), go to Step size reduction.

1.

Pattern search

3.
Let x_{k + 1} = y_{m + 1} and let y_{1} = x_{k + 1} + α(x_{k + 1}x_{k}). Replace k by k + 1, let j = 1, and go to Step 1.

3.

Step size reduction

4.
If Δ ≤ ε, stop. Otherwise, replace Δ by Δ/2. Let y_{1} = x_{k}, x_{k + 1} = x_{k}, j = 1, replace k by k + 1, and repeat Step 1

4.
The method has three parameters: initial steps size (Δ), stopping criterion (ε), and α. In this study, the initial step size was different for different decision variables, and it was equal to 0.1 times the range of variation allowed for the decision variable. The stopping criterion was equal to 0.01 times the initial step size. The search was terminated once the step size was smaller than the stopping criterion for every decision variable. Parameter α was taken as 1.0.
Rights and permissions
About this article
Cite this article
Peng, W., Pukkala, T., Jin, X. et al. Optimal management of larch (Larix olgensis A. Henry) plantations in Northeast China when timber production and carbon stock are considered. Annals of Forest Science 75, 63 (2018). https://doi.org/10.1007/s1359501807391
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s1359501807391
Keywords
 Growth and yield model
 Net present value
 Multifunctional forestry
 Multiobjective optimization