How to maximize the joint benefits of timber production and carbon sequestration for rural areas? A case study of larch plantations in northeast China

Abstract

Background

Implementing large-scale carbon sink afforestation may contribute to carbon neutrality targets and increase the economic benefits of forests in rural areas. However, how to manage planted forests in China to maximize the joint benefits of timber production and carbon sequestration is still unclear. Therefore, the present study quantified the effects of different rotation lengths, thinning treatments, site quality (SCI), stand density (SDI), and management costs on the joint benefits of carbon sequestration and timber production based on a stand-level model system developed for larch plantations in northeast China.

Results

The performances of the different scenarios on carbon stocks were satisfactory, where the variations in the outcomes of final carbon stocks could be explained by up to 90%. The joint benefits increased significantly with the increases of SDIs and SCIs, regardless of which rotation length and thinning treatments were evaluated. Early thinning treatments decreased the joint benefits significantly by approximately 131.53% and 32.16% of middle- and higher-SDIs, however longer rotations (60 years) could enlarge it by approximately 71.39% and 80.27% in scenarios with and without thinning when compared with a shorter rotation length (40 years). Discount rates and timber prices were the two most important variables affecting joint benefits, while the effects of carbon prices were not as significant as expected in the current trading market in China.

Conclusions

The management plans that promote longer rotations, higher stand densities, and no thinning treatments can maximize the joint benefits of carbon sequestration afforestation and timber production from larch plantations located in northeast China.

Graphical Abstract

Highlights

• A compatible model of stand volume and carbon for larch plantations was developed.

• Thinning decreased the joint benefits of higher SDIs, but longer rotations enlarged it.

• Discount rates and timber prices predominated the joint benefits.

• The effects of carbon prices on the joint benefits were unconspicuous.

Background

Global surface temperatures have increased by approximately 0.6 ℃ during the last few decades and will continue to grow until at least mid-century under all emissions scenarios considered [1]. Thus, deep reductions in CO₂ emissions are essential actions that may keep global warming well below 2 ℃ or the more aggressive goal of 1.5 ℃. As a part of nature-based solutions for mitigating climate change, implementing sustainable forest management has drawn much attention from policymakers, forest landowners, and other stakeholders. Planted forests may not only provide a large and persistent wood and non-wood product supply for human society, but also they may also absorb a large amount of CO₂, thus contributing significantly to mitigating global climate change. It has been estimated that China might be the most promising Clean Development Mechanism (CDM) market [2], which may account for about 40–50% of the global market. Thus, the carbon trading market under the framework of China Certified Emission Reduction (CCER) started again in 2021; this market is considered an essential route to achieving carbon–neutral goals. The world's total planted forest area approached about 795.43 million hectares in 2020, representing a quarter of the global greening increases [3]. However, how to manage these forests to maximize the joint benefits of timber production and carbon sequestration is still unclear.

Forest growth simulators are powerful tools that assist in developing valuable guidelines for implementing sustainable forest management. Similar to the development of estimates of forest stand (a collection of similar trees) volume, stand-level carbon stocks are also a function of site quality, development phase, stand density, and management treatments applied. The most popular methods for estimating stand carbon stocks utilize a biomass conversion factor (BCF), biomass expansion factors (BEF), and stand volumes. However, the commonly used BCFs and BEFs are often fixed values for specific tree species, and the estimates of carbon stocks are unsatisfactory [4,5,6]. To overcome the problems encountered, some continuous variable BCF and BEF models have recently been developed [7, 8]. A second strategy to estimate stand-level carbon stocks involves combining commonly used tree- [7, 9] or stand- [4, 8] level biomass models with a detailed forest inventory. The variables in tree-level models include diameter at breast height (DBH), tree height (HT), and combinations of these (e.g., DBH²HT). The variables in stand-level models commonly include mean DBH, mean HT, stand basal area (BAS) and stand age. However, the second strategy can only estimate the carbon stocks at a particular point in time, and do not provide a series of stand carbon stocks with stand age. In addition, some forest stand models (e.g., FORMIND) [10] or forest landscape models (e.g., LANDIS-II) [11] can also be used to simulate the development of forest ecosystems. However, some uncommonly used terminologies (e.g., “patch” in FORMIND and “cohort” in LANDIS-II) and uneasy measured parameters may prevent foresters from using them properly. Thus, developing suitable prediction models of stand stocks is urgently needed for the development of the forest carbon sequestration projects.

Assigning suitable management plans to planted forests to maximize the joint benefits of carbon and timber is one of the essential tasks in the carbon sequestration afforestation project under the framework of either CDM or CCER. The length of optimal rotation periods (the time from tree planting to final harvest) for planted forests has been discussed in previous studies, especially when the benefits of carbon sequestration were considered. However, the conclusions of these works were not always consistent. For example, Hoel et al. [12] reported that a social cost of carbon implied longer optimal rotation periods, and that if the social cost of carbon exceeded a specific threshold value, the forest should not ever be harvested. However, Zhou and Gao [13] and Dong et al. [14] indicated that the effects of increasing carbon prices on optimal rotation lengths were not as remarkable as expected. The reasons may be that an age-dependent volume growth model and a fixed BCF and BEF were used by Hoel et al. [12], while a set of stand-level growth and yield models were employed by Dong et al. [14], where the derived BCF and BEF were both alterable. Meanwhile, carbon prices, discount rates, timber prices [12, 15, 16], and carbon accounting strategies [14, 17], may significantly affect the optimal management plans. Thus, managing carbon sequestration afforestation projects is a very complex issue that must be approached with caution.

Larch species, which include Larix gmelinii, Larix olgensis, and Larix kaempferi, are some of the most critical planted tree species in northeast China. The ninth National Forest Inventory of China reported that the current areas of larch plantations accounted for approximately 31,630 km² of land [18]. Some management techniques, such as cultivations for large-sized and pulp timber, have been put forward [19,20,21], while composite management technologies involving carbon sequestration and timber production are still highly lacking. Thus, this study focused on (1) developing compatible stand volume and carbon stocks of larch plantations in northeast China; (2) simulating different thinning treatments and rotation lengths on the joint benefits of timber production and carbon sequestration; (3) analyzing the sensitivity of different costs and benefits on the profitability of a carbon sequestration afforestation project.

Materials and methods

Study area

This study was conducted in the three provinces of northeast China (Fig. 1; 118° 50′–135° 05′ E, 38° 43′–43° 25′ N), namely Heilongjiang, Jilin, and Liaoning. These three provinces cover an area of approximately 787,300 km², and forest coverage is about 42.4%. The average elevation of the region is 320 m, and weather conditions are dominated by a temperate continental monsoon climate, where the mean annual temperature is about 3℃, and the mean annual precipitation is about 650 mm. Soil conditions in this region are predominantly composed of a mountain of dark brown soil, accompanied by a small amount of meadow and white pulp soil. Most of these planted forests are located in rural areas, where the mean annual incomes of forest dwellers are only about half of that for non-forest residents (about 16,400 RMB/person≈2300 $/person), thus improving the incomes and well-being of forest dwellers is a worthwhile endeavor. Therefore, developing carbon sequestration afforestation projects have received widespread attention from forest dwellers [22].

Fig. 1

The locations of study area (A) and the age frequency (B) and spatial distribution (C) of sample plots

Field data

The data used in this study were obtained from field plots associated with the 7th National Forestry Inventories (NFI) of China in the three provinces. The size of the square field plots was 0.067 ha. In each plot, all trees with a diameter at breast height (DBH) larger than 5 cm were measured. Meanwhile, the heights of at least three intermediate (crown position) trees were measured using the ultrasonic altimeter to calculate the mean stand height. The stand ages were determined by consulting the afforestation archives or counting tree rings from a core extracted at breast height. Other variables measured included stand (e.g., canopy density, historical disturbances) and topography (e.g., elevation, aspect, slope, and soil types) characteristics. To avoid complications that may arise during the analysis of larch forests, particularly where other tree species may regenerate naturally within the plots, only larch plantation plots where the larch component of volume was larger than 65% were selected for this analysis [23]. In total, 342 plots were obtained across the three provinces: 90 plots in Heilongjiang, 202 plots in Jilin, and 50 plots in Liaoning (Fig. 1).

Calculating forest carbon stocks

For larch trees, the carbon stocks of each component (e.g., root, stem, branch, and leaf) were jointly calculated using the compatible biomass models of Peng et al. [24] and the component-specific carbon concentration [25], namely 0.4617 of root, and 0.4610 of stem, and 0.4736 of branch, and 0.4734 of leaf, respectively. The carbon stocks of other species (e.g., Betula platyphylla, Populus davidiana, and Fraxinus mandshurica) were estimated using a similar process. However, the biomass models were extracted from the work of Dong [26], and the carbon concentrations were a fixed value (0.5) for all four components. Finally, the total carbon stocks of entire plots were calculated using the accumulation stocks of all the trees and the area of the plots, and the values were expanded to a per-hectare value. In forestry, the stand density index (SDI) can combine the effects of stand density and mean DBH together, while the site class index (SCI) can also synthesize the combined impacts of terrain, soil, climate, and other factors on forested land. Since SDI and SCI have huge effects on the final yields (Fig. 2), they have been widely used in forest growth and yield models [17, 27, 28]. Thus, both SDI and SCI of each plot were further calculated using the functions from the work of Wang [29], which were fitted using 1140 plots.

$$\text{SDI}=N\cdot {\left(\frac{{D}_{0}}{{D}_{g}}\right)}^{-1.605}$$

(1)

$$\text{SCI}=\text{HT}\frac{{(1-Exp(-0.0231{t}_{0}))}^{0.8365}}{{(1-Exp(-0.0231t))}^{0.8365}}$$

(2)

where N is the number of trees per hectare; Dg is the mean DBH of the stand; ${t}_{0}$ and ${D}_{0}$ are respectively the base age and base DBH of larch plantation, which was set 30 years and 20 cm in this analysis [29]. Further, HT is the mean height of the stand, and t is the actual current age of the trees in each plot. The descriptive statistics of the stand variables used for modeling are shown in Table 1.

Fig. 2

The relationships between stand age and mean stand height (A), and between stand age and stand basal area (B), and between mean stand height and stand carbon stocks (C), and between stand density index and stand basal area (D), and between stand basal area and stand carbon stocks (E), and between stand volume and stand carbon stocks (F), where the blue lines represent the smoothed conditional means, and the grey ribbons represent the 95% confidence interval

Table 1 Descriptive statistics of the stand variables for larch plantations in northeast China, where DBH and HT represent the diameter at breast height and tree height, and Std and CV represented the variable's standard deviation and variation coefficient, respectively

Carbon stock prediction models

Stand basal area (BAS) and mean stand height (HT) are not only two critical factors that significantly affect the development of stand carbon stocks, but they are also straightforward to measure in practice. Therefore, BAS and HT were used as bridges to estimate stand carbon stocks of larch plantations. Scatter plots indicate that both BAS and HT increase significantly with ages in the early stages, and then the increments decrease gradually with further increases in stand ages (Fig. 2), which suggests a typical monotonic increasing tendency. Thus, the Mitscherlich function was selected to simulate the growth processes of BAS and HT.

$$y=a\left(1-\text{Exp}(-bt)\right)$$

(3)

where y is the dependent variable, namely either BAS or HT in this analysis; Exp() is an exponential function with the base of natural numbers; a and b represent the potential maximum values and growth rates of BAS or HT; and t represents the current age of the trees in each plot.

Numerous studies have demonstrated that the potential maximum values of HT highly depend on site quality, whereas the effects of stand densities and management treatments were both inconspicuous [30,31,32]. Thus, the parameter a of Eq. 3 was associated with SCI using a power form.

$$HT={a}_{0}{SCI}^{{a}_{1}}\left(1-\text{Exp}(-{a}_{2}t)\right)$$

(4)

where${a}_{0}$, ${,a}_{1}$ and ${a}_{2}$ are the estimated parameters. BAS is a composite variable between the mean DBH and the number of trees per hectare (N), thus both stand density and site quality have significant effects on the growth rates of BAS, namely parameter b of Eq. 3. Since the effects of SCI on HT have been quantified in Eq. 4, it was not considered again in the model of BAS. The function of BAS was formulated as:

$$\text{BAS}={b}_{0}\left[1-Exp(-{b}_{1}{\left(\frac{SDI}{1000}\right)}^{{b}_{2}}t)\right]$$

(5)

where b₀, b₁, and b₂ are all the estimated parameters. Based on BAS and HT, the stand volume could then be estimated as [28, 33]:

$$\text{VOL}=\text{BAS}\left[\frac{{c}_{0}HT}{(HT+{c}_{1})}\right]$$

(6)

where ${c}_{0}$ and ${c}_{1}$ are the estimated parameters. Since significant linear relations between stand volume and stand carbon stocks could be observed (Fig. 2), they were thus formulated as:

$$\text{CAR}={d}_{0}+{d}_{1}VOL={d}_{0}+{d}_{1}\left\{\text{BAS}\left[\frac{{c}_{0}HT}{(HT+{c}_{1})}\right]\right\}$$

(7)

where ${d}_{0}$ and ${d}_{1}$ are the estimated parameters.

However, one thing to note was that the variables of HT, BAS, and VOL were not only the dependent variables of their respective models (namely Eqs. 4, 5, and 6), but also the independent variables of Eq. 7; thus, they were solved as a set of simultaneous equations.

$$\left\{\begin{array}{c}HT={a}_{0}{SCI}^{{a}_{1}}\left(1-\text{Exp}(-{a}_{2}t)\right)\\ \begin{array}{c}\text{BAS}={b}_{0}\left[1-Exp(-{b}_{1}{\left(\frac{SDI}{1000}\right)}^{{b}_{2}}t)\right]\\ \text{VOL}=\text{BAS}\left[\frac{{c}_{0}HT}{(HT+{c}_{1})}\right]\end{array}\\ CAR={d}_{0}+{d}_{1}VOL\end{array}\right.$$

(8)

Due to the characteristics mentioned above of simultaneous equations, several methods, e.g., seemingly unrelated regression (SUR), two-stage least squares (2SLS), and three-stage least squares (3SLS), have been put forward to estimate the parameters synchronously [7, 8]. SUR, proposed by Parresol [34], has been widely used in recent years [9, 35]. Thus, the method of SUR embedded in the R package of “systemfit” was employed to estimate the parameters of the simultaneous equations (Eq. 8).

Since developing a stand density prediction model requires at least two periods of monitoring data, which are not available in our study, we decided to use a published model for this species by Chen (2010), who developed it using a total of 285 field measurement plots that were visited five times (1986, 1990, 1995, 2000, and 2005). The coefficient of determination for this model was as significant as about 0.9223.

$${N}_{2}={N}_{1}\text{exp}\left(-\left(0.0103+0.0003SCI\right)\left({t}_{2}-{t}_{1}\right)\right)$$

(9)

where N₁ is the trees per hectare at time t₁, and N₂ is the trees per hectare at time t₂; SCI is the site condition index. The coefficients within the functions indicate that the number of trees per hectare decrease gradually with increased SCIs, which is consistent with biological knowledge [27].

The models' accuracies and performances were evaluated using the k-fold cross-validation method. This approach involves randomly dividing the set of observations into k folds of approximately equal size. Then, the first fold was treated as a validation set, while the remaining k-1 folds were used to fit the models [36]. The processes of fitting and validating would be repeated k times, and the final results of k-fold cross-validation would be summarized with the means of various statistics on the model. Following the suggestions of James et al. [36] and Zeng et al. [7], the value of k was set as 5 in this analysis. The statistics used in this analysis included the adjusted coefficient of determination (${R}_{a}^{2}$; Eq. 10), root mean square error (RMSE; Eq. 11) and mean absolute percent error (MAPE; Eq. 12), respectively.

$${R}_{a}^{2}=1-\frac{{\sum }_{i=1}^{N}{({y}_{i}-{\widehat{y}}_{i})}^{2}}{{\sum }_{i=1}^{N}{({y}_{i}-\overline{y })}^{2}}\cdot \left(\frac{N-1}{N-p-1}\right)$$

(10)

$$\text{RMSE}=\sqrt{\frac{{\sum }_{i=1}^{N}{({y}_{i}-{\widehat{y}}_{i})}^{2}}{N-p}}$$

(11)

$$\text{MAPE}=\frac{1}{n}\sum_{i=1}^{N}\left|\frac{{y}_{i}-{\widehat{y}}_{i}}{{y}_{i}}\right|\times 100\%$$

(12)

where ${y}_{i}$ and ${\widehat{y}}_{i}$ are respectively the measured- and predicted-values for i-th plot; $\overline{y }$ is the mean values of response variable; p is the number of model parameters.

Application of carbon stock models

Following the frameworks of CDM and CCER [37], the monitoring period of afforestation projects for carbon sequestration should be fixed at 5 years after the first monitoring time is determined. Thus, we assumed that each of the commitment periods was five years long to be consistent, namely the endings of the corresponding commitment periods were respectively 5th, 10th, 15th, etc. Based on the conservative principles of CDM, carbon sequestration considered in this analysis is only related to the carbon stocks of above- and below-ground living biomass. This analysis did not consider carbon stocks related to soil organic matter, dead trees, and litter. The costs and benefits of a carbon sequestration afforestation project could then be determined using the net present value (NPV) of forest management activities.

$${NPV}_{5t}=\sum_{t=1}^{N}\frac{{W}_{5t}+{{C}_{5t}-R}_{5t}}{{\left(1+r\right)}^{5t}}$$

(13)

where t is the t-th 5-year period; N is the number of total 5-year periods, which was calculated as N = T/5; T is the rotation length of a larch plantation; ${NPV}_{5t}$ is the joint net present value of timber and carbon together for the t-th 5-year period; ${R}_{5t}$ is the present value of the management costs for t-th 5 years; ${W}_{5t}$ and ${C}_{5t}$ are respectively the present values on the benefits of the t-th 5-year periods for timber and carbon; r is the discount rate.

Following the project design document (PDD) of a carbon sequestration afforestation project for Shibazhan Forestry Bureau that was certificated in 2019 (https://www.ccchina.org.cn/, Accessed on Dec 23, 2023), the afforestation costs, which included the site preparation, seedling, and tending, were $280 ha⁻¹. The annual maintenance costs, including activities to prevent pests and fires, were $12 ha⁻¹ yr⁻¹. The net prices of commercial wood were assumed to be $120 m⁻³. The certification costs, which included monitoring, validation, and certification, were assumed to be $10 per hectare per period, and incurred every 5 years. The carbon prices were extracted from the history data of Carbon Trading Network (http://www.tanpaifang.com/, Accessed on Dec 23, 2023), and were assumed to average $5 ton⁻¹, while the discount rates followed the commonly used values (namely 3%) in forestry in China. To further quantify the sensitivity of different management costs and benefits on the joint NPV of timber and carbon together for a carbon sequestration afforestation project, the benefits and costs mentioned above were increased or decreased by 50%, respectively.

Since both SCIs and SDIs have significant effects on the growth and final yield of carbon stocks (Fig. 2), three different classes on the site quality of SCI and the stand density of SDI that were inspired by the statistics in Table 1 were employed in this analysis: 100, 300 and 500 trees ha⁻¹ of SDI, and 10 m, 14 m and 18 m of SCI, respectively. The harvest plans for larch plantations followed the Technical Schedule of Fast-growing Larch Plantations [19], where the thinning treatments were assigned at the 17th, 25th, and 33rd years, and with a volume-based intensity of 18.0%, 34.1%, and 25.0%, respectively. Since the minimum rotation lengths of larch plantations were 40 years in northeast China, two crediting periods were considered in this analysis. However, some studies, e.g., Dong et al. [14], Holtsmark et al. [15], and Hoel et al. [12], have argued that the rotations could be extended significantly, especially when the benefits of carbon sequestration were considered, thus lengthening the operating period from two (40 years) to three (60 years) crediting periods would also be simulated in this analysis. The average outturn percentages of commercial wood were assumed to be 30% for thinning treatment following the Technical Regulation of Commercial Timber Ratio [38]. In comparison, they were considered to be 70% for final clearcutting. The dynamics of various stand variables under different scenarios were simulated using the developed models following the procedure described in the Appendix.

Results

Accuracy of carbon stock models

The simultaneous equations' parameter estimates differed significantly from zero (p < 0.01; Table 2). The goodness-of-fits indicated all four models were fitted well with the measured values, during which the ${R}_{a}^{2}$ of HT model was the largest (${R}_{a}^{2}$=0.9993), while the lowest was observed for BAS model (${R}_{a}^{2}$=0.8800). The values of RMSE for HT, BAS, VOL, and CAR were 0.1046 m, 1.1746 m²ha⁻¹, 1.7400 m³ ha⁻¹, and 0.5640 ton ha⁻¹, respectively, which accounted only for approximately 0.89%, 12.91%, 2.69% and 2.88% of the average amounts of these. The MAPEs of the four models were all relatively small, during which the largest was observed for the BAS model (15.36%), and the smallest was found for the HT model (0.85%).

Table 2 Parameter estimated values and goodness-of-fit values for the simultaneous equations (Eq. 8), where HT, BAS, VOL, and CAR represent the mean stand height, stand basal area, stand volume, and stand carbon stocks; ${R}_{a}^{2}$, RMSE and MAPE are the adjusted coefficient of determination, and root mean square error, and mean absolute percent error, respectively

The predicted values of stand volume and carbon stocks were highly correlated with their observations (Fig. 3), where the slopes were both near 1.0 (namely 1.0037 and 1.0033 of stand volume and carbon stocks), and the correlation coefficients were also as large as 0.98. The five-fold cross-validation statistics also highlighted that the accuracy and performance of the four models were relatively higher (Table 3). Therefore, these models could meet the carbon sequestration afforestation project requirements.

Fig. 3

Predicted vs measured values of stand volume and carbon stocks for larch plantations in northeast China, where the red dashed lines represent the 1:1 line between measured- and predicted-values of stand volume and carbon stocks, respectively

Table 3 Statistics regarding the five-fold cross-validations for the simultaneous equations (Eq. 8), where HT, BAS, VOL, and CAR represent the mean stand height, stand basal area, stand volume, and stand carbon stocks; ${R}_{a}^{2}$, RMSE and MAPE are the adjusted coefficient of determination and root mean square error, and mean absolute percent error, respectively

Simulation of carbon stock developments

Based on the estimated parameters in Table 2, the development of stand carbon stocks and carbon sequestrations with stand ages were simulated for different combinations between SCIs and SDIs (Fig. 4), which highlighted that SCIs have considerable effects on the final yields of carbon stocks, while SDI affects the rates of stand carbon sequestration. The differences in stand carbon stocks were as large as 17.72 ton ha⁻¹ between a higher (18 m) and a lower (10 m) SCI when evaluated for an average level of SDI (300 trees ha⁻¹). The gaps on the corresponding stand volume also reached about 60.41 m³ ha⁻¹. For all simulated SCIs, the carbon sequestrations for lower SDI (100 trees ha⁻¹) were less than that of the middle (300 trees ha⁻¹) and higher (400 trees ha⁻¹) SDI before the first 60 years and 80 years, respectively. After that, the situations were entirely reversed, but the differences were not remarkable. The corresponding times of the peak carbon sequestrations were 40 years of lower SDI, and 20 years of middle SDI, and 10 years of higher SDI, respectively. The carbon stocks of lower SDI were substantially less than those of middle and higher SDI values; however, the differences in the final yields between 300 and 500 trees ha⁻¹ were relatively inconspicuous.

Fig. 4

Simulation on the development of stand carbon stocks (left) and carbon sequestration (right) with stand ages for different site class index (SCI) and stand density index (SDI), where the values of SCI were 10 m, 14 m, and 18 m, and the values of SDI were 100 trees ha⁻¹, 300 trees ha⁻¹ and 500 trees ha⁻¹, respectively

Benefits of carbon sequestration afforestation

The benefits and costs of timber production and carbon sequestration for alternative combinations among different thinning scenarios, rotations, SCIs, and SDIs are shown in Tables 4 and 5, respectively. For 40-year rotations, the total benefits of larch plantations both increased dramatically with the increases of SDIs and SCIs whether the thinning treatments were implemented or not, however, the average increments of total benefits between any consecutive SDIs (about 3.04 times) were notably more extensive than that of SCIs (about 1.46 times) for scenarios with thinning treatments, whereas the differences for that of scenarios with no thinning treatments were not remarkable (1.36 vs 1.34 times). One thing to note is that the total benefits with lower SDIs were always negative or negligible, regardless of which combinations among alternative thinning scenarios, SCIs, and SDIs were considered. The marginal effects, defined as the ratios between total benefits and carbon stocks when the forest was clearcutting, also increased significantly with the increases of SDIs and SCIs. Compared with the scenarios with thinning treatments, the total benefits, carbon stocks, and stand volumes were increased by about 110.72%, 164.83%, and 88.43%, respectively, while the marginal effects were decreased by $24.36 ton⁻¹ on average.

Table 4 The total benefits and costs of carbon sequestration afforestation project for larch plantation with alternative combinations among SCIs and SDIs when the different thinning treatments were considered for 40-year rotations, where SCI, SDI, and NPV represent site class index, stand density index, and net present value, respectively

Table 5 The total benefits and costs of carbon sequestration afforestation project for larch plantation with alternative combinations among SCIs and SDIs when the different thinning treatments were considered for 60-year rotations, where SCI, SDI, and NPV represent site class index, stand density index, and net present value, respectively

The results, as mentioned earlier, could also be observed for 60-year rotations (Table 5). However, the results further indicated that the total benefits could be further increased by about 71.39% and 80.27% of scenarios with and without thinning, respectively, when the rotations were extended by 20 years (namely 60 years). Similarly, the increments on the amount of commercial timber and carbon stocks were also very substantial, 35.74% and 21.82% for timber production, and 49.15% and 21.36% for carbon stocks of scenarios with and without thinning. The average values of marginal effects for longer rotations were also significantly more extensive than those of lower rotations, which increased by 14.78% of scenarios with thinning treatments and 48.71% without thinning treatments.

For 40-year rotations with thinning treatments (Table 6), the variations of discount rates had the most significant effects on the total benefits (+ 99.09% and − 57.10% when the discount rates were 1.5% and 4.5%), followed by timber prices (± 87.37%), while the effects of carbon prices (± 2.01%) and certification costs (± 3.08%) could be almost negligible. Some differences in the amounts of variations could be observed among alternative scenarios. However, the rankings of the seven factors were never changed. No matter which kind of rotation, the effects of different factors on the total benefits of scenarios with thinning treatments were always larger than those without thinning treatments. Meanwhile, the effects of different scenarios with 40-year rotations were also larger than those with 60-year rotations, except for the effects of rotations, when the same thinning treatments were evaluated.

Table 6 Sensitivity analysis of carbon prices, discount rates, timber prices, afforestation costs, maintenance, and certification costs on the total benefits of larch plantations for different rotations in northeast China

Discussion

Afforestation is usually considered one of the least expensive solutions to address climate change, but a series of management decision problems also make it not as simple as it seems. In the development of carbon forestry, predicting and simulating the developments of carbon stocks with stand ages is a prerequisite for making management plans for forest managers. Different from the traditional methods that directly link the stand volume, BCF, and BEF together, the present study developed a compatible model system of stand volume and carbon stock for larch plantations based on the same set of data. To overcome the intrinsic correlations among different variables (e.g., HT, BAS, VOL in Eq. 8), the models were developed using the method of solving simultaneous equations. The obtained results suggest that the accuracy of the proposed models was as large as 95% on the evaluation indicator ${R}_{a}^{2}$, demonstrating the high reliability of the proposed models.

Since the system integrated both stand density and site quality variables, the models can be used in a broader stand environment. The simulations indicated that the differences on the carbon stocks and stand volume were about 17.72 ton ha⁻¹ and 60.41 m³ ha⁻¹ for different SCIs (10–18 m) and 28.74 ton ha⁻¹ and 97.97 m³ ha⁻¹ for different SDIs when evaluated for an average SCI (14 m) and SDI (300 trees ha⁻¹) levels. These results imply that site enhancement techniques such as irrigation, fertilization, and harvest residue removal strategies could increase carbon stocks [39]. Meanwhile, higher planting density is also beneficial in increasing carbon sequestration on sites of the same quality [40]. As emphasized here, the site quality variable (SCI) could represent a portion of the climate characteristics, but it did not directly capture the current warming and humidification effects in the study area, thus integrating climate variables, e.g., mean annual temperature and mean average rainfall, into stand growth models might be of value [28].

As the analysis of carbon sink afforestation is more complex than traditional afforestation, e.g., including the process of accounting, payments, and reissuance, understanding the impact of different costs on the ultimate profitability is critical. The sensitivity analysis we employed indicated that discount rates and timber prices were the two most important variables affecting the joint benefits, which have also been demonstrated by Zhou and Gao [13] and Dong et al. [14]. However, the effects of carbon prices in our analysis were not as large as expected [12, 17, 40], mainly because the base scenario of carbon prices used in our analysis was low ($5 ton⁻¹). However, this carbon price is consistent with the actual prices reflected in China's carbon trading market. The World Bank [2] also reported that about 51% of emissions covered by carbon trading or carbon are priced at less than $10 ton⁻¹ CO₂. Thus, profitability would increase, and rotation lengths could be extended, if carbon prices were increased to ranges ($40–$80 ton⁻¹ CO₂) recommended by the IPCC [41].

The calculation of total benefits we used is based on the UNFCCC rules [37], with an assumption that verification and certification of carbon sequestration are carried out every 5 years. However, the time between verification and certification may significantly affect the final benefits, and may depend on biological, market, and administration factors. For example, Juutinen et al. [42] found that an annual carbon payment mechanism might be feasible only with very high carbon prices, mainly because of the relatively high associated transaction costs. Hou et al. [17] and West et al. [16] also found that using temporary certified emission reduction accounting methods to value carbon credits for fast-growing species incentivized landowners more to participate in CDM projects than long-term certified emission reduction accounting methods. Thus, the monitoring time and administration details might also be essential decision parameters from an optimization perspective.

Managing planted forests for the joint benefits of carbon sequestration and timber production is also a complex problem. Some previous studies have indicated that the rotations of planted forests could be extended significantly when the benefits of carbon sequestrations were considered [12, 14]. This conclusion has also been confirmed by our results, where the total benefits could be increased by approximately 71.39% and 80.27% of scenarios with and without thinning. Thinning from below (thinning smaller diameter trees) is a critical technique often used for adjusting the pressures of competition at early stages. Compared with no thinning treatments, the results highlighted that commercial thinning decreased the amount of wood production, carbon stocks, and joint benefits significantly, especially for lower SDIs. The average decreases on the total benefits were as large as 131.53% of middle SDIs, but were only 32.16% of higher SDIs, which was in line with the results of Peng et al. [24]. Regardless of whether the thinning treatments were implemented or not, the average marginal effects of longer rotation lengths ($106.13 ton⁻¹) were still more extensive than that of shorter rotation lengths ($82.63 ton⁻¹), indicating the rotations of larch plantations could be further extended. However, one thing that needs to be noted is that the thinning treatments tested in this analysis may not be optimal for composite management. Peng et al. [24] and Pukkala [43] argued that plans for thinning treatments should be adjusted when the benefits of carbon are considered. Meanwhile, the impact of climate would not be negligible if the rotation lengths of planted larch forests were extended from 40 to 60 years or much longer. Lei et al. [28], for example, reported that the periodic annual increments of DBH of larch plantations were 12.23%, 10.43%, and 0.11% higher, and the mortality of trees was also 16.62%, 13.00%, and 4.17% higher under RCP2.6, RCP4.5 and RCP8.5, respectively, when compared with current climate assumptions. Since the processes of growth and accumulation may change significantly in the future, optimizing the joint benefits of timber and carbon may generate the best combinations of proposed management activities under changing climate scenarios. This hypothesis comprises our next area of investigation.

Conclusions

The present study developed a compatible model system for estimating the stand volume and carbon stock of larch plantations synchronously based on data from the NFI of China and simultaneously solved biometric equations. The fitness of the models was evaluated using fivefold cross-validation, where the statistics on the ${R}_{a}^{2}$, RMSE and MAPE of the final carbon stocks were 0.9665, 0.5635 and 3.5228 respectively. The effects of different rotation lengths, thinning treatments, site qualities, stand densities, and management costs on the joint benefits of carbon sequestration and timber outcomes were quantified. The results suggest that the total outcomes of carbon and timber from larch plantations increased significantly with an increase in site quality and stand density, regardless of which combination of rotation length and thinning treatment was assumed. Early thinning treatments decreased the joint benefits significantly by approximately 131.53% and 32.16% of middle- and higher-stand densities, however longer rotations (60 years) could enlarge the outcomes by approximately 71.39% and 80.27% of scenarios with and without thinning when compared with shorter assumed rotation ages (40 years). The sensitivity analysis indicated that the discount rates and timber prices were the two most important variables affecting the joint benefits, however the effects of carbon prices were not as large as expected under the current trading market in China. Thus, management plans with longer rotations, higher stand densities, and no thinning treatments are recommended for maximizing the joint benefits of carbon sequestration afforestation and timber production from larch plantations in northeast China.

Appendix

For a specific stand, the value of SCI can be predicted using Eq. 2, which was fixed consistently during the entire rotation. The dynamics of various stand variables (HT, DBH, BAS, NUM, VOL, and CAR) can be predicted using the following procedure. N₁ and N₂ are the stand density at time t₁ and t₂, and D₂, BAS₂, and SDI2 are the stand mean DBH, basal area, and stand density index at time t₂.

1) Predicting HT at time t₂ (Eq. 4):

$${HT}_{2}={a}_{0}{SCI}^{{a}_{1}}\left(1-\text{Exp}(-{a}_{2}{t}_{2})\right)$$

(A1)

2) Estimating N at time t₂ (Chen’s model, [20]):

$${N}_{2}={N}_{1}\text{exp}\left(-\left(0.0103+0.0003SCI\right)\left({t}_{2}-{t}_{1}\right)\right)$$

(A2)

3) Estimating DBH at time t₂ (Eq. 5):

The formulations for BAS and SDI can be defined as (Chen, 2010; West, 2015; Lei et al., 2016):

$$BA{S}_{2}=\left(\pi /40000\right){N}_{2}{D}_{2}^{2}$$

(A3)

$$SD{I}_{2}={N}_{2}\times {\left({D}_{0}/{D}_{2}\right)}^{-\beta }$$

(A4)

Thus, combining the above formulations (A3, A4) and the developed BAS model (Eq. 5), we can have:

$$\frac{\pi }{40000}{N}_{2}{D}_{2}^{2}-{b}_{0}\left[1-\text{Exp}\left(-{b}_{1}{\left(\frac{{N}_{2}\times {\left({D}_{0}/{D}_{2}\right)}^{-\beta }}{1000}\right)}^{{b}_{2}}\left({t}_{2}-{t}_{1}\right)\right)\right]=0$$

(A5)

In this equation, N₂ can be obtained from Chen’s model (Eq. A2), and the parameters (${b}_{0}$, ${b}_{1}$, ${b}_{2}$ and $\beta$) have been estimated using the dataset. Then, the Mean DBH (D₂) at time t₂ can be solved from this equation using a bisection method. This prediction procedure involves the theory behind whole-forest growth modeling, which has been widely used by Chen [20], Hong et al. [33], and Lei et al. [28].

4) Estimating SDI at time t₂ (Eq. A4):

With the predicted D₂ from step 3 and N₂ from step 2, the SDI at time t₂ (SDI₂) can be predicted using Eq. A4.

5) Estimating BAS at time t₂ (Eq. 5).

Using the new stand age t₂ and the estimated SDI₂ from step 4, then we can get the BAS at time t₂ (BAS₂) by using Eq. 5.

6) Estimating VOL at time t₂ (Eq. 6):

Using the estimated BAS₂ from step 5 and the estimated HT₂ from step 1, the VOL₂ at time t₂ can be generated using Eq. 6.

7) Estimating CAR at time t₂ (Eq. 7):

Based on the estimated VOL₂ from step 6, the CAR at time t₂ (CAR₂) can be obtained using Eq. 7.