Exploring the Effects of Thinning on Cunninghamia lanceolata Lamb. Carbon Allocation in Southwestern China Using a Process-Based Model

Abstract

We investigated the effects of thinning intensity on the carbon allocation of Cunninghamia lanceolata Lamb. Hook by analyzing the stand growth and carbon content of a plantation under three thinning intensities (I: 70%; II: 50%; III: 30%) and with no thinning treatment. Using the carbon balance framework of the CROwn BASe (CROBAS) model and multi-source inventory data, we calibrated the parameters of the CROBAS-C. lanceolata (CROBAS-CL) model to simulate the carbon content in the plantation. We validated the CROBAS-CL model by comparing the predicted stand diameter at breast height (DBH) and stand height (H) with the measured values. Finally, the predicted stand carbon was compared with the soil carbon to assess the dynamics and allocation of ecosystem carbon content. Overall, our findings suggest that the predicted stand carbon of CROBAS-CL satisfies the statistical test requirements: the deviation of height and DBH predicted by the CROBAS-CL model from the measured height and DBH are less than 0.087 m and 0.165 cm, respectively. These results confirm that the model is useful for a dynamic prediction of stand carbon in C. lanceolata plantations. Based on the results of the proposed model, we determine that Thinning III (30% thinning intensity) is beneficial for the growth of C. lanceolata plantations and improving soil carbon sequestration. Additionally, the simulated carbon storage of an individual tree in the C. lanceolata plantation gradually increased with the tree age. Our study provides a strong reference for the efficient operation and management of C. lanceolata plantations in southwestern China.

1. Introduction

The rapid transformation of our planet due to global climate change is an indisputable truth. Afforestation and reforestation in the largest carbon pool of terrestrial ecosystems—forests have made tremendous contributions to mitigating global climate change [1]. Specifically, the sustainability of forest management is being used to improve the carbon sink function of forest ecosystems. Forests play a pivotal role in the potential effects of climate change on terrestrial carbon sequestration and carbon stocks. Approximately 52% of global forests are managed with varying intensity [2]. Thus, understanding if, how, and to what extent different forest management practices may modify the processes that control carbon dynamics during undisturbed stand development and in response to climate change is crucial in improving our understanding of land-based climate mitigation capacity [3,4].

It is widely believed that the biomass allocation of forest ecosystems is closely related to carbon storage. The biomass in forest ecosystems constitutes the largest carbon pool and accounts for 90% of the entire terrestrial ecosystem [5]. The vegetation biomass is often increased by afforestation to improve the carbon accumulation of the terrestrial ecosystems. However, several unknown factors can affect forest management as well as the carbon cycle and carbon sink function of the plantation ecosystems, suggesting an urgent need for further investigation [6]. Afforestation plays an important role in accumulating soil carbons and promoting the combination of organic carbons and micro-aggregates to enhance the stability of the carbon. This is one of the reasons for higher soil carbon storages in plantations than in farmland and grassland ecosystems [7]. Moreover, due to the slow turnover of the soil carbon pool in afforested areas, carbon storages can be maintained for an extended period. Fang et al. [8] suggested that the carbon sink of Chinese forests has mainly originated from plantations over the last 20 years. Considering the significant role plantations in reducing CO₂ discharge, they are being widely investigated globally to reduce greenhouse gas emissions. In recent years, several studies [9,10] have focused on the carbon storage of afforested ecosystems using different tree species, forest age and densities, and their allocation patterns, suggesting that carbon storage in plantations increases with the age of the forests [11,12]. However, the conclusions obtained by studying the impact of stand densities on forest carbon storage are contradictory and primarily based on natural forests. Wei and Blanco [12] suggested that carbon storages in plantations in the subtropics of southwestern China might increase with the age of forests. Conversely, Scott et al. [13] studied Pinus radiata D. Don plantations in New Zealand and concluded that carbon storages in afforestation might decrease with the age of the forests. Although these studies have enhanced the understanding of carbon sinks and sequestration in forests, the correlation between carbon storages of stand and other parts (soil, understory plants, and litter) are still unclear.

Computer simulation models have emerged as effective research tools for the quantitative measurement of ecosystem productivity at global and regional scales. Several ecosystem models have been reported globally. Based on the principle of simulations, these models can be categorized into two groups: empirical models and process-based models [14]. The empirical models are primarily based on statistical methods and utilize climatic factors including temperature and precipitation to estimate the primary productivity of plants. These methods are relatively accurate in a specific range. Unfortunately, they exhibit major limitations in terms of applicable conditions and range. Furthermore, they are not sufficient to reveal the mechanism of energy flow and material circulation [15]. In contrast, process-based forest growth models are valuable tools to evaluate forest dynamics, development, management, and changing climate and to assess the long-term effects on forest carbon cycling [16]. Compared to the empirical models, process-based models have complicated variables, numerous parameters, and complex forms. Furthermore, the number of parameters they employ is gradually increasing with the constant evolution of the models [17]. There is, in fact, a long-standing interest and a pressing need for including a detailed representation of forest management in coupled land-climate models for scenario analyses. Under the circumstance of insufficient validation data, it is always beneficial to use the process-based models to predict stand growth [15].

The preliminary process-based model, CROBAS, is a carbon balance model that was developed for Pinus sylvestris L. by Mäkelä [18], which uses functional relationships and carbon balance modeling to describe tree growth and biomass allocation. This model, which contains 39 parameters, is based on the principles of light entrapment, photosynthesis, and carbon allocations. Several studies have successfully validated the reliability of this model for simulating or predicting stand growth. In particular, its accurate estimation of carbon storages in individual trees or even-aged stands provides an effective method for investigating the correlation of carbon storages. The CROBAS model has been applied to several coniferous tree species, such as Picea abies L. Karst., Pinus banksiana Lamb., and Chinese pine Pinus tabuliformis [19,20]. The differential evolution (DE) algorithm offers a new and effective approach for meeting the reasonable localization of model parameters. The possible minimum nonlinearity and continuous non-differentiable functions of the population-based DE algorithm enable it to be used by an optimized carbon balance model for forest stand management [21]. The algorithm utilizes an iterative technique where, starting from a randomly generated initial population (or individuals), the values of newly generated optimal individuals are used to replace the entire population (or a part). A relatively simple convergence criterion is the major advantage of this approach [21].

Cunninghamia lanceolata is an important and fast-growing tree species in southwestern China with a wide distribution and long-standing artificial cultivation. Moreover, it plays a vital role in carbon sequestration. Consequently, domestic scholars have investigated the carbon storage properties of C. lanceolata in great detail. However, most studies [22] rarely analyze the carbon allocations of other components in forest ecosystems, including barks, leaves, and roots. Therefore, their results might not reveal the accurate carbon sink features of the entire forest ecosystem. Furthermore, researchers have utilized forestry areas to obtain the forest biomass and then combined it with carbon content and other parameters to calculate the carbon storages. Although this biomass calculation method is simple, it requires huge workloads for investigation, consumes large amounts of time and energy, and is relatively expensive. Therefore, this method is only appropriate for small-scale research [23].

In recent years, several highly accurate and strongly applicable process-based models have been reported [18,20]. These process-based models are gradually transforming from traditional statistical models into ecosystem-mechanism models and provide new possibilities to evaluate the carbon storages and allocations in C. lanceolata plantations. It is important to develop materials of C. lanceolata with the major diameter in the current scenario, which primarily depends on density control [24]. Therefore, there is no doubt that reasonable thinning can effectively promote stand growth [25]. In fact, the selection of an optimal thinning intensity is vital to forest management. However, the effects of thinning on the dynamic response of stand growth and carbon storages in C. lanceolata have been rarely reported.

To evaluate the effects of thinning on the carbon storage in C. lanceolata plantations, we selected sample plots of C. lanceolata plantations with a 6-year survey under different thinning intensities (I: 70%; II: 50%; III: 30%; CK: no thinning) for the carbon balance model analysis. This thinning treatment mainly considers the gradient of thinning intensities, with I, II, and III corresponding to weak thinning, moderate thinning, and heavy thinning, respectively [26]. The DE algorithm was adopted to localize the key parameters of the process-based CROBAS model. Furthermore, a CROBAS-C. lanceolata (CROBAS-CL) model was constructed to simulate the carbon allocations in each part (including stems, leaves, branches, and roots) of a plantation. Based on the relationships between measured stand growth and soil carbon storage, the carbon allocations in the stand were studied, and the correlations and differences between carbon storage in forests and litters were also investigated to reveal the dynamic variation of carbon content in the plantations. Our study provides a sound theoretical reference for the efficient operation and management of C. lanceolata plantations.

2. Materials and Methods

2.1. Study Site

This study was conducted at the Guangxi Youyiguan Forest Ecosystem Research Station, which is a part of the Experimental Center of Tropical Forestry, Chinese Academy of Forestry. The study site is located in the center of Pingxiang City, Guangxi Province (22°02′–22°04′ N, 106°51′–106°54′ E), as shown in Figure 1, and has a southern subtropical humid climate. The annual average temperature ranges from 20.5 to 21.7 °C. In January and July, the average temperature reaches 12.5 °C and 26.0 °C, respectively. The average annual precipitation over the last 5 years is 1386 mm, with 80% of that precipitation occurring from April to September. The annual evaporation capacity is 1275 mm. The altitude ranges from 250 to 800 m. The highest mountain is the Daqing Mountain with an elevation of 1045 m and a gradient of 25–30°. The main types of landforms are hills and foothills. The soil is principally composed of laterite and red soil based on the Chinese soil classification; this soil texture is classified as ferralsols in the World Reference Base for Soil Resources. The soil thickness is generally higher than 100 cm. The main vegetation in this area includes Pinus massoniana, C. lanceolata, and other local broad-leaf tree species (Erythrophleum fordii, Castanopsis hystrix, Betula alnoides, etc.) [27]. We chose C. lanceolata plantations for this investigation based on long-term and continuous research, as the average age of C. lanceolata before thinning is 14 years.

2.2. Data Collection and Sampling Methods

We chose to use C. lanceolata plantations under different thinning intensities as our research object. According to IPCC-Good Practice Guidance for Land Use, Land-Use Change and Forestry (GPG-LULUCF) [28], the sections with favorable and consistent site conditions for tree growth were chosen to conduct the thinning test from 2007 to 2012. The sample plots were square and each plot had an area of 20 × 20 m². Four sample plots were selected to investigate and measure the stand diameter at breast height (DBH), height, and X–Y coordinates. A sample of 123 trees in total was used for analysis in the thinning test at the study site. Specifically, 15, 25, 35, and 48 trees were selected for measurements in each sample plot managed by thinning intensity (heavy, medium, and low) and no thinning, respectively. The C. lanceolata plantation was thinned with the three intensities based on the number of trees removed in each plot of the same area (I: 70%; II: 50%; III: 30%) and one control group (CK, no thinning and density = 1201 stem/ha). Thinning I, II, and III correspond to the remaining tree densities of 375, 625, and 875 stem/ha, respectively.

The crop tree management practices are paramount in realizing close-to-nature forest management. Crop tree is a type of forest management technology that promotes the quality of an individual tree by reducing the canopy competition from interference trees and increasing the growth space of the crop trees. According to the competitive features of individual tree growth, forests were subdivided into crop trees, special crop trees, interference trees, and general trees, which were marked and numbered [29]. Crop trees were indicated with permanent markers. This subdivision method was beneficial in determining the density of crop trees and forest growth and confirmed which of the crop trees were evenly distributed in each sample plot after thinning. All thinning treatments used the traditional thinning-from-below approach [30]. The thinning operation was conducted from October 2007, and all interference trees and some general trees were cut down. Subsequently, information on stand growth was recorded every two years. Meanwhile, soils, litters, and understory plants (shrub and herb) in the sample plots were collected in successive years at a fixed time. Variations in the carbon storage of soils, litters, and understory plants were recorded and analyzed to explore the relationships between carbon content in the stands, understory plants, litter, and soils under different thinning intensities. A typical sampling method was adopted to set the layout of sample plots [31]. All the trees were considered for the analysis based on the measurement of all individuals in the population at different ages. DBH, height, stand densities, and tree growth conditions in the forests were investigated and the number of trees was also recorded [31].

After the field investigation and measurement, the stand characteristics of C. lanceolata plantations obtained under different thinning intensities were recorded and are listed in

Table 1. Stand characteristics of C. lanceolata plantations under different thinning intensities.

Item	Age (Year)	Thinning-I	Thinning-II	Thinning-III	CK Treatment
DBH (cm)	14	19.40 ± 0.59 ab	19.08 ± 1.80 ab	20.33 ± 0.67 a	17.78 ± 1.39 a
	15	21.13 ± 0.43 bc	20.15 ± 1.65 ab	21.98 ± 0.66 b	18.40 ± 1.09 b
	17	23.03 ± 1.79 bc	21.15 ± 2.44 bc	24.33 ± 1.15 c	19.30 ± 1.17 b
	19	27 ± 1.64 abc	25.58 ± 2.55 c	28.08 ± 1.08 d	20.38 ± 1.29 b
Height (m)	14	12.45 ± 1.02 a	12.23 ± 0.40 a	12.30 ± 0.59 a	13.19 ± 1.06 a
	15	13.13 ± 0.70 a	13.28 ± 0.40 b	13.30 ± 0.47 b	13.90 ± 0.75 ab
	17	13.90 ± 0.69 ab	14.58 ± 0.51 c	14.13 ± 0.41 b	14.85 ± 0.70 bc
	19	15.23 ± 0.67 c	16.25 ± 0.81 d	15.50 ± 0.91 c	15.6 ± 0.79 cd
Single-tree biomass (kg per tree)	14	113.05 ± 11.44 a	145.60 ± 9.80 a	188.72 ± 31.04 ab	100.02 ± 21.28 a
	15	107.80 ± 22.77 a	133.88 ± 26.58 a	165.97 ± 38.83 a	113.78 ± 14.22 ab
	17	123.95 ± 10.97 a	160.96 ± 8.50 a	216.02 ± 18.04 ab	136.42 ± 18.53 bc
	19	134.74 ± 35.74 a	161.90 ± 42.14 a	223.72 ± 42.75 b	162.75 ± 19.32 c

DBH, diameter at breast height; Data represent mean ± standard error. Stand characteristics with the same lowercase letters indicate a non-significant difference between them (p = 0.05).

.

According to the collection standards of soil nutrient samples in forests (LY/T 1210–1990), the distribution of each soil’s natural hierarchy and main root system was fully considered. A soil auger was used to drill topsoil samples of 0–20 cm depth. A composite sample (about 1000 g) was synthesized from the samples in each hierarchy. Three soil-sampling points were randomly selected in each sample according to the inverted triangle rule. After mixing them uniformly, the samples were brought to the lab, dried, and granulated for measurement. Soil samples were analyzed according to the International standard methods published by Wilke [32]. The content of organic carbons in the soil was measured using a volumetric method with potassium dichromate oxidation [31].

We chose representative sections with relatively uniform distributions in four corners and central positions for each fixed sampling plot to set up five shrub layers, field layers, and litter layers. Quadrat sampling was used in the harvest method. All shrubs and herbs in the shrub layers (2 m × 2 m), herb layers (1 m × 1 m), and litter layers (1 m × 1 m) were harvested and transported to the lab to be dried at 70 °C. Subsequently, they were weighed to obtain the biomass in each unit area and each sampling plot. To investigate the litters, all those in the herb quadrat were harvested, weighed, and then sampled [33]. All samples were weighed and taken to the lab. The carbon content of each sample was measured using a volumetric method with potassium dichromate oxidation as described by Singh [34]. Additionally, we computed a correlogram using a ggplot2 package (an R-based visual framework) [35] for correlation analysis between carbon content in stand, understory plant, litter, and soil.

2.3. Model and Statistical Analysis

2.3.1. CROBAS-CL Model

CROBAS is a process-based model established on three important concepts. (1) The lateral crown surface area is related to the abnormal velocity of leaves and foliage. (2) Adopting the idea of functional balance between leaves and foliage, it is assumed that the fine root weight is proportional to the foliage weight. (3) The number of leaves present is linearly correlated with the stem’s sapwood cross-sectional area; this is the famous Pipe model theory [18,36]. These concepts can be expressed as follows:

$$w_f=\epsilon A_c^z$$

(1)

$$w_{\gamma }={\alpha }_{\gamma }{w}_f$$

(2)

$$A_i={\alpha }_i{W}_f$$

(3)

where W_f and W_r represent the dry mass of leaves and roots, respectively, A_c denotes the crown surface area, A_i denotes the stem’s sapwood cross-sectional area, z, ε, ${\alpha }_{\gamma }$, and α_i represent the parameters of the model, and i indicates stems, branches, and transportation roots.

The core objective of this model is to estimate the biomass allocations of five carbon reserve sites in the average tree number of a stand. These five parts contain leaves (f), roots (r), stems (s), branches (b), and transportation roots (t). Based on the aforementioned concepts and other assumptions, the biomass of stems, branches, and transportation roots can be expressed as:

$$W_s=p_s{\alpha }_s\left({\varphi }_s{H}_s+{\varphi }_c{H}_c\right)W_f$$

(4)

$$W_b={\rho }_b{\alpha }_b{\varphi }_b{H}_c{W}_f$$

(5)

$$W_t={\rho }_t{\alpha }_t{\varphi }_t\left(H_s+H_c\right)W_f$$

(6)

where W_s, W_b, and W_t represent the dry mass of stems, branches, and transportation roots, respectively, H_s is the height of bare trunk or stalk, H_c is the crown height, ρ_s, ρ_b, and ρ_t represent the densities of stems, branches, and transportation roots, respectively, and α_s, Φ_s, Φ_c, α_b, Φ_b, α_t, and Φ_t are the parameters. The above equations can be used to obtain the biomass allocation in each part of the trees. Furthermore, the entire CROBAS-CL model was constructed by combining these equations with models for photosynthesis, respiration, self-pruning rate, and withered losses, which can be expressed as follows:

$$G=\sum G_i=Y^{-1}\left(P-R\right),\left(i=f,r,s,b,t\right)$$

(7)

$$P=P_0\left(1-e^{-kl}\right)/N$$

(8)

$$R_m=r_1\left(W_f+W_r\right)+r_2\left(W_s+W_b+W_t\right)$$

(9)

where G is the total growth of trees, G_i represents the growth of leaves, roots, stems, branches, and transportation roots, P denotes the photosynthetic capacity, R_m is the respiration rate needed to maintain life processes, R refers to O₂ consumption in respiration, Y is the transformation factor, P₀ is the maximum photosynthetic rate in the unit area, N is the number of trees per square meter, k is the extinction coefficient, l represents the leaf area index, and r₁ and r₂ are empirical parameters.

2.3.2. Localization Parameters

Based on the CROBAS modeling framework of carbon equilibrium, multivariate data, including continuously observed data in sampling plots, tree analysis information, and literature data, were integrated. Regression analysis using SPSS22.0 software (IBM Inc., Chicago, IL, USA) and target optimization technology were used for parametric localization in the C. lanceolata plantations. The relevant details are available at

Table A1. Localization parameters of the CROBAS-CL model.

Parameter	Value	Methods and References
φs	1	Theoretical value based on the pipe model assumption [18]
φc	0.75	Based on the conical form [18]
φb	1	As suggested by Sun and Sheng [52]
φt	1	No heartwood in coarse roots [18]
cb	0.23	From data used in this study (unpublished data)
ct	1	As suggested by Sun and Sheng [52]
ρs, ρb, ρt (kg/m3)	300	Duan et al. 2016 [53]
αs (m2/kg)	0.001	Parameter estimation based on research plots
αb (m2/kg)	0.00068	Parameter estimation based on research plots
αt (m2/kg)	0.00035	Parameter estimation based on research plots
αr	0.0175	Tian et al., 2011 [54]
2z	2.499	Theoretical value based on Mäkelä and Sievänen, 1992 [55]
ξ	0.0492	Calculated from data in Mäkelä and Albrektson, 1992 [56] to give reasonable foliage weight when 2z = 2.499
Y (kg C/kg DWa)	0.65	Y = fc + rg where fc = carbon content of DW = 0.45 [57] and rg = specific growth respiration rate = 0.20 [58]
r1 (kg C/kg DW yrb)	0.2	As suggested by Mäkelä [18]
r2 (kg C/kg DW yrb)	0.02	As suggested by Mäkelä [18]
sf	0.25	Based on needle lifetime of 4a
sr	1	Based on fine-root lifetime of 1a
ds0, db0, dt0	1	By definition as presented in Mäkelä [18]
ds1, db1, dt1	0.01	As estimated by Mäkelä [18]
Ψs	1	Theoretical value according to pipe model
Ψc	0.5	As estimated by Mäkelä [18]
Ψb	0.9	As estimated by Mäkelä [18]
Ψt	0	No heartwood in transport roots
an (m2/kg)	4	Barclay and Trofymow, 2000 [59]
P0	2.134	X(1), parameter optimization within the range 2–4
aσ	0.02	Yoder et al., 1994 [60]
k	0.2	Oker-Blom 1986 [61]
q	1	Trial and error
p	9.7914	X(2), parameter optimization within the range 0–10
aq	0.03912	X(3), parameter optimization within the range 0–1
m0	0.001	As suggested by Mäkelä [18] and confirmed through trial and error
m1	0.01	As suggested by Mäkelä [18] and confirmed through trial and error

φ_s, form factor of sapwood in stem below crown; φ_c, form factor of sapwood in stem within crown; φ_b, form factor of sapwood in branches; φ_t, form factor of sapwood in transport roots; c_b, ratio of crown radius to crown length; c_t, ratio of transport root length to stem length; ρ_s, ρ_b, ρ_t, density of wood; α_s, sapwood area: foliage weight ratio in stem; α_b, sapwood area: foliage weight ratio in branches; α_t, sapwood area: foliage weight ratio in transport roots; α_r, fine root: foliage weight ratio; 2z, “fractal dimension” of foliage in crown; ξ, “surface area density” of foliage; Y, carbon use efficiency; r₁, specific maintenance respiration rate of foliage + fine roots; r₂, specific maintenance respiration rate of wood; s_f, specific senescence rate of foliage; s_r, specific senescence rate of fine roots; d_s0, d_b0, d_t0, specific sapwood area turnover rate per unit relative pruning; d_s1, d_b1, d_t1, specific turnover rate of sapwood area in case of no pruning; Ψ_s, form factor of senescent swapwood in stem below crown; Ψ_c, form factor of senescent swapwood in stem inside crown; Ψ_b, form factor of senescent swapwood in branches; Ψ_t, form factor of senescent swapwood in transport roots; a_n, specific leaf area; P₀, maximum rate of canopy photosynthesis per unit area; a_σ, decrease of photosynthesis per unit crown length; k, extinction coefficient; q, degree of control by crown coverage of self-pruning; p, degree of control by crown coverage of mortality; a_q, parameter related to self-pruning; m₀, specific mortality rate independent of density; m₁, density-dependent mortality parameter.

of the Appendix A.

Based on the methods and references in Table A1 of the Appendix A, a DE optimization model was constructed for three parameters, where X(1), X(2), and X(3) represent P₀, p, and a_q. respectively. P denotes the photosynthetic capacity; P₀ is the maximum rate of canopy photosynthetic in the unit area; α_q is a parameter related to self-pruning from a function of crown coverage obtained by the CROBAS-CL model.

The stand factor of C. lanceolata plantations after thinning at 14 years of stand age was used as the input for 5000 iterations. These three parameters were calibrated within the given parameter range. The error between the height and DBH obtained by the CROBAS-CL model and that obtained by the continuously observed sampling plots was minimized. The optimization model for the three parameters can be expressed as follows:

$$Minxf\left(x,y_0\right)={\displaystyle \sum }_{j=1}^{\left(2\right)}{e}_j$$

(10)

$$s.t.x=\left(x_1,x_2\right)$$

(11)

$$y_0=\left(h_0,d_0\right)$$

(12)

In Equation (10), $e_j=\left\{{\displaystyle {\sum }_{i=1}^{(2)}{[C(t+5i)]}^2}\right./{\left.3\right\}}^{0.5}/{\displaystyle {\sum }_{i=1}^{(2)}R(t+5}i{)}^2/3$, where e represents the predicted error for each variable, j is the number of the predicted variables, i denotes the predicted period, C is the predicted value of each variable obtained by the CROBAS model, R is the measured value, and t represents the initial age of the stand. In Equation (12), y₀ represents the initial state of the stand at t = 14 and h₀ and d₀ represent the mean height and DBH in the initial stand, respectively.

2.3.3. Error Test

We employed root mean square error (RMSE) and mean absolute deviation (MAD) to evaluate the prediction accuracy of the CROBAS model for the measurement of height and DBH, which are defined as follows:

$$RMSE(\%)=100\%\times \frac{\sqrt{{\displaystyle \sum {(H_c-H_Q)}^2/(n-1)}}}{{\displaystyle \sum H_c/n}}$$

(13)

$$MAD=\frac{{\displaystyle \sum \left|H_c-H_Q\right|}}{n}$$

(14)

where H_c and H_Q are the value predicted by the model and measured value, respectively, and n is the number of parameters.

2.3.4. Flowchart of the Process for Exploring the Effects of Thinning on C. Lanceolata Carbon Allocation

The ideas and methods for analyzing the carbon allocation proposed by Mäkelä et al. [18] and Liao et al. [20] were used to determine the effects of thinning on C. lanceolata carbon allocation, as shown in Figure 2. First, based on the CROBAS modeling framework of carbon equilibrium, a DE optimization model was constructed for three parameters, where X(1), X(2), and X(3) represent P₀, p, and a_q. respectively. Then, the trends in average height and DBH between the predicted and measured values for various thinning intensities were built to evaluate the MAD and RMSE of prediction accuracy and reliability obtained by the CROBAS-CL model. DBH and height provided input to the CROBAS-CL model to calculate individual tree volumes for various thinning intensities. Therefore, the carbon allocation in the various organs of C. lanceolata plantations as a function of tree age were calculated under the different thinning intensities obtained by the CROBAS-CL model. Meanwhile, the correlations and differences between carbon storage in forests and litters were also investigated to further reveal the dynamic variations of carbon content in the plantations.

3. Results

3.1. Optimum Parameters

Table 2. Optimum parameters obtained by the CROBAS-CL model.

Fitting Results	Parameters			Object Value
	X(1)	X(2)	X(3)
CK	2.134	9.791	0.039	0.320
I	3.787	8.295	0.089	0.261
II	2.671	8.998	0.098	0.422
III	2.671	8.999	0.099	0.413

X(1)–X(3) are the optimized parameters obtained by implementing the DE algorithm in the localized CROBAS model, which is available in Appendix A.

shows the optimized values of parameters X(1)–X(3), which reveal that the parameters varied with respect to the thinning intensities. The target values (minimized error values) ranged between 0.032 and 0.422. The target values increased with the thinning intensities with an overall surging trend.

3.2. Verification of Accuracy and Reliability of CROBAS-CL

Figure 3 shows the variation of the average height and DBH as a function of forest age for various thinning intensities. These graphs show that after thinning, stand DBH increased with the forest age, but the enhancement rate under different thinning intensities was slightly different. After six years of thinning, these different intensities caused a remarkable variation in DBH. In other words, thinning promoted an increase in the mean DBH of the stand. The DBH for the control group also exhibited a slight increase. The increase in the mean DBH referred to the measured value of thinning intensities following the order Thinning III > Thinning I > Thinning II > CK over time. Compared to the DBH for the control group without thinning, the mean DBH corresponding to Thinning III, II, and I was increased by 43.62%, 30.54%, and 41.36%, respectively. This indicates that the stand densities and strains were adjusted by thinning. The general protocol of thinning was to cut inferior trees and reserve superior trees. Consequently, we observed an obvious increase in the mean DBH for Thinning III and I, especially for Thinning III. The increase corresponding to Thinning II was minimal. The height growth of the C. lanceolata plantations after thinning exhibited a different trend than that exhibited by DBH, as shown in Figure 3. Here, the increase in the mean height with different thinning intensities was manifested as Thinning II > Thinning III > Thinning I over time; however, after six years of thinning, the height for Thinning I was slightly lower than that that for CK, implying that the impact of thinning on height was not obvious. However, compared to the height for CK, the height corresponding to Thinning II and Thinning III was increased by 5.92% and 1.11%, respectively, which implies that the growth of mean height in the sampling plots with different thinning intensities did not exhibit any significant difference.

Table 3. Results of model reliability.

Thinning Type	MAD		RMSE (%)
	H	DBH	H	DBH
CK	0.053	0.056	0.004	0.003
I	0.059	0.159	0.005	0.008
II	0.087	0.134	0.007	0.007
III	0.067	0.165	0.006	0.008

shows the MAD and RMSE (obtained by the CROBAS-CL model) of the predicted height and DBH under different thinning intensities. The MAD of predicted height is less than 0.087, while its RMSE is less than 0.007%. The MAD of predicted DBH is less than 0.165, while its RMSE is less than 0.008%. This indicates that the deviation of height predicted by the CROBAS-CL model from the measured height is less than 0.087 m. Similarly, the deviation between the predicted and measured DBH is less than 0.165 cm. These results confirm that our model is consistent with the precision requirements.

3.3. Carbon Storage and Distribution in Trees

Figure 4 shows that the effect of thinning intensities on the calculated stand volume was not uniform. After implementing the high-intensity Thinning I, the predicted stand volume of C. lanceolata was reduced by 119.62–126.62 m³/hm², as compared to that in the control group CK without thinning. For Thinning III, the stand volume started increasing in the third year after thinning began. Although Thinning II did not have an obvious impact on the stand volume, Thinning levels I and III affected it significantly (p < 0.05).

Figure 5 shows the carbon allocation in the various organs of C. lanceolata plantations as a function of tree age under the different thinning intensities obtained by the CROBAS-CL model. The carbon storage of a single strain in the plantation clearly exhibited an increasing trend. By comparing with the observations in the control group CK, it is evident that thinning promoted the growth of the individual trees. This might be related to the increases in biomass and soil organic content due to stand growth. The carbon storage in various organs decreased in the following order: stems, roots, branches, barks, and leaves. Stems made up the majority of carbon storage in the C. lanceolata plantation, accounting for the biggest allocation proportion (54.59%–55.29%) in carbon storage based on an individual tree. Moreover, Thinning III caused the maximum carbon storage in all the organs except branches, showing that carbon storage in each organ was closely linked with the physiological features of trees. Under the varied thinning intensities, the carbon allocation size of each organ in C. lanceolata followed a similar trend as that for the variation of DBH, i.e., Thinning III > Thinning I > Thinning II > CK. This implies that thinning had a strong impact on the carbon allocations in the stand and promoted the growth of plants. Furthermore, these findings indicate that the thinning intensities had no obvious impact on the growth of stems, branches, leaves, and roots (p = 0.144, 0.111, 0.169, and 0.140, respectively). The carbon content of each organ in arbors was not affected by the thinning intensities.

3.4. Dynamic Variation of Soil Carbon

Figure 6 shows the dynamic variation of soil carbon of the C. lanceolata plantation under different thinning intensities. It is clear here that the soil carbon storage in the absence of thinning remained almost constant, with a minor increase over the years. Under the three levels of thinning, it decreased initially, followed by an increasing trend with respect to the age of trees. In the sixth year of thinning, soil carbon storage in C. lanceolata under the three thinning treatments was obviously higher than that under the control group CK. Compared to the results before thinning, the soil carbon content was increased by 6.3%, 10.3%, and 7.2% for Thinning I, II, and III, respectively. Therefore, the soil carbon content followed the following order: Thinning II > Thinning III > Thinning I > CK. However, the difference in the impact of various thinning intensities on the dynamic variation of soil carbon content was not obvious (p-value = 0.513; F-test = 0.783), and this should be further studied.

3.5. Relationship between Carbon Content in the Stands, Understory, and Soil

Similar to the carbon storage in each organ of the C. lanceolata plantation (Figure 5), the carbon content in the forest vegetation grew for 6 years, and for Thinning I, II and III, it reached 221.39, 205.96, and 233.94 t/hm², respectively, which are significantly higher than 139.50 t/hm² reached by the control group CK. Therefore, carbon storage in the forest vegetation decreased in the following order: Thinning III > Thinning I > Thinning II > CK.

Figure 7 shows the relationship between carbon contents in the stand, understory, and soil. It is clear here that the carbon storages in various locations were positively correlated with other characteristics, except those in litters, which were negatively correlated with other carbon storages. In particular, carbon storage in forest vegetation exhibited extremely high positive correlation with height and DBH (r = 0.71 and 0.97, p < 0.05). Furthermore, soil carbon also revealed extremely high correlation with the stand volume and height (r = 0.89 and 0.75, p < 0.05). However, soil carbon had no obvious correlation with the carbon content in litters and DBH.

4. Discussion

4.1. Effects of Thinning on Stands and Soil

The use of plantations for wood production and carbon absorption to increase carbon storage has garnered significant research attention globally [37]. DBH and height are the chosen indicators of model accuracy verification, because these two factors are often selected and easily measured in forestry surveys. CROBAS-CL model uses functional relationships and carbon balance modeling to describe tree growth and biomass allocation, and simulate the carbon allocations in each part for a plantation by DE algorithm. Therefore, the result of model prediction is the carbon allocations in the various organs of C. lanceolata plantations under different thinning intensities during the following six years of thinning. However, the impact of thinning on carbon allocations of plantation ecosystems needs to be further achieved through correlations and differences between carbon storage in forests, litters and soils, in order to reflect the dynamic variations of carbon content in the plantations.

Among each organ of the C. lanceolata, the stem’s mean carbon storage value under different reserve densities accounts for 55.17% of the total average in the tree phytomass layer. Furthermore, roots also occupy a high proportion, accounting for 33.67% of carbon storages in the arbor layer. After cutting, the preservation of roots in soils can reduce carbon drainage. Branches, leaves, and roots exhibit a high proportion of carbon content. Therefore, logging residue treatment and forest recovery and renewal have a significant impact on carbon storages. After cutting, several branches, leaves, and barks can be used as fuel wood or “tempering mountains”, resulting in the decomposition of various organic matters. The organic matter in the soil surface rapidly conducts oxygen lysis, resulting in the discharge of CO₂ [38]. In this study, the allocation order of carbon storages in each organ of the C. lanceolata arbor layers was observed as stems > roots > branches > leaves. This result might be similar to the allocation order in the majority of plantations globally [39,40].

However, under different thinning levels, the allocation proportion of the carbon storages in the same organ of C. lanceolata as the stand exhibited a significant difference. Generally, the age of trees has a critical impact on the carbon storages of the stand [41,42]. It not only affects the total carbon storages of the stand, but also the allocations of each component within them. This dynamic variation of stand carbon storages as a function of tree age provides an important indicator for forest ecosystems. Particularly, the investigation of dynamic variations of the long-term biomass in ecosystems is of great significance for forest health and sustainable management of plantations [43]. The litter layers have lower carbon storages and the mean value is 1.9–4.4 g/kg. As the carbon storage of litters is primarily determined by regional hydrothermal conditions, higher latitude corresponds to worse decomposition conditions [10]. Furthermore, more litters will accumulate with time. Multiple studies have indicated that soil carbons exhibit a positive correlation with carbon storages of arbors [44]. Our results confirmed this conclusion. Our study could only reveal the carbon storages of each stand component in 2007 and the stand carbon storages during the following six years of thinning because the carbon storage results corresponded to different stages and could not reveal the effect of natural transformation on stand carbon storages. Consequently, using the process-based model, the variation of the stand’s carbon reserve, which was simulated in the next period, was used to predict its growth dynamics.

4.2. Model Evaluation and Uncertainty Analysis

This study is based on the framework of the CROBAS carbon equilibrium model, we conducted parametric localization, parametric calibration, and optimization of the C. lanceolata plantation in Guangxi Province. The CROBAS-CL model meets the requirements of an accurate estimation for height and DBH; thus, the predicted results of carbon reserves matched the accuracy and reliability of the results obtained by statistical test, which implies that the model is useful for the dynamic estimation of carbon reserves (Table 3). The MAD of the predicted height is less than 0.087, while its RMSE is less than 0.007%. The MAD of the predicted DBH is less than 0.165, while its RMSE is less than 0.008%. These results inferred that the model calibration and error tests meet the CROBAS accuracy requirements; the smaller the error value, the higher the model’s accuracy. Theoretically, the disadvantage of model calibration is that it makes model output as dependent on output data as statistical models. On the other hand, it also makes model predictions more reliable, at least under conditions similar to the data. This aspect highlights the importance of datasets for model calibration; a comprehensive dataset covering a wide range of processes and environmental conditions will lead to more robust calibration and the model will be more generally applicable. Model prediction is theoretically possible. However, accurate prediction of the CROBAS-CL model may take more considerations in practice, for example, the life cycle of trees, the relationship between tree species differentiation and naturally sparse process, and the upper limit of C. lanceolata age should also be important influencing factors.

Several studies [45,46] have been conducted over the past decade using “trial and error” or statistical fitting of selected parameters against output variables to calibrate process-based models based on forest inventory and measured data. Blanco et al. [47] and Jyske et al. [48] performed a calibrated process-based forest model on sites over large regions in Canada and Finland, respectively. Then, the effects of thinning on stand growth and diameter distribution were tested against data from permanent growth experiments. Accurate information in forest carbons cannot be obtained from most of the continuously observed sampling plots, unless a destructive experiment is conducted. This is obviously inappropriate under the “theme of protection” [45]. The process-based CROBAS model effectively simulated the carbon dynamics of the stand in the C. lanceolata under different thinning intensities. Moreover, the results indicated that such a prediction was feasible, and it could provide at least some reference values for forest management. Compared to the stand growth model based on relatively traditional biological statistics, the CROBAS-CL model exhibits the advantages of a process-based approach. It not only contains the features of the empirical models and functions of the mechanism models but also explains the reasons for stand growth and allocations from the process mechanism [49,50], which provides an effective means for carbon estimation, as demonstrated in this study. Due to the regulating factors of thinning in plantations, the stand densities were assumed to be unchanged within the considered period in this study. In other words, there was no natural thinning phenomenon, which made it hard to predict the effect of ecological disturbance such as destructive events, including climatic pressure and fire [51]. These factors often cause uncertain impacts in practice. In addition, because C. lanceolata is an endemic species in China, there are few related studies and even less research on carbon distribution. Some model parameters are based on the original parameters from the initial CROBAS model or the data obtained in this experiment. Therefore, further localization of the model parameters is required in future studies to improve the model accuracy.

In this study, we employed the DE algorithm in the CROBAS-CL model for localizing growth module parameters and parameter optimization. The convergence of optimized results depended on the initial value. The application of the DE algorithm improved the validity and accuracy of parameter estimation in the model [17].

5. Conclusions

This study investigated the effects of thinning intensity on the carbon allocation of C. lanceolata plantations. We localized the parameters of the CROBAS model and employed the DE algorithm to optimize the three most important parameters related to total canopy photosynthesis (CO₂ assimilation). Furthermore, we simulated the carbon allocations of each part (stems, leaves, branches, and roots) of the C. lanceolata plantations, and analyzed the correlation and difference in the carbon content in forests and litters between the plantation and measured soil carbon reserves, which revealed the dynamic variation of carbon storages under different thinning intensities. We proved that the CROBAS-CL model accurately predicted height and DBH, and the predicted carbon storages satisfied the requirements of the statistical inspection. Additionally, the model can be useful for the dynamic prediction of carbon reserves in the stand of C. lanceolata plantations.

Compared to C. lanceolata in the contrast group CK, thinning promoted the growth of the remaining trees. The carbon storages of a single strain in the C. lanceolata plantations gradually increased with tree age. Stems had the maximum allocation proportion of carbon storages in the stand of the C. lanceolata and accounted for 54.59%–55.29% of its carbon storages. Moreover, a 30% thinning intensity (Thinning III) corresponded to the maximum carbon storages in other organs, with the exception of branches, showing that carbon storages in each organ were closely related to the physiological features of the tree species. Under different thinning intensities, the allocation of carbon in each organ of the C. lanceolata plantations followed a similar order as the DBH variation, i.e., Thinning III > Thinning I > Thinning II > CK. This implies that thinning has a strong impact on carbon allocation. After thinning, soil carbon storages in the control group CK reduced initially, and then followed an increasing trend. This was related to the environmental conditions and stand growth process. Forest vegetation carbon storage exhibited a remarkable correlation with height and DBH. Furthermore, we correlated soil carbons with the stand volume and height. Overall, we propose that Thinning III (30% thinning intensity) should be used in the study site for C. lanceolata plantations under 20 years or immature timber, as it can promote their growth in southwestern China.