Optimal management of larch (Larix olgensis A. Henry) plantations in Northeast China when timber production and carbon stock are considered

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³ ha⁻¹ a⁻¹, 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 fast-growing and high-yielding 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³ ha⁻¹ for SI 12, 8.2 m³ ha⁻¹ for SI 16, and 11.7 m³ ha⁻¹ 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 fast-growing and high-yielding 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³, 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 sub-alpine areas of China (Chen 2010). It plays an important role in the establishment of fast-growing and high-yielding 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 multi-objective 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 multi-objective optimization of forest management, mainly focusing on the joint production of wood and non-wood 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 (Pasalodos-Tato 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 population-based 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 stand-level models, diameter distribution models, and individual tree models. Individual-tree models may be divided into distance-dependent and distance-independent 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 stand-level growth models whereas nonlinear programming methods can use also other types of growth models. For example, Brodie et al. (1978) used dynamic programming with stand-level growth models, de Miguel et al. (2014) employed nonlinear programming with distance-independent individual-tree growth models, and Muchiri et al. (2002) used nonlinear programming with distance-dependent individual-tree growth models. Since individual-tree 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 optimization-based 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, individual-tree 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 individual-tree 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 5-year 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²).

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.

Table 1 Description of the dataset (205 sample plots) used to model dominant height development

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 5-year intervals. Only one 5-year interval per plots was used for modeling. The number of individual-tree observations for modeling the 5-year diameter increment was 7114, and 7838 observations were available for modeling the 5-year survival (Table 2).

Fig. 1

Left panel: Distribution of DBH and height for trees that were used for height modeling. Right panel: Quadratic mean diameter and average height on plot level for the plots that were used for modeling diameter increment, survival, and tree height

Table 2 Summary for the dataset (148 sample plots) that was used for modeling individual-tree diameter increment, survival, and height

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

Table 3 Descriptive statistics of the dataset (75 trees) used for taper and biomass modeling

Of the 75 trees felled for the taper modeling, 68 trees were also measured for biomass (Table 3). The stems were cut into 1-m sections and each section was weighed in the field. A 2–3-cm-thick 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 oven-dried 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):

$$ {HD}_{quide}={b}_1\left(1-\exp \left({b}_2T\right)\right) $$

(1)

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₅) and beginning (d₀) of a 5-year period was used:

$$ \ln \left({d}_5^2-{d}_0^2\right)=f\left(d,G, BAL, SI\right) $$

(2)

Binary logistic regression analysis was used to model the probability that a tree will survive for the coming 5-year 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:

$$ s=\frac{1}{1+\exp \left(-\left(f\left(d,G, BAL, SI\right)\right)\right)} $$

(3)

The Näslund function (Näslund 1937) was selected as the height-diameter 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:

$$ h=1.3+\frac{d^2}{{\left(a\times d+b\right)}^2} $$

(4)

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

$$ h=1.3+\frac{d^2}{{\left(\left({a}_1+{a}_2\ln \kern0.5em HD\right)\times d+{b}_1+{b}_2\ln \kern0.5em HD\right)}^2} $$

(5)

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

$$ {d}_i={b}_0{d}^{b_1}{h}^{b_2}{\left[\frac{1-{q}^{1/3}}{1-{t}^{1/3}}\right]}^{\left({b}_3{q}^4\left[1/{e}^{d/h}\right]+{b}_5{\left(\frac{1-{q}^{1/3}}{1-{t}^{1/3}}\right)}^{0.1}+{b}_6\left(1/d\right)+{b}_7{h}^{1-{q}^{1/3}}+{b}_8\left(\frac{1-{q}^{1/3}}{1-{t}^{1/3}}\right)\right)} $$

(6)

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₀–b₈ 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:

$$ \left\{\begin{array}{l}{w}_r=\exp \left({\beta}_{10}\right)\cdot {d}^{\beta_{11}}+{\varepsilon}_r\\ {}{w}_s=\exp \left({\beta}_{20}\right)\cdot {d}^{\beta_{21}}+{\varepsilon}_s\\ {}{w}_b=\exp \left({\beta}_{30}\right)\cdot {d}^{\beta_{31}}+{\varepsilon}_b\\ {}{w}_f=\exp \left({\beta}_{40}\right)\cdot {d}^{\beta_{41}}+{\varepsilon}_f\\ {}{w}_t=\exp \left({\beta}_{10}\right)\cdot {d}^{\beta_{11}}+\exp \left({\beta}_{20}\right)\cdot {d}^{\beta_{21}}+\exp \left({\beta}_{30}\right)\cdot {d}^{\beta_{31}}+\exp \left({\beta}_{40}\right)\cdot {d}^{\beta_{41}}+{\varepsilon}_t\end{array}\right. $$

(7)

where w is the oven-dry 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 β_i0 and β_i1 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⁻¹, 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.

Table 4 Characteristics of initial stands

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.

Table 5 Timber assortments and prices used in the study

2.6 Objective function

Both single-objective and multi-objective optimizations were used in this study. The following utility functions were used in the multi-objective optimization:

$$ {MOF}_1=0.45\left( NPV/{NPV}_{\mathrm{max}}\right)+0.45\left( WP/{WP}_{\mathrm{max}}\right)+0.1\left(C/{C}_{\mathrm{max}}\right) $$

(8)

$$ {MOF}_2=0.4\left( NPV/{NPV}_{\mathrm{max}}\right)+0.4\left( WP/{WP}_{\mathrm{max}}\right)+0.2\left(C/{C}_{\mathrm{max}}\right) $$

(9)

$$ {MOF}_3=0.333\left( NPV/{NPV}_{\mathrm{max}}\right)+0.333\left( WP/{WP}_{\mathrm{max}}\right)+0.333\left(C/{C}_{\mathrm{max}}\right) $$

(10)

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⁻¹). Similarly, WP_max (m³ ha⁻¹ a⁻¹) and C_max (t ha⁻¹) 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₂ 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₁ and MOF₃) 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₁ of the thinning intensity curve (see below)

• Parameter a₂ 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):

$$ TI(d)=\frac{1}{1+\exp \left({a}_1\left({a}_2-d\right)\right)} $$

(11)

where TI(d) is the proportion of harvested trees when DBH is d cm; and a₁ and a₂ are parameters to be optimized. Parameter a₂ gives the diameter at which thinning intensity is 0.5, and a₁ defines the type of thinning. When a₁ is negative, the thinning is low thinning (thinning from below), a₁ > 0 corresponds to high thinning (thinning from above), and a₁ = 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; Pasalodos-Tato 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):

$$ {HD}_{guide}=30.4965\left(1-\exp \left(-0.0332\times T\right)\right) $$

(12)

Fig. 2

The dominant height and site index model. The black dots are the data points for fitting the guide curve model (dotted line). The solid lines with different colors give the dominant height development for site indices 12–24 m

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:

$$ SI=\frac{HD(T)}{HD_{\mathrm{guide}}(T)}\times {HD}_{\mathrm{guide}}(30) $$

(13)

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

$$ HD(T)=\frac{SI}{HD_{\mathrm{guide}}(30)}\times {HD}_{\mathrm{guide}}(T) $$

(14)

After fitting, we obtained the following model for diameter increment:

$$ \ln \left({d}_5^2-{d}_0^2\right)=y=1.945+0.995\ln {d}_0-0.018{d}_0+0.028 SI-0.023G-0.082 BAL/\ln \left({d}_0+0.01\right) $$

(15)

$$ {I}_d={d}_5-{d}_0=\sqrt{\ \exp (y)+{d}_0^2}-{d}_0 $$

(16)

where I_d is the 5-year diameter increment (cm). The RMSE of the model was 0.41 and its R² 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.

Fig. 3

Relationships between DBH, stand basal area (G), site index (SI), the basal area of larger trees (BAL), and 5-year diameter increment

The model for the probability that a tree survives for 5 years (s) was as follows:

$$ s=\frac{1}{1+\exp \left[-\left(-3.9749+2.8665\ln d+0.1131 SI-0.2782 BAL/\ln \left(d+0.01\right)\right)\right]} $$

(17)

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.

Fig. 4

The relationships between DBH, site index (SI), the basal area of larger trees (BAL), and 5-year survival probability. Note that the y axes of the diagrams are cut

The following height-diameter model was fitted to calculate the individual tree height:

$$ h=1.3+\frac{d^2}{{\left(\left(0.4520-0.0817\;\ln HD\right)\times d+1.5310-0.2988\kern0.22em \ln HD\right)}^2} $$

(18)

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.

Fig. 5

Tree height-diameter curve for different dominant heights (Hdom)

The DBH-height curve “climbed up” when the stand developed and its dominant height increased.

The Kozak model (Kozak 2004) fitted for stem taper was as follows:

$$ {d}_i=0.8753{d}^{0.9939}{h}^{0.0551}{\left[\frac{1-{q}^{1/3}}{1-{t}^{1/3}}\right]}^{\left(0.5350{q}^4-0.2409\left[1/{e}^{d/h}\right]+0.3399{\left(\frac{1-{q}^{1/3}}{1-{t}^{1/3}}\right)}^{0.1}+0.801\left(1/d\right)-0.0091{h}^{1-{q}^{1/3}}+0.0811\left(\frac{1-{q}^{1/3}}{1-{t}^{1/3}}\right)\right)} $$

(19)

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.

Table 6 Parameters of biomass models and carbon content of each component for larch plantations

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 multi-objective utility function (MOF₂) 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.

Fig. 6

Effect of number of thinnings and discount rate (1%, 2%, 3%) on relative utility (MOF₂) in SI 16 m for initial stand densities 2800 trees/ha (top) and 2125 trees/ha (bottom). The highest utility with a certain discount rate has relative utility equal to 1 and the other utilities are relative to the highest value. Note that the y axis is cut at 0.85

3.3 Effect of management objective on optimal management

The results of the single-objective (NPV or wood production) and multi-objective 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⁻¹when WP was maximized instead of NPV. Maximizing the baseline multi-objective utility function (MOF₂) 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 multi-objective utility function extended the rotation length.

Table 7 Optimization results of single-objective (NPV, WP) and multi-objective (MOF_1–3) optimization with three thinnings for three different site indices when the stand density at 10 years is 2800 trees/ha

Table 8 Optimization results of single-objective (NPV, WP) and multi-objective (MOF_1–3) optimization with two thinnings for three different site indices when the stand density at 10 years is 2125 trees/ha

The trade-offs 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 trade-off. Maximizing NPV led to a 30–35% decrease in the average carbon stock of the forest, when compared to maximizing MOF₃, where the weight of carbon stock was 0.33. The baseline objective function, MOF₂, was good in all other outputs except carbon stock, where it was 5–10% worse than in the management schedules that maximized MOF₃.

Fig. 7

Sensitivity of the amounts of different outputs to the objective function used in optimization. The highest amount of each output gets a value equal to 1, and the amounts of outputs in the other solutions are expressed as proportions of the highest value

When MOF₂ 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 multi-objective schedule was a few years longer than in the schedules that maximized NPV.

Fig. 8

Development of growing stock volume for SI 16 m and initial stand density 2800 trees/ha when multi-objective utility function (MOF₂) is maximized. The right-hand-side diagram shows the optimal harvest percentage in the first, second, and third thinning for different DBH classes

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.

Table 9 Recommended management models for larch plantations

Fig. 9

Development of stand volume in the optimal management schedule (for MOF₂) in different site indices for two initial stand densities (2800 or 2125 trees/ha at 10 years)

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

Fig. 10

Impact of management constraints on NPV, wood production, carbon stock, and utility when the objective function is MOF₂ and the initial stand density is 2800 trees/ha. “Max thinning intensity 45%” indicates that the thinning intensity could not exceed 45% of stand basal area, and “Forced low thinning” indicates that the maximum thinning intensity is 45% and only thinning from below is allowed

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 large-sized timber assortments. Large-sized 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 long-term carbon stores. A gradual shift towards longer rotations, higher growing densities, and saw log-oriented management would increase the carbon stocks for many decades, both in forests and in different wood-based 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; Pasalodos-Tato and Pukkala 2007; Pasalodos-Tato 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 5-year 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 science-based 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 clear-cutting 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 clear-cutting 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 medium-sized 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 medium-sized logs. Chen (2010) concluded that site index 16 m is suitable for the production of medium-sized 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 large-diameter 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 multi-objective 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 multi-objective 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.

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:

$$ {\displaystyle \begin{array}{l}\operatorname{Max}\kern0.5em U=f\kern0.1em \left(\mathbf{x}\right)\\ {}\mathrm{subject}\ \mathrm{to}\\ {}\mathbf{x}\in {\mathrm{\mathbb{R}}}^n\end{array}} $$

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.

Fig. 11

Optimization of stand management employs two software components: optimization algorithm and simulation program

The HJ method starts from an initial vector of decision variables (x₁). 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₁,…,d_m be the coordinate directions and let y₁ = x₁ 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.

• Pattern search

  • 3. Let x_k + 1 = y_m + 1 and let y₁ = x_k + 1 + α(x_k + 1-x_k). Replace k by k + 1, let j = 1, and go to Step 1.

• Step size reduction

  • 4. If Δ ≤ ε, stop. Otherwise, replace Δ by Δ/2. Let y₁ = x_k, x_k + 1 = x_k, j = 1, replace k by k + 1, and repeat Step 1

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.