Optimization Forest Thinning Measures for Carbon Budget in a Mixed Pine-Oak Stand of the Qingling Mountains, China: A Case Study

Abstract

Forest thinning is a silviculture treatment for sustainable forest management. It may promote growth of the remaining individuals by decreasing stand density, reducing competition, and increasing light and nutrient availability to increase carbon sequestration in the forest ecosystem. However, the action also increases carbon loss simultaneously by reducing carbon and other nutrient inputs as well as exacerbating soil CO₂ efflux. To achieve a maximum forest carbon budget, the central composite design with two independent variables (thinning intensity and thinning residual removal rate) was explored in a natural pine-oak mixed stand in the Qinling Mountains, China. The net primary productivity of living trees was estimated and soil CO₂ efflux was stimulated by the Yasso07 model. Based on two years observation, the preliminary results indicated the following. Evidently chemical compounds of the litter of the tree species affected soil CO₂ efflux stimulation. The thinning residual removal rate had a larger effect than thinning intensity on the net ecosystem productivity. When the selective thinning intensity and residual removal rate was 12.59% and 66.62% concurrently, the net ecosystem productivity reached its maximum 53.93 t·ha⁻¹·year⁻¹. The lower thinning intensity and higher thinning residual removal rated benefited the net ecosystem productivity.

1. Introduction

Forest thinning is one of the most efficient tending measures and it is also an effective management technique [1,2]. Selective thinning is one type of high thinning, which removes dominant and co-dominant trees [3]. It improves the vigor of residual trees as they benefit from the water, nutrient, and light resources no longer exploited by the felled trees [4]. Although the objectives of thinning vary, they typically include increasing the share of large diameter trees, improving the quality of timbers, increasing yield value, improving stand stability and influencing tree species composition [3] by decreasing stand density, reducing competition, and increasing the light and nutrient availability for the remaining trees primarily to promote growth of the remaining individuals [5]. Efforts to optimize carbon sequestration in forest ecosystems have mainly focused on enhancing stand biomass productivity and density by adapting thinning intensity and tree species composition [6]. However, timbers and thinning residuals (branches and foliage) are usually removed from stands after thinning which reduces carbon and other nutrient inputs [7,8,9] and the soil microclimate [10] is changed. For this reason, forest thinning may reduce carbon stocks in forest soil and vegetation [11] due to the increase of soil CO₂ efflux [7]. Thinning intensity and residual removal rate are important for carbon budget in forest systems. Litter decomposition is an important process in the global carbon cycle [12]. Its decomposition affects the soil carbon content, carbon dioxide emissions and is closely related to the chemical quality of litter types and climatic conditions [13] in forest ecosystems, and in particular for litter.

A considerable amount of literature has addressed the effects of forest thinning on the forest carbon cycle. Based on these observations, some studies examined the effects of forest thinning on soil CO₂ efflux with equivocal results ranging from positive effect on soil carbon release [14], negative effect [5,15] and no effect [16,17,18]. Others reported that intensive biomass harvesting may negatively impact carbon stocks in forest soil and vegetation [11], forest thinning did not have a significant impact on carbon stocks or fluxes [19], and remaining residues after harvesting increased carbon storage [20]. Unfortunately, few studies have examined the concurrent effects of two or more factors simultaneously, e.g., thinning and residual removal (branches, needle/leaf and twigs removed with thinning) on the forest net ecosystem productivity (NEP). In addition, how to balance the tree biomass increase and carbon loss as a consequence of forest thinning is still uncertain.

To attain a preliminary combination of thinning intensity and thinning residual rate for the largest carbon budget in thinned stands, a central composite design with two independent variables (thinning intensity and thinning residual removal rate) was explored in a natural pine-oak mixed stand in the Qinling Mountains, China.

The objectives of our study were (i) to examine the effect of a chemical compound of litter on the stimulation of soil CO₂ efflux; (ii) to estimate the effects of selective thinning and residual removal on the carbon budget of the pine-oak mixed stand; and (iii) to optimize the combined thinning intensity and residual removal rate to achieve the highest total carbon budget.

2. Material and Methods

2.1. Site Description

The Qinling Mountains (32°30′–34°45′ N, 104°30′–112°45′ E) in central China constitute a substantial physical obstacle for northward and southward movement of air masses, because of their high elevation and east-to-west arrangement. Therefore, these mountains are critical in stabilizing the distribution of climate and life zones in eastern China [21]. The pine-oak mixed stand is extensive in the middle of the altitudinal gradient of the Qinling Mountains and it plays important ecological roles in water purification and carbon sequestration.

Experiments were conducted at the Qinling National Forest Ecosystem Research Station (QNFERS), located on the southern slope of the Qinling Mountains, in Huoditang, Ningshan County, Shaanxi Province (32°18′ N, 108°20′ E). The altitude of the study area is from 1500–2500 m above sea level. The area experiences a subtropical climate, with annual mean temperatures around 8–10 °C, annual mean precipitation around 900–1200 mm, and annual mean evaporation around 800–950 mm. The main soil type is mountain brown soil, developed from granite material, with depths ranging from 30 to 50 cm. The total forest area is 2037 hectares in the station. Natural forest occupies 93% of the total forest area in QNFERS, with various vegetation types distributed in this region along the altitudinal gradient, such as evergreen deciduous mixed forest (pine-oak mixed forest), deciduous broad-leaved forest (oak, red birch), temperate coniferous forest (Chinese red pine, Armand pine) and cold temperate coniferous forest (spruce, fir). The most dominant forest type is pine-oak mixed forest with an average stand age of the stands of 42 years and an average height of 9.2 m. Common tree species include Pinus tabulaeformis, Pinus armandii, and Quercusaliena var. acuteserrata. The understory species are abundant.

2.2. Experimental Design

Our experimental plots were located on steep slopes (average slope gradient 30°) with thin soil depth (<50 cm), and in fragmented terrain. These characteristics make it extremely difficult to obtain the amount of replications for a randomized block design or orthogonal experiment design.

The Central Composite design (CCD) is the most popular of the many response surface methodology (RSM) classes, and is widely used for estimating second-order response surfaces [22]. The application of these statistical techniques in experiments has the advantages of requiring fewer resources (time, numbers of duplication, and amount of experimentation), but can also reduce process variability [23]. RSM is a collection of statistical and mathematical techniques useful for developing, improving, and optimizing processes [24]. These techniques relate a response variable to predictors that have multiple levels. The coded levels and the actual levels of the independent variables were calculated according to Equations (1)–(3):

$$X_{0j}=\frac{X_{1j}+X_{2j}}{2}$$

(1)

$${\mathrm{\Delta }}_j=\frac{X_{2j}-X_{0j}}{\alpha }$$

(2)

$$x_{\mathsf{\alpha }j}=\frac{X_{\mathsf{\alpha }j}-X_{0j}}{{\mathrm{\Delta }}_j}$$

(3)

where X_j is the real value of the independent variable, X₀ is the real value of an independent variable at the center point, ∆_j is the step change value, and x_j is the coded value of an independent variable [25].

The CCD consists of 2^k full or 2^k−1 half-replicate (k = number of independent variables) factorial points (±1, ±1, …, ±1); 2k axial or star points of the form (±α, 0, …, 0), (0, ±α, …, 0); and a center point (0, 0, …, 0). The axial points are replicated one and two times, and allow for the efficient estimation of pure quadratic terms. The center points are replicated one and three times, and provide information about the existence of curvature. The number of center runs can be altered, providing flexibility to improve error estimates and power. Finally, the factorial points allow estimation of the first-order and interaction terms. The CCD can be summarized with the following equation (Equation (4)):

$$N=f+\left(2k\right)\mathsf{\alpha }+n_0$$

(4)

where N is the total number of experimental runs, f is the number of factorial points, 2k is the axial point, α is the number of times the axial point is replicated, and n₀ is the center point. The axial distance, α is chosen based on the region of interest. Selecting the appropriate values of α specifies the CCD type, with $\mathsf{\alpha }=\sqrt{k}$ being a spherical CCD [22].

The relationship between response and predictor levels can be approximated with a second-order response surface mode [22] (see Equation (5)):

$$y={\mathsf{\beta }}_0+\sum _{i=1}^k{\mathsf{\beta }}_i{x}_i+\sum _{i=1}^k{\mathsf{\beta }}_{ii}{x}_{ii}^2\sum _{j=i+1}^k\sum _{i=1}^{k-1}{\mathsf{\beta }}_{ij}{x}_i{x}_j+{\mathsf{\epsilon }}_{ij}$$

(5)

where y is the measured response; β₀, β_i, β_j, β_ii, and β_ij are parameter coefficients; x_i, x_j are the input variables; and ε_ij is an error term.

Analysis of Variance (ANOVA) is used to optimize Equation (5) and analyze the interaction effect of the input variables on measured response and the single effect of input variables on measured response by differentiating variable j on variable i and vice versa [25].

The contributions of the controlled variables to the dependent variable, as F > 1 were estimated following the method of Tang [25] (see Equations (6) and (7)).

$${\mathrm{\Delta }}_j=S_j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\sum S_{ij}+S_{jj}$$

(6)

$$S=1-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$F$}\right.$$

(7)

where Δ_j is the contribution of controlled variable j to the dependent variable; S_j is the linear term for the controlled variable j; S_ij is the interaction term for the controlled variables i and j; S_jj is the quadratic term for the controlled variable j; F is the F-value in the ANOVA.

Our experiment was based on CCD, generated using Data Processing System (DPS) version 14.50 [25].

In a preliminary investigation conducted over the 15–25 August 2012, we selected 13 plots (20 m × 20 m) with similar slope gradients, canopy cover, tree species composition, and soil depths (

Table 1. General information of plots. Composition of tree species was calculated by their basal area.

Plot No.	Gradient (°)	Canopy Density		Tree Species Composition			Soil Depth (cm)
		2012	2013	2012	2013	2014
1	35	0.75	0.72	7 Pt 2 Pa 1 Q	9 Pa 1 Pt	9 Pa 1 Pt	30
2	30	0.8	0.73	5 Pt 3 Pa 2 Q	5 Pa 1 Pt 4 Q	5 Pa 1 Pt 4 Q	45
3	30	0.8	0.72	7 Pt 2 Pa 1 Q	7 Pa 2 Pt 1 Q	7 Pa 2 Pt 1 Q	40
4	30	0.8	0.73	6 Pt 2 Pa 2 Q	2 Pa 4 Pt 4 Q	2 Pa 4 Pt 4 Q	50
5	25	0.85	0.81	4 Pt 3 Pa 3 Q	1 Pa 5 Pt 4 Q	1 Pa 5 Pt 4Q	46
6	30	0.8	0.7	1 Pt 7 Pa 2 Q	1 Pa 8 Pt 1 Q	1 Pa 8 Pt 1 Q	35
7	30	0.75	0.73	2 Pt 7 Pa 1 Q	8 Pt 2 Q	8 Pt 2 Q	40
8	30	0.9	0.81	2 Pt 5 Pa 3 Q	2 Pa 5 Pt 3 Q	2 Pa 5 Pt 3 Q	44
9	25	0.9	0.84	2 Pt 7 Pa 1 Q	1 Pa 7 Pt 2 Q	1 Pa 7 Pt 2 Q	33
10	30	0.8	0.72	1 Pa 8 Pt 1 Q	5 Pa 1 Pt 4 Q	6 Pa 4 Q	38
11	30	0.9	0.83	6 Pt 2 Pa 2 Q	1 Pa 3 Pt 6 Q	1 Pa 3 Pt 6 Q	41
12	25	0.8	0.72	6 Pt 3 Pa 1 Q	7 Pa 2 Pt 1 Q	8 Pa 2 Q	45
13	28	0.8	0.73	8 Pt 1 Pa 1 Q	9 Pa 1 Q	8 Pa 1 Pt 1 Q	45

Pa, Pt, and Q in the table is Pinus armandi, Pinus tabulaeformis, and Quercus aliena var. acutesserata respectively.

) for the current experiment. For each plot, we surveyed the soil depth, as well as the height (m), diameter at breast height (DBH, cm), and canopy cover (%) of the tree species. The intensity of selective thinning (ST, %) and thinning residual removal rate (TRR, %) were calculated using Equations (8) and (9), respectively:

$$ST=\frac{A_f}{A_T}$$

(8)

$$TRR=\frac{Q_i}{Q_t}$$

(9)

where A_f, A_T, Q_i, and Q_t are the basal area of logged trees, total basal area of trees, fresh weight of removed residue, and fresh weight of total residue in all plots.

The design consisted of two independent variables (X₁ = thinning and X₂ = thinning residual removal), each with five intensity/rate gradients (

Table 2. Experiment design runs in DPS v14.50 (Data Processing System).

Variables	Levels
	(−α) − 1.414	−1	0	1	(+α) + 1.414
X1 (%)	5.00	7.93	15.00	22.07	25.00
X2 (%)	0.00	14.64	50.00	85.36	100.00

X₁ and X₂ in the table represents actual thinning intensity and thinning residual removal rate respectively. Levels (−α, −1, 0, +1, and +α) in the table are coded values of variables (thinning intensity and thinning residual removal rate) generated by the software Data Processing System (version 14.50).

). For the controlled factor (independent variable) in the current study, the value α was $\sqrt{2}$ = 1.414. We set the +α and −α level thinning intensity to 25% and 5% respectively according to the Regulation for Tending of Forest [26], and 100% and 0% for the residual removal rate. To explore the effects of the thinning operation on NEP, zero-treatment of thinning intensity was excluded.

For a central composite design with two independent, five-level variables, 13 experimental runs are required, with four factorial points from treatment I to treatment IV, four axial points from treatment to treatment VIII, and five center points treatment IX (

Table 3. Experiment design runs in DPS v14.50.

Plot No.	Treatment	Design Code		Thinning Factors		Block
		x1	x2	X1 (%)	X2 (%)
1	I	1	1	22.07	85.36	Factorial points
2	II	1	−1	22.07	14.64
3	III	−1	1	7.93	85.36
4	IV	−1	−1	7.93	14.64
5	V	−α	0	5	50	Axial points
6	VI	+α	0	25	50
7	VII	0	−α	15	0
8	VIII	0	+α	15	100
9	IX	0	0	15	50	Center points
10	IX	0	0	15	50
11	IX	0	0	15	50
12	IX	0	0	15	50
13	IX	0	0	15	50

). The factorial points were a combination of controlled variables at ±1 levels (a thinning intensity at ±1 level represent 22.07% and 7.93% respectively; a thinning residual removal rate at ±1 levels represent 85.36% and 14.64% respectively) in our study. Similarly, the star points were a combination of controlled variables at ±α and 0 levels (a thinning intensity at ±α and 0 levels represented 25%, 5%, and 15% respectively; a thinning residual removal rate at ±α and 0 levels represented 100%, 0%, and 50% respectively). The center point was a combination of controlled variables at 0 levels. Different thinning factors were applied in each plot, except for the plots categorized as center points (Table 3).

The dependent variable in this study was the average of NEP in 2013 and in 2014 in post-treatments. The experimental results were fitted to a second-order polynomial model, and the regression coefficients were determined. The quadratic model for predicting the optimal combination of thinning intensity and removal rate to reach the highest value of NEP (Y_k) is described by Equation (10):

$$Y_k=b_{k0}+\sum _{i=1}^2{b}_{ki}{x}_i+\sum _{i=1}^2{b}_{kii}{x}_i^2+\sum _{i=1}^2{b}_{kij}{x}_i{x}_j$$

(10)

where b_k0, b_ki, b_kii, and b_kij are the constant regression coefficients of the model, and x_i, x_j are codes of the independent variables (x_i = thinning intensity and x_j = thinning residual removal rate).

2.3. Thinning, Residual Removal, and Dynamics of Tree Growth and Litterfall

All trees in the 13 plots with DBH > 10 cm were numbered and tagged between 28 and 30 August 2012. Crop trees in plots were selected. Stem quality, crown size, vitality, spatial distribution of potential crop trees, diameter, and tree damages were taken into account when selecting the crop trees. Competing trees were marked and cut. The intensity of thinning in each plot was determined as shown in Table 2. After log harvesting, the total fresh weight of residuals (branches, needle/leaf and twigs) in each plot was measured. Then, the branches and twigs were cut into pieces of length 60–80 cm and mixed with needle/leaf. According to the experiment design (Table 2), a part of the residuals was removed and the rest was thrown on the forest ground with a similar thickness in each plot. The actions of thinning and residual removal were done manually for steep slopes of plots in September 2012.

The height and DBH of the remaining trees were monitored between 20 and 28 September 2013 and 2014.

Along the slope from bottom to top, each plot was mechanically partitioned in three sections. Nine circular litter traps (diameter 30 cm) were set equidistantly in each plot to monitor the dynamics of litterfall since 20 September 2012. At the end of December 2013 and 2014, littler in each plot was mixed separately in the field and then was brought to the laboratory to dry and measure dry weight.

2.4. Soil CO₂ efflux

For the high spatial variability of the soil carbon stocks and the high uncertainty in their changes [27], it is difficult, laborious and expensive to measure soil CO₂ efflux directly [13]. Models are needed to estimate the dynamics of carbon in forest soils [28]. Comparing existing soil carbon models Century [29], Q-model [30], ROMUL [31], RothC [32] and DECOMP [33] to Yasso07, the significant advantage of Yasso07 is that the parameters to operate the model are easily accessible [34]. Therefore the Yasso07 model was applied in the current study to estimate CO₂ efflux in post-treatments.

2.5. Chemical Analysis

A mixed litter sample 0.50 kg was collected from nine litter traps in each plot. Samples were separated by conjunction of litter types (twig, needle, leaf) and tree species. Chemical compound groups, ESC (ethanol soluble compound), WSC (water soluble compound), ASC (acid soluble compound) and NSC (non-soluble compound) of litter [28] were analyzed by Bai et al. [35]. The analysis processing was the following. (1) ESC: Dried and ground litter sample (diameter is 0.074 mm) of 1.00 gram was put into a cylinder made by filter paper with diameter 4 cm. Then, the cylinder with the sample was removed to a reflux line of a Soxhlet extractor (Yuming Instrument Co. Ltd., Shanghai, China) (BSXT-02-250). Next, 100 mL mixed benzene-alcohol (1:1) solution was added into a fat-wax bottle connected with a reflux line and a condenser pipe and the bottle was placed in a 80 °C thermostat water bath for 20 h until the color of the solution in the reflux line disappeared completely. Afterwards, the cylinder was taken out and put into a draught cupboard until the benzene-alcohol solution volatilized entirely and the first residue was attained. Finally, the residue was dried in a 50 °C oven for 4 h until its weight remained unchanged. The weight of ESC was calculated from the difference between the weight of the sample and the first dried residue; (2) WSC: The first dried residue was moved to a 250 mL beaker and 150 mL distilled water was injected in, and the residue was stirred and broken into pieces with a glass rod. Then, the beaker was covered and placed in a 100 h hydrolysis pot for 3 h. Next, a sand core funnel (type G3) was used to suck filter the contents and the second residue was attained. Afterwards, the second residue was washed in 30 °C pure water until the eluate was without color. After that, the colatuie and the eluate were injected into a new 250 mL beaker and its volume maintained at 200 mL. Finally, the second residue was removed from the sand core funnel to another beaker and dried in a 50 °C oven for 4 h. The weight of WSC was the difference between the weight of dried residues obtained from the first and the second time; (3) ASC: The 150 mL and 2% hydrochloride resolution was injected in the beaker filled with the second dried residue. Then, the residue was stirred and broken into pieces with a glass rod. Next, the beaker was put into a 100 h hydrolysis pot for 5 h and the solution was suck filtered with a G3 sand core funnel and the third residue was attained. After this, the third residue was placed in a porcelain crucible and was washed in 30 °C pure water until sulfate radical could not be detected by 5% barium chloride solution. Finally, the third residue was dried in a 105 °C oven for 4 h. The weight of ASC was the difference between the weight of the dried residues obtained from the second and third time; (4) NSC: The third dried residue was placed in a porcelain crucible and the crucible was put on a 45 °C electric stove to carbonize the third residue for 3 h. Then, the crucible was placed on a 450 °C electric stove to burn the third residue for 8 h. Later, the residue was cooled naturally until room temperature. Finally, the cooled residue was weighed to obtain the fourth residue. The weight of NSC was the difference between the weight of the third dried residue and the fourth.

All these parameters were inputs to run the Yasso07 model.

2.6. Data Processing and Analysis

Ratio of tree species (R_i) was calculated as Equation (11).

$$R_i=\frac{b_i}{B}\times 100\%$$

(11)

where b_i and B is basal area of the tree species i and total basal area of tree species in an identical treatment.

Composition of tree species is the proportion R_i: R_m: R_n.

Where R_i, R_m and R_n is ratio of tree species i, m and n respectively.

Diameter classes were used to describe the DBH dynamics of tree species (

Table 4. Classification of diameter classes.

Diameter Class (cm)	DBH (cm)
4	≤6
8	6.1 ≤ DBH ≤ 10
12	10.1 ≤ DBH ≤ 14
16	14.1 ≤ DBH ≤ 18
20	18.1 ≤ DBH ≤ 22
24	22.1 ≤ DBH ≤ 26
28	26.1 ≤ DBH ≤ 30
32	30.1 ≤ DBH ≤ 34
36	34.1 ≤ DBH ≤ 38
40	38.1 ≤ DBH ≤ 42
44	42.1 ≤ DBH ≤ 46
48	46.1 ≤ DBH ≤ 50

). DBH ratio of tree species (D_ij) was calculated as follows:

$$D_{ij}=\frac{n_{ij}}{N_i}\times 100\%$$

(12)

where n_ij is number of tree species i in diameter calss j and N_i is number of total tree species in all diameter classes of 13 plots.

Chemical compound groups (ESC, WSC, ASC, and NSC) of litterfall were calculated by mass weighted average of tree species. We inputted the quality of litterfall from each post-treatment, the measured chemical compound groups of tree species, and the data of annual precipitation and air temperature from QNFERS into the Yasso07 model to estimate soil CO₂ efflux (Rs) [28].

Based on monitoring DBH and the height of remaining trees in the thirteen plots in 2013 and 2014, living biomass (Mg·ha⁻¹) of whole remaining trees was estimated as Equations (13)–(15) [36].

P. tabulaeformis:

$$Y=15.525+0.6269x$$

(13)

P. armandi:

$$Y=54.280+0.4048x$$

(14)

Q. aliena var. acutesserata:

$$Y=13.394+1.0564x$$

(15)

where Y is the living biomass of trees and x is the stand growing stock (m³·ha⁻¹).

The stand growing stock was calculated by the stems and volume of each tree. The volume of a single tree was calculated as follows [37] (see Equations (16)–(18)):

P. armandi:

$$\mathrm{lnV}=0.95697\ln \left({\mathrm{D}}^2\mathrm{H}\right)-9.95738$$

(16)

P. tabulaeformis:

$$\mathrm{lnV}=0.99138\ln \left({\mathrm{D}}^2\mathrm{H}\right)-10.20211$$

(17)

Q. aliena var. acutessera:

$$\mathrm{lnV}=0.96884\ln \left({\mathrm{D}}^2\mathrm{H}-10.07352\right)$$

(18)

where V (m³), D (cm), and H (m) are volume, DBH, and height of a tree respectively.

Net primary productivity (NPP) of the current forest was the living tree biomass increment in two consecutive years multiplied by the carbon ration in plants (0.50 in this study).

NEP was calculated as following.

$$\mathrm{NEP}=\mathrm{NPP}-{\mathrm{R}}_{\mathrm{s}}$$

(19)

Figures were plotted using Origin8.0 (OriginLab Corporation, Northampton, MA, USA) software. DPS v14.50 software (Zhejiang University, Hangzhou, China) was used to fit models, analyze the data, determine the effects of a single independent variable and the interaction of independent variables on NEP and optimize the combination of thinning intensity and residual removal rate for the highest NEP.

3. Results

3.1. Dynamics of Tree Species Composition

Selective thinning decreased the canopy density and changed the composition of tree species (Table 1). The proportion of Quercus aliena var. acuteserrata of the total tree species in each plot was not more than 30% before thinning in 2012 and increased in most plots after thinning in 2013 and 2014 (Table 1).

3.2. DBH Dynamics of Tree Species

Selective thinning and thinning residual removal promoted DBH of the remaining tree species increase in post-treatments (Figure 1).

3.3. Dynamics of Net Primary Productivity of Living Trees

Based on height and DBH of the tree species monitored, the net primary productivity of living trees was estimated. The net primary productivity of living trees decreased between pre-treatments and post-treatments (Figure 2). Comparing to the start point in 2012, the average increment of biomass carbon of living trees decreased respectively 20.46 t·ha⁻¹·year⁻¹·and·15. 64 t·ha⁻¹·year⁻¹·after the first year thinning in 2013 and after the second year thinning in 2014 (Figure 2). With tree growing after thinning, the net primary productivity of the living trees gradually increased after thinning for two years in 2014 (Figure 2).

3.4. Soil CO₂ Efflux Stimulation

Based on the chemical compound groups of the litter of the tree species we measured and data of annual precipitation and air temperature in 2013 and 2014 from QNFERS, soil CO₂ efflux was estimated. The results stimulated by Yasso07 demonstrated that selective thinning and residual removal had a hysteretic effect on soil CO₂ efflux. Soil CO₂ efflux in post-treatments in 2013 was lower than that in 2014, except for treatment III (Figure 3).

3.5. Model Fitting

Effects of selective thinning intensity and thinning residual removal rate on net ecosystem productivity were analyzed by fitting a quadratic model. The final quadratic model was obtained for each response and was expressed by the following second-order polynomial equation after optimization.

$$\mathrm{NEP}=50.10-1.87x_1+2.85x_2-2.00x_1^2-2.54x_2^2+2.48x_1{x}_2$$

(20)

Regression coefficient, standard error, and ANOVA for the regression model of NEP are presented (

Table 5. Regression coefficient, standard error, and ANOVA (Analysis of Variance) for the regression model of net ecosystem productivity (NEP).

Parameter	Degrees of Freedom	Sum of Squares	Mean Square	F-Value	Pr > F
x1	1	27.9913	27.9913	5.5179	0.0512
x2	1	65.1148	65.1148	12.8360	0.0089
x12	1	27.9793	27.9793	5.5155	0.0512
x22	1	44.9868	44.9868	8.8682	0.0206
x1x2	1	24.5520	24.5520	4.8399	0.0637
Model	5	182.4717	36.4943	7.1941	0.0226
Lack of fit	3	13.3448	4.4483	0.8028	0.5308
Residual	7	35.5097	5.0728
Error	4	22.1649	5.5412
Total	12	217.9815
R2 = 0.837
Adjusted R2 = 0.849

). Under the condition p < 0.10, the quadratic model was a better fit relationship of NEP in conjunction with the selective thinning rate and thinning residual removal rate.

3.6. Effects of Forest Management on Carbon Budget

The model Equation (16) demonstrated that both variable x₁ and x₂ affected NEP and they also had an interaction effect on NEP. According to Equations (6) and (7), and the data shown in Table 5, the thinning removal rate (Δx₂ = 2.21) had a larger effect than selective thinning intensity (Δx₁ = 2.03) on NEP.

Effects of single factor and their interaction on NEP were analyzed by software DPS v14.50 respectively. The results indicated that selective thinning intensity and thinning residual removal rate were positively related to NEP when the range of the independent variable was x₁∈ [−1.414, −0.5] and x₂∈ [−1.414, 0.5] respectively (Figure 4). In contrast, as independent variables ranged in x₁∈ [−0.5, 1.414] and x₂∈ [0.5, 1.414], both selective thinning intensity and thinning residual removal rate were negatively related to NEP (Figure 4). The modeled effects of thinning and residual removal on NEP are illustrated in a 3D-contour plot (Figure 5). The effect of increasing residual removal rate on NEP was conspicuous when selective thinning intensity was at a low level x₁∈ [−1.414, −0.25] (Figure 5). Thereafter, with an increase of selective thinning intensity, NEP showed a slow increase with increased residual removal rate. Thinning intensity and residual removal rate thus interacted negatively in their effects on NEP. The outcomes suggested that lower selective thinning intensity and higher thinning residual removal rate benefited NEP.

3.7. Optimization Forest Management Measures for NEP

Each independent variable (thinning intensity, residual removal rate) was individually increased or decreased in an attempt to find the maximum response in NEP. Once the optimal value was found separately for thinning intensity and residual removal rate, the value was selected as the condition for obtaining the overall maximum NEP. Our analysis demonstrated that the maximum NEP (53.93 t·ha⁻¹·year⁻¹) was achieved when the independent variables were x₁ = −0.17 and x₂ = 0.48. To verify these index values, the codes for the independent variables were incorporated into Equations (5)–(7). For NEP, these codes yielded selective thinning intensity and residual removal rate of 12.59% and 66.62%, resulting in the predicted maximum NEP of 53.93 t·ha⁻¹·year⁻¹.

4. Discussion

4.1. Effect of Chemical Compound Groups of Litterfall on Soil CO₂ Efflux

Chemical compound groups of litter on Euro-American tree species are provided in the Yasso07 manual which are more convenient for the model users in those countries [12]. Whether the parameters of tree species in China affect soil CO₂ efflux stimulation is uncertain. We analyzed chemical compound groups of litter types (leaf/needle, fine root, twig, and coarse root) of three tree species (

Table 6. Chemical compound groups of litter among tree species.

Tree Species	Ethanol	Ethanol Std	Water	Water Std	Acid	Acid Std	Non Soluble	Non Soluble Std	Litter Types
Pt	0.125a	0.004	0.153b	0.006	0.433w	0.003	0.289c	0.006	Leaf/needle
Pa	0.122a	0.005	0.151b	0.006	0.418w	0.005	0.309c	0.011	Leaf/needle
Q	0.199b	0.004	0.252c	0.004	0.445w	0.005	0.104d	0.002	Leaf/needle
Pt	0.120a	0.005	0.142d	0.002	0.477w	0.003	0.261e	0.007	Fine root
Pa	0.119a	0.002	0.130ad	0.002	0.476w	0.002	0.276e	0.003	Fine root
Q	0.154c	0.003	0.187e	0.009	0.429w	0.009	0.229f	0.015	Fine root
Pt	0.086k	0.004	0.088kf	0.004	0.525m	0.027	0.301c	0.005	Twig
Pa	0.061e	0.006	0.097f	0.005	0.523m	0.005	0.321c	0.005	Twig
Q	0.038u	0.002	0.081f	0.008	0.562m	0.005	0.319c	0.32	Twig
Pt	0.073g	0.002	0.091f	0.003	0.530m	0.002	0.306c	0.004	Coarse root
Pa	0.072g	0.004	0.085f	0.004	0.524m	0.004	0.319c	0.006	Coarse root
Q	0.111h	0.003	0.139g	0.013	0.473w	0.007	0.278f	0.007	Coarse root

Ethanol, water, acid and non-soluble in the table represents content of ESC, WSC, ASC, and NSC in the litter respectively. The acronym std is standard error. Different lowercase letters indicate differences among tree species and litter types (leaf/needle, fine root, twig, and coarse root) at p < 0.05 level by LSD. The acronym of tree species in Table 6 is as same as in Table 1.

). The results indicated that content of ESC, WSC, ASC, and NSC among different tree species with the identical litter type varied significantly (Table 6). A similar trend was also found among litter types with the same tree species (Table 6). To examine the effect of chemical compound groups of litterfall on soil CO₂ efflux, the values we measured (scenario 1) and the global from Yasso07 manual (scenario 2) [13] with litterfall quality, precipitation and air temperature in 2013 in each treatment were adopted in running the Yasso07. The stimulated soil respiration was conspicuously underestimated in scenario 2 than that in scenario 1 with identical treatments (Figure 6). The average underestimated soil CO₂ was 2.55 t·ha⁻¹·year⁻¹ and the maximum even reached 3.50 t·ha⁻¹·year⁻¹ (treatment IX) (Figure 6).

Biological characteristics of tree species might lead to varieties of chemical compound groups of litter [38]. Soil CO₂ efflux results from chemical compound groups of litter transforming into different soil carbon compartments [13]. Chemical compound groups of litter differed from each other between tree species even in the same genera. ESC in needles of Pinus sylvestris was 10.79 times of that in Pinus pinaster, 8.56 times of that in Pinus pinea and 7.24 times of that in Pinus banksiana respectively [12]. In addition, NSC in leaf of Quercus robur was 5.41 times of thath in Quercus garryana and ASC in leaf of Quercus garryana was 1.54 times of that in Quercus robur [38,39]. Cellulose in litter is usually shielded by lignin and the litter decomposition rate decreases with its NSC for high nitrogen and lignin content [40]. ASC is composed of cellulose and hemicellulose, and its decomposition requires cellulase and other special conditions [12]. High NSC and ASC in litter may decrease the litter decomposition rate [12]. In contrast, high ESC and WSC in litter will promote litter decomposition resulting in the main components of dissolved fat, pigment and oil in ESC, and sugars in WSC, with low nitrogen [12].

The current study suggested that the chemical groups of litterfall apparently affected the result of soil CO₂ efflux. We recommended that the chemical compound groups of the litter should be measured before applying the Yasso07 model to stimulate soil CO₂ efflux.

4.2. The Response of Soil CO₂ Efflux to Management Measures

We observed that selective thinning and thinning residual removal increased soil CO₂ efflux with a hysteretic effect. Soil CO₂ efflux in 2014 was generally higher than that in 2013 with identical treatments (Figure 3). Most studies have indicated that temperature and moisture are the main factors positively influencing soil respiration over various climate regions [41,42,43]. With timber harvesting, thinning residual removal, canopy openness, and residual decomposition, more light with rain having reached the forest ground, increases soil temperature and moisture, and activates soil microorganisms. Thinning also increased soil respiration [14]. Others reported that the prescribed thinning had a negligible effect on soil respiration [16,17,18]. The likely reasons were that thinning induced increases in shrub abundance might be responsible for commensurate increases in fine root turnover which contributed substantially to the increased light use efficiency of thinned plots [16]. Sullivan et al. indicated that thinning may reduce soil respiration by killing trees, by altering the soil environment, or by changing the amounts and sources of below ground carbon for microbial metabolism [15]. Based on a short-term data, this study reported a preliminary effect of forest thinning on soil CO₂ efflux. Future monitoring will assist in clarifying the relationship between soil CO₂ efflux and forest thinning.

5. Conclusions

Our two years observation preliminary demonstrated that the chemical compound groups of the litter of tree species should not be ignored in analyzing soil CO₂ efflux. Thinning intensity and thinning residual removal rate have different effects on NEP. The thinning removal rate had a larger effect than selective thinning intensity on NEP. When selective thinning intensity and residual removal rate was 12.59% and 66.62% concurrently, the NEP reached its maximum 53.93 t·ha⁻¹·year⁻¹. A lower thinning intensity and higher thinning residual removal rate benefited NEP.