An Alternative Method for Estimation of Stand-Level Biomass for Three Conifer Species in Northeast China

Abstract

Accurate large-scale biomass prediction is crucial for assessing forest carbon storage and dynamics. It can also inform sustainable forest management practices and climate change mitigation efforts. However, stand-level biomass models are still scarce worldwide. Our study aims to introduce the generalized additive model (GAM) as a convenient and efficient approach for forest biomass estimation. Data from 311 sample plots of three conifer species in northeastern China were used to evaluate the performance of the GAM model and compare it with traditional nonlinear seemingly unrelated regression (NSUR) models in predicting stand biomass, including total, aboveground, and component biomass. The results indicated that the goodness of fit of GAM was better than that of NSUR in two model systems. In the majority of cases, the scatter plots and prediction performance revealed that the stand total and component biomass models utilizing GAM outperformed those based on NSUR. Disregarding heteroscedasticity and requiring fewer statistical assumptions provide additional support for the replacement of NSUR-based models with GAM-based models. This study implies that the GAM approach has greater potential for developing a system of stand biomass models.

1. Introduction

Forests play a critical role in terrestrial ecosystems, covering a substantial portion of the terrestrial surface of the Earth. They are crucial in mitigating the human-induced greenhouse effect by absorbing carbon dioxide from the atmosphere [1,2,3]. Since forest biomass is a fundamental characteristic of forest ecosystems, its accurate estimation is essential for experts and policymakers interested in the exchange and storage of carbon on a global scale, nutrient cycling, and energy flow [4,5,6,7]. The global focus on carbon storage estimation was further heightened with the implementation of the Kyoto Protocol under the United Nations Framework Convention on Climate Change. To quantify carbon, assessing forest biomass is a prerequisite. As a result, researchers have been striving to enhance the precision of biomass prediction at the individual tree level [8,9,10,11], and at a regional scale [12,13,14,15].

Generally, forest biomass is estimated using the scaling-up approach, which obtains stand biomass or local scale biomass by adding up the estimated tree biomass [16]. The scale-up strategy uses easily measured variables of the tree, i.e., diameter at breast height, height of the tree, age, and crown factors [17,18]. However, individual tree information is occasionally available for broader management landscapes. An efficient and direct approach to estimating forest biomass involves establishing a correspondence between forest monitoring data and biomass projections at a comparable spatial level. A link provides easily accessible stand variables to develop the total stand and component biomass models [19]. The models are usually divided into two types, including stand variables and volume-derived biomass. The volume-derived biomass models use the stand volume as a single predictor, while multiplying the obtained stand volume with biomass expansion factors allows for the calculation of a stand’s total or component biomass [20,21,22]. The establishment of additive individual tree biomass models in different regions and stand conditions worldwide has been crucial in supporting the development of stand biomass models [23,24,25,26,27,28,29].

Furthermore, national monitoring of forests is based on national forest inventories (NFIs), whose data generally summarize stand attributes, such as basal area and stand density. This approach offers rapid and convenient conditions to facilitate the development of stand biomass models. Nonetheless, compared to studies on individual tree biomass models, there are still relatively few tree species and areas available [30,31].

An essential and logical feature for predicting biomass is additivity, i.e., when multiple component models are constructed for individual trees or stands, the total predicted biomass should be equal to the sum of the predicted values of each component, ensuring consistency in the estimation process [32]. Therefore, various estimation methods have been proposed by researchers to ensure the additivity and statistical efficiency of biomass system parameters [19,33,34,35,36]. Among the methods suggested, nonlinear seemingly unrelated regression (NSUR) is a popular parameter estimation method because it can ensure additivity by fitting the system simultaneously and considering the constraints of linear correlation between the components. The flexibility and general nature of NSUR allow the component model to use different independent variables, and each component model may contain a particular weight function for heteroscedasticity in the system. As the result, the total biomass produces a low variance [37].

Additionally, if NSUR is expected to estimate subtotal biomass (crown, aboveground, ect.), Parresol [35] proposes the addition of one or more constraints to maintain the additivity among the components while guaranteeing the additivity to the total biomass. Therefore, NSUR has been widely applied in crown width models, individual tree models, and stand biomass models [11,38,39,40].

In cases where prediction is the primary goal of developing a model, parametric models can deliver accurate future predictions [41]. However, sometimes the relationship between explanatory and response variables might be too complex to find an appropriate mathematical function [42]. Semiparametric models, e.g., generalized additive models (GAM), provide an attractive substitute for such scenarios.

The advantage of GAM is the automatic identification of the amicable relationships between the dependent and predictor variables [43]. Thus, GAM does not face the problem of comparing several potential alternative model forms. Another merit of GAMs is their flexibility and robustness [44], as GAMs are more relaxed toward the basic assumptions required in developing linear or nonlinear models (i.e., independent, normality, and equal variance). Therefore, GAMs have been used extensively in forestry studies like taper modeling [45], insect pests [46], increment of basal area [47], taper functions [41,48], and tree biomass [9].

Korean pine (Pinus koraiensis Siebold & Zucc.), a multipurpose tree species, has substantial economic value. It mainly produces good-quality timber and edible nuts. The nuts of Korean pine are highly priced in the international market [49]. Thus, Korean pine contributes to the maintenance of ecosystem diversity as well. Korean larch (Larix olgensisi A. Henry) is highly valued for its luxuriant growth and resistance to diseases and insects, making it an important commercial tree species in the north, northeast, and south subalpine regions of China as well as a valuable indigenous species. The wood is highly valued for its excellent mechanical properties and is widely used for various applications, including furniture, flooring, housing, and plywood. It is among the four pulp-producing coniferous species in China [50]. Mongolian pine (Pinus sylvestris var. mongolica Lity.), is highly resistant to soil infertility, drought, and cold. It has been extensively planted in the successful afforestation project in the Three-North regions of China for soil and water conservation, windbreaks, and timber production [51,52]. The proportion of planted forests in Jilin Province is around 12%, covering an area of about 1,759,400 hectares, and continues to expand [53]. Planted forests have relatively high productivity, which is also essential for increasing forest cover and fixing carbon dioxide in the atmosphere. Therefore, accurately estimating the biomass of plantation forests is important for understanding the mechanisms of the global carbon cycle and formulating policies to mitigate global warming [54].

The aims of this study were (1) to introduce GAM to establish stand biomass prediction systems, i.e., stand biomass models including observed stand variables (Model 1) and volume-derived models (Model 2) for three types of plantation species in the northeastern region of China; and (2) to evaluate the predictive ability of two modeling techniques (GAM and NSUR) for stand total and component biomass prediction.

2. Materials and Methods

2.1. Experimental Site and Data Source

The research was fulfilled in a forest plantation that is considered representative (125°3′–125°50′ E, 42°18′–43°14′ N), administered by the Dongfeng County Forestry Bureau, Jilin Province, Northeast China. The study area covers 48,486.79 hectares of forest plantation, dominated by Korean pine (4396.73 ha), Korean larch (27,021.1 ha), and Mongolian pine (17,068.96 ha). Dongfeng County is situated in the southern region of the Haddaling Mountains, with an elevation ranging from 300 to 914 m. The topographic features are mainly hilly terraces, accounting for 61.4% of the county’s area. The region is situated in the monsoon zone and classified as having a humid, mid-temperate climate. The mean annual temperature ranges from −37 °C to 35 °C, with a mean of 4.5 °C. The annual precipitation fluctuates from 451.9 mm to 867.5 mm, and the frost-free period is about 128 days. There are twelve forest farms in the Dongfeng Forestry Bureau. The study utilized data collected from seven forest farms, namely: Daxing, Hengdaohe, Renhe, Shahezhen, Yangmulin, Yimianshan, and Zhongyu. The experiment site and the spatial distribution of the sample plots are shown in Figure 1.

In the summer and autumn of 2021, a total of 311 temporary plantation plots were established using a stratified sampling method, with site types as the stratification variable. Specifically, 121 plots for Korean pine, 90 for Korean larch, and 100 for Mongolian pine. Out of the total number of plots, 214 were classified as pure plantation forests, with the dominant species accounting for 100%. The 97 plots are mixed forests, with 20 plots having dominant tree species accounting for 50%–80% of the composition, while the remaining 77 plots have dominant tree species accounting for over 80% of the composition. The main species included in the mixed forest plots are Quercus mongolica, Fisch. ex Ledeb., Betula platyphylla Sukaczev, Fraxinus mandshurica Rupr., Ulmus pumila L., Juglans mandshurica Maxim., and Acer elegantulum W.P. Fang and P.L. Chiu. The plot size varied from 100 to 900 m².

Site information, including elevation, longitude, latitude, slope, and aspect, was recorded for each plot. All tree information was measured and recorded in each plot, except for those with a diameter at breast height less than 5 cm, a height less than 1.3 m, or those that were dead. The forked trees were considered a single tree when the fork’s location was above 1.3 m or two trees when the fork’s location was below 1.3 m [55]. The measurements taken for sampled trees in each plot included diameter at breast height, tree height, the spatial location of the tree, the height of the first living whorl, and crown width in four directions (east, west, north, and south). In a total of 311 sample plots, we measured 13,140 trees, with 4976 trees measured in Korean pine plots, 3518 trees measured in Korean larch plots, and 4646 trees measured in Mongolian pine plots.

2.2. Stand Characteristics and Biomass Estimation

Stand characteristics were calculated using the individual tree variables. For example, the basal area (G, m²·ha⁻¹) and stand quadratic mean diameter (Dg, cm) were calculated through $\frac{\sum \frac{\pi }{4}{D}^2}{A}$ and $\sqrt{\sum D^2/n}$. Where D is diameter at breast height, A is the area of the plot, and n is the number of trees per plot, respectively. The stand mean height (Hm, m) was calculated as mean tree height measurements within the plot, while stand density (N, trees·ha⁻¹) was obtained dividing the number of trees by the plot area. The stand volume was derived through adding up individual volume of all trees within a plot, and then dividing that sum by the area of plot. The individual stem volume was estimated employing the model of the tree volume table of China (Table S1), which contains tree volume models for three species in Dongfeng County. The next step involved calculating the total and component biomass of the stand. Firstly, the biomass of each tree was obtained using a previously reported tree biomass model (Table S2). To obtain the total biomass of a plot, we sum up the individual tree biomass within the plot. Secondly, the total biomass of a plot was determined by dividing the plot area by the stand biomass.

Table 1. Descriptive statistics for variables and biomass components.

Variables	Korean Pine				Korean Larch				Mongolian Pine
	Min.	Max.	Mean	SD	Min.	Max.	Mean	SD	Min.	Max.	Mean	SD
G (m2·ha−1)	18.82	42.52	29.65	5.18	11.00	35.76	24.71	5.19	20.18	45.80	31.20	5.41
Dg (cm)	7.61	35.78	21.07	4.86	8.42	27.50	18.99	4.09	11.95	23.90	18.61	2.80
Hm (m)	5.08	18.18	13.29	2.41	7.02	24.90	17.39	3.63	8.92	18.96	14.80	2.18
N (trees·ha−1)	350.00	4375.00	988.77	561.95	400.00	2625.00	981.34	452.92	450.00	3200.00	1244.33	481.30
V (m3·ha−1)	78.97	266.49	172.57	36.38	55.94	284.11	185.83	48.66	109.92	264.86	183.61	32.54
Age (a)	19.00	65.00	45.72	11.12	14.00	58.00	39.73	14.60	18.00	49.00	33.96	7.58
Bt (Mg·ha−1)	65.82	224.97	140.84	32.25	35.10	235.74	144.53	44.32	83.21	208.13	143.72	26.06
Ba (Mg·ha−1)	50.28	181.12	111.91	26.58	28.58	185.76	114.45	34.63	65.83	171.71	118.67	22.12
Br (Mg·ha−1)	15.54	43.85	28.92	5.82	6.53	49.97	30.09	9.70	15.84	36.42	25.04	4.40
Bs (Mg·ha−1)	40.33	117.92	77.55	15.61	24.02	167.15	101.97	31.73	51.36	140.13	95.39	18.17
Bb (Mg·ha−1)	4.87	48.49	22.41	8.62	4.36	15.63	10.72	2.38	8.49	21.20	14.61	2.67
Bn (Mg·ha−1)	4.32	20.13	11.95	3.24	0.19	3.92	1.76	0.74	5.21	12.84	8.67	1.55

Note: G represents basal area (m²·ha⁻¹); Dg represents stand quadratic mean diameter (cm); Hm represents stand mean height (m); N represents stand density (trees·ha⁻¹); V represents stand volume (m³·ha⁻¹); Age represents stand age (a); and B_t, B_a, B_r, B_s, B_b, and B_n represent the component biomass of total, aboveground, belowground, stem, branch, and needle (Mg·ha⁻¹).

summarizes the statistical information for the variables as well as the total and component biomass.

2.3. Stand Biomass Model Specification

Some previous researchers have indicated a strong relationship between stand biomass and stand variables [56,57,58]. Furthermore, Castedo-Dorado et al. [12] included the variable of stand volume to explore efficient stand biomass modeling for large management units. Therefore, we evaluated two alternate models for stand biomass prediction: (1) The system of stand total and component biomass models that used observed stand variables (Model 1), and (2) The model system that used stand volume as a predictor (Model 2). The base models are as follows:

$$B={\alpha }_0{X}_1^{{\alpha }_1}{X}_2^{{\alpha }_2}{X}_3^{{\alpha }_3}\cdots X_n^{{\alpha }_n}+\epsilon$$

(1)

$$B={\alpha }_0{V}^{{\alpha }_1}+\epsilon$$

(2)

where B is the total or component biomass, ${\alpha }_0$−${\alpha }_n$ are the estimated parameters of the models, X₁ − X_n are stand variables (stand age, stand density, basal area, and stand mean height); V is standing volume; and $\epsilon$ is error term of the model.

2.3.1. NSUR-Based Stand Biomass Model

An attempt was made to develop a stand biomass model using stand variables, and scatter plots of stand variables and stand biomass were plotted (Figure 2). The results showed that the scatter trend of biomass with G, Hm, and N was more obvious than with other variables. The modeling attempts using stand variables showed that the former two were the best predictor variables for estimating stand total and component biomass. To ensure the additivity of the models, an additive model system was applied, containing six equations that simultaneously fit the system with cross-equation constraints by Nonlinear Seemingly Unrelated Regression (NSUR) [35]. Therefore, Model 1 employs the basal area and stand mean height as predictor variables, as shown in Equation (3), while Model 2 uses stand volume as the only predictor variable, as shown in Equation (4).

$$\{\begin{array}{l}{B}_i=c_{ij}{G}^{k_{ij}}H{m}^{m_{ij}}+{\epsilon }_i\\ B_a=B_s+B_b+B_n+{\epsilon }_a\\ B_t=B_r+B_s+B_b+B_n+{\epsilon }_t\end{array}$$

(3)

$$\{\begin{array}{l}{B}_i=c_{ij}{V}^{k_{ij}}+{\epsilon }_i\\ B_a=B_s+B_b+B_n+{\epsilon }_a\\ B_t=B_r+B_s+B_b+B_n+{\epsilon }_t\end{array}$$

(4)

where B_i is the biomass of i components; G is the basal area; Hm is stand mean height; V is standing volume; c_ij, k_ij, and m_ij are parameters of the model system; i indicates r, s, b, or n (belowground, stem, branch, or needle); a and t are aboveground and total; and j indicates Korean pine, Korean larch, or Mongolian pine.

The NSUR-based stand biomass models were developed by the SAS/ETS MODEL procedure.

2.3.2. GAM-Based Stand Biomass Model

The generalized additive model (GAM) is a semi-parametric expansion of the generalized linear model [48,59,60]. GAM, a flexible semi-parametric method, explores complex nonlinear and non-monotonic relationships between the response and predictor variables. The model is generally expressed as:

$$Y=\beta +{\displaystyle \sum _{n=1}^p{S}_n\left(X_n\right)}+\epsilon$$

(5)

where Y is the response variable with an exponential family distribution; X_n (n = 1…p) are predictor variables; S_n (·) (n = 1…p) are smoothing functions; and $\beta$ is the intercept. The GAM follows the basic assumption that smoothing functions are additive to each other. It implies that the implementation of GAM is similar to the addition of multivariate linear regression functions. The use of smooth functions may fit better, tapping the potential of the data better than a purely parametric model because it is sometimes data-driven rather than model-driven. The GAM is usually fitted by a numerical algorithm, i.e., backfitting. The minimization of the smoothing function is sought through Equation (6):

$${\epsilon }^2={\left[Y-\left(\beta +{\displaystyle \sum _{n=1}^{p}S\left(X_n\right)}\right)\right]}^2$$

(6)

Various smooth functions are available for GAM, allowing different smooth functions to be combined by the additive. Different smooth functions include thin-plate regression spline (TP), B-spline (BS), Gaussian process smooths (GP), cubic spline (CS), polynomial cubic spline (PS), and locally weighted scatter plot smoothing (LOESS) [61]. If the spline function is not set at modeling, then the default thin-plate regression spline will be used [62]. This ensures consistency with the simultaneous fitting of the stand total and component biomass models using NSUR. Two types of GAM-based stand biomass modeling systems were established. We chose the same stand variables, i.e., basal area (G), stand mean height (Hm), and stand volume (V), to develop the semi-parametric stand biomass model systems. Therefore, the prediction systems (Model 1) utilizing basal area (G) and stand mean height (Hm) are established according to Equations (7)–(9) and have the following format:

$${\widehat{B}}_i={\widehat{\beta }}_{ij}+S\left(G\right)+S\left(Hm\right)$$

(7)

$${\widehat{B}}_a={\widehat{B}}_s+{\widehat{B}}_b+{\widehat{B}}_n$$

(8)

$${\widehat{B}}_t={\widehat{B}}_r+{\widehat{B}}_s+{\widehat{B}}_b+{\widehat{B}}_n$$

(9)

The second system (Model 2) used stand volume (V) and comprised the following Equations (10)–(12).

$${\widehat{B}}_i={\widehat{\beta }}_{ij}+S\left(V\right)$$

(10)

$${\widehat{B}}_a={\widehat{B}}_s+{\widehat{B}}_b+{\widehat{B}}_n$$

(11)

$${\widehat{B}}_t={\widehat{B}}_r+{\widehat{B}}_s+{\widehat{B}}_b+{\widehat{B}}_n$$

(12)

where ${\widehat{B}}_i$ is i component biomass estimates, ${\widehat{\beta }}_{ij}$ are estimated parameters, S (·) are smoothing functions; other variables have been previously defined.

The GAM-based stand biomass models were fitted using the R software (Version 4.2.0, Vienna, Austria) package mgcv [61].

2.4. Heteroskedasticity

Heteroscedasticity is often present in the tree or stand biomass modeling process, which reflects a particular pattern in the residuals of the estimated base model, i.e., the variance of the residuals increases or decreases as the predicted values increase [63]. The issue of heteroscedasticity does not hinder the development of GAM-based stand biomass models since the underlying assumptions of GAM are not constrained by this problem. However, the presence of heteroscedasticity can compromise parameter estimates [64], lead to inaccurate confidence intervals [65], and introduce systematic errors in predictions [66,67]. The weight functions were introduced to overcome the problem. Firstly, stand biomass models were initially fitted using NSUR [68]. In the current state, the residuals with heteroscedasticity were obtained for each model. The variance of residuals could be modeled with one or more stand-alone variables [16,35]. Accordingly, the assumption was as follows:

$$E\left({\epsilon }_i^2\right)={\sigma }^2\left(X_1^{{\phi }_1}\cdots X_n^{{\phi }_n}\right)$$

(13)

Secondly, the estimated residuals squared (${\widehat{\epsilon }}_i^2$) of each model were transformed using the natural logarithm, as well as the stand variables. Then they were fitted by stepwise regression of the following form:

$$\ln \left({\widehat{\epsilon }}^2\right)=\ln \left({\sigma }^2\right)+{\phi }_1\ln \left(X_1\right)+\cdots +{\phi }_n\ln \left(X_n\right)$$

(14)

where $\widehat{\epsilon }$ represents unweighted residuals for each model, ${\phi }_1$− ${\phi }_n$ are parameters, ${\sigma }^2$ represents the variance of residuals, and X₁ − X_n represents stand variables.

Thirdly, the parameters of Equation (14) were selected based on a significance level of α = 0.05. The weight functions were $1/\left(X_1^{{\widehat{\phi }}_1}\cdots X_n^{{\widehat{\phi }}_n}\right)$ for different stand biomass models. Finally, the model system was refitted using the NSUR in SAS Procedure, and the $resid. B_i=resid.B_i/\sqrt{\left(X_1^{{\widehat{\phi }}_1}\cdots X_n^{{\widehat{\phi }}_n}\right)}$ (where resid. B_i represents the residual of ith model) was added in the process [69,70].

2.5. Evaluation of Systems of Stand Biomass Models

The performance of Model 1 and Model 2 based on GAM and NSUR methods was evaluated using coefficient of determination (R²) and root mean square error (RMSE) as metrics for evaluating goodness-of-fit. Additionally, the models were validated using mean error (ME), RMSE, and relative root mean square error (RRMSE). The mathematical expressions for these statistics are as follows:

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n{\left(B_i-{\widehat{B}}_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(B_i-\overline{B}\right)}^2}}$$

(15)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(B_i-{\widehat{B}}_i\right)}^2}}{n-1}}$$

(16)

$$ME=\frac{1}{n}{\displaystyle \sum _{i=1}^n(B_i-{\widehat{B}}_i)}$$

(17)

$$RRMSE=\frac{\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(B_i-{\widehat{B}}_i\right)}^2}}{n-1}}}{\overline{B}}$$

(18)

where B_i represents observed values, ${\widehat{B}}_i$ represents predicted values, $\overline{B}$ represents the mean of the observed values, and n represents the number of samples. The units for RMSE and ME are Mg ha⁻¹, while R² and RRMSE are dimensionless.

To validate the optimal model system, we employed the leave-one-out cross-validation (LOOCV) technique. This technique leaves individual sample plots in each step. We deleted a sample plot from the data and executed the model fitting with the remaining dataset. In each stage, the calculated parameters were employed to forecast the biomass in the removed sample plot. The same procedure was carried out for all 311 sample plots. Then, we used the measured and predicted values to calculate validation statistics and verify the performance of the models [48,71,72].

In addition, to evaluate the performance of the NSUR and GAM models in predicting the outcome, we utilized scatter plots to compare the predicted values with the observed values. Furthermore, we employed a simple linear model to fit the predicted values against the observed values, facilitating a clear comparison between the two methods. The linear model adopted the form y = kx + b, where y is the observed value, x is the predicted value, k is the slope, and b is the intercept. When the intercept of the linear model approaches zero and the slope approaches 45 degrees, it indicates a closer match between the predicted values and the observed values, signifying a more accurate prediction by the model. This approach offers the advantage of visually identifying the differences in predictions between the methods, particularly for subtle variations. As a result, it has found practical applications among forestry researchers [73,74].

3. Results

3.1. Stand Biomass Models Fitting

While fitting Model 1 and Model 2 using the NSUR method, the inherent correlation among stand biomass components was considered. Additionally, two constraints were considered: The total biomass is the combined biomass of its individual components, whereas the aboveground biomass is comprised of the collective biomass of its stems, branches, and needles. Thus, three 6 × 6 matrices were assumed for three conifer species. Figure 3 displays the matrices showing the correlation of residuals between the total biomass and the individual component biomasses for Model 1 and Model 2. The residual correlations of stand totals and components of Korean larch were higher than those of the other two species. For Korean pine, significant correlations existed between the total biomass and its aboveground biomass, as well as between the total and branch biomass. For Mongolian pine, the correlations were stronger between the stand total and aboveground biomass, stem biomass, and root biomass.

By employing all available data, Model 1 and Model 2 were fitted to ensure that the aggregate of the predicted values for biomass components equaled the predicted values of the total biomass. In order to address the heteroscedasticity in Model 1, different weight functions were employed (

Table 2. Goodness-of-fit statistics and weight functions of stand biomass models for three species.

Model Type	Method	Components	Korean Pine			Korean Larch			Mongolian Pine
			R2	RMSE	Weight Function	R2	RMSE	Weight Function	R2	RMSE	Weight Function
Model-1	NSUR	Total	0.9428	7.7168	G2.6765	0.9261	12.0525	G1.1228	0.9250	7.1377	G4.4215
		Aboveground	0.9299	7.0347	G2.7411	0.9302	9.1551	G0.9810	0.9002	6.9891	G4.3542
		Belowground	0.9738	0.9429	G4.2025Hm−3.3574	0.9104	2.9025	G0.2832	0.9595	0.8846	G−2.2436Hm5.7752
		Stem	0.9515	3.4382	G2.0790	0.9248	8.7041	G0.9381	0.8774	6.3642	G4.8456
		Branch	0.8058	3.7971	G4.9491	0.9899	0.2387	G2.6160Hm−2.6810	0.9260	0.7265	G4.3753
		Needle	0.8839	1.1051	G2.9662	0.6451	0.4435	G3.7972Hm−3.1957	0.9680	0.2777	G−3.0660Hm4.9095
Model-1	GAM	Total	0.9510	7.1376	—	0.9278	11.9050	—	0.9304	6.8725	—
		Aboveground	0.9400	6.5088	—	0.9318	9.0403	—	0.9079	6.7126	—
		Belowground	0.9767	0.8895	—	0.9124	2.8698	—	0.9603	0.8766	—
		Stem	0.9554	3.2972	—	0.9263	8.6145	—	0.8865	6.1239	—
		Branch	0.8345	3.5058	—	0.9901	0.2363	—	0.9317	0.6981	—
		Needle	0.8953	1.0491	—	0.6618	0.4329	—	0.9683	0.2766	—
Model-2	NSUR	Total	0.9853	3.9155	V3.3461	0.9740	7.1482	V1.8583	0.9895	2.6644	V0.1614
		Aboveground	0.9752	4.1908	V2.6208	0.9773	5.2144	V1.5499	0.9708	3.7786	V0.3408
		Belowground	0.9793	0.8390	V−0.4416	0.9600	1.9403	V1.6494	0.8783	1.5342	V−0.0474
		Stem	0.9765	2.3933	V−0.0536	0.9714	5.3676	V1.3832	0.9512	4.0169	V0.3475
		Branch	0.7986	3.8668	V−2.4178	0.9521	0.5207	V−0.4525	0.9760	0.4135	V−0.6777
		Needle	0.9197	0.9197	V−1.5208	0.7228	0.3920	V−0.6429	0.8533	0.5946	V−0.4920
Model-2	GAM	Total	0.9879	3.5504	—	0.9749	7.0253	—	0.9906	2.5300	—
		Aboveground	0.9791	3.8407	—	0.9781	5.1306	—	0.9732	3.6238	—
		Belowground	0.9798	0.8285	—	0.9616	1.9011	—	0.8790	1.5297	—
		Stem	0.9774	2.3481	—	0.9725	5.2644	—	0.9548	3.8641	—
		Branch	0.8203	3.6524	—	0.9555	0.5017	—	0.9760	0.4136	—
		Needle	0.9249	0.8887	—	0.7303	0.3866	—	0.8535	0.5942	—

).

Table 2 presents indicators of model fitness for stand biomass models developed with both NSUR and GAM for three species. The GAM-based model systems showed superior fitting to the NSUR-based model systems. As a whole, Model 2 reflected larger values of R² (0.7303–0.9906) and smaller values of RMSE (0.4136–7.1482). Although Model 1 showed competitive values of R² (0.6618–0.9901), it reflected significant variations in terms of RMSE (0.2363–12.0525). However, the models behaved differently in terms of the component biomass of the species analyzed. With both NSUR and GAM methods, the performance of Model 1 was better for branch biomass of Korean pine and Korean larch. Similarly, Model 1 showed better results for the belowground and needle biomass of Mongolian pine. The results also revealed that Model 1 and Model 2 delivered the highest RMSEs for total, aboveground, belowground, and stem biomass of Korean larch.

3.2. Stand Biomass Model Validation

Table 3. Validation of stand biomass models for three conifer species.

Model Type	Method	Components	Korean Pine			Korean Larch			Mongolian Pine
			ME	RMSE	RRMSE	ME	RMSE	RRMSE	ME	RMSE	RRMSE
Model-1	NSUR	Total	−0.5504	7.8563	5.5783	0.9524	12.4345	8.6033	0.1155	7.2955	5.0763
		Aboveground	−0.4118	7.1647	6.4020	0.7250	9.4467	8.2544	0.0340	7.1503	6.0252
		Belowground	−0.1386	0.9548	3.3011	0.2274	2.9931	9.9482	0.0814	0.8951	3.5742
		Stem	−0.0273	3.5093	4.5252	0.7044	8.9834	8.8101	0.0689	6.5092	6.8235
		Branch	−0.2745	3.8622	17.2346	0.0043	0.2467	2.3020	−0.0154	0.7437	5.0897
		Needle	−0.1099	1.1231	9.3948	0.0162	0.4564	25.8887	−0.0194	0.2812	3.2436
Model-1	GAM	Total	−0.0099	7.5668	5.3728	0.0222	12.3306	8.5314	−0.0229	7.2490	5.0439
		Aboveground	−0.0102	6.9119	6.1761	0.0178	9.3666	8.1844	−0.0212	7.0872	5.9721
		Belowground	−0.0003	0.9283	3.2095	0.0044	2.9693	9.8691	−0.0017	0.9077	3.6247
		Stem	−0.0085	3.5047	4.5194	0.0154	8.9198	8.7477	−0.0198	6.4655	6.7768
		Branch	−0.0011	3.7245	16.6200	0.0026	0.2506	2.3384	−0.0033	0.7375	5.0470
		Needle	−0.0006	1.1019	9.2177	−0.0002	0.4477	25.3993	0.0019	0.2895	3.3393
Model-2	NSUR	Total	−0.0173	3.9802	2.8261	0.2464	7.2727	5.0319	−0.0536	2.7190	1.8919
		Aboveground	−0.0298	4.2552	3.8022	0.1848	5.3045	4.6350	−0.0479	3.8606	3.2532
		Belowground	0.0124	0.8540	2.9528	0.0616	1.9748	6.5636	−0.0056	1.5688	6.2644
		Stem	−0.2630	2.4398	3.1461	0.2447	5.4619	5.3565	−0.0529	4.1005	4.2986
		Branch	0.1904	3.9274	17.5256	−0.0388	0.5306	4.9521	−0.0001	0.4258	2.9139
		Needle	0.0429	0.9354	7.8252	−0.0211	0.3982	22.5869	0.0051	0.6103	7.0403
Model-2	GAM	Total	−0.0016	3.7532	2.6649	0.0044	7.1727	4.9627	−0.0116	2.6935	1.8742
		Aboveground	−0.0025	4.0415	3.6113	0.0036	5.2381	4.5770	−0.0160	3.8290	3.2265
		Belowground	0.0010	0.8429	2.9142	0.0007	1.9412	6.4519	0.0044	1.5829	6.3207
		Stem	0.0074	2.4107	3.1085	0.0029	5.3750	5.2713	−0.0151	4.0687	4.2652
		Branch	−0.0097	3.8200	17.0463	0.0011	0.5128	4.7852	−0.0036	0.4371	2.9913
		Needle	−0.0003	0.9207	7.7021	−0.0003	0.3943	22.3655	0.0026	0.6141	7.0838

exhibits the assessment measures for the stand biomass models through the application of LOOCV. The values of ME for GAM-based models were noticeably lower than those for NSUR-based models. The only exception was the branch component in NSUR-based Model 2 for Mongolian pine. As a whole, the GAM-based stand total and the majority of the component biomass models delivered lower values of RMSE and RRMSE compared with the NSUR-based models. However, for the NSUR method, Model 1 performed slightly better in the branch component for Korean larch. Additionally, the RMSE and RRMSE of Model 1 and Model 2 (NSUR) were slightly lower in belowground and needle and belowground and branch components, respectively, for Mongolian pine. It illustrated that Model 1 and Model 2 systems based on GAM could provide reliable prediction accuracy for stand total and component biomass. Like the fitting outcomes, the overall prediction accuracy of Model 2 was noticeably better than Model 1. Moreover, the accuracy of the models for predicting the stand total and aboveground, belowground, and stem components was lower for Korean larch than for Korean pine and Mongolian pine. This results also showed larger values of RRMSE for stand branch or needle components of Model 1 and Model 2 with both NSUR and GAM methods. The larger values of RRMSE were generated because of smaller average values of stand branch and needle than other components, according to Equation (18) above.

3.3. Prediction Accuracy of GAM and NSUR Methods

The NSUR and GAM modeling techniques were evaluated using scatter plots of the predicted values against the observed values for Model 1 and Model 2 (Figure 4 and Figure 5). In terms of the stand biomass models for Korean pine, the GAM-based Model 1 outperformed the NSUR method, with an intercept closer to 0 and a slope closer to 45 degrees (Figure 4). For Korean larch, the GAM-based Model 1 showed better performance for the needle and belowground components. Although the NSUR method had a slightly better slope for the remaining components, the GAM method maintained its superiority in terms of the intercept. In the case of Mongolian pine, the GAM method exhibited slightly poorer performance, except for the belowground and needle components. When considering Model 2, the GAM-based models outperformed the NSUR method in terms of both the intercept and slope for all three species (Figure 5).

The predicted biomass values of stand total and aboveground acquired through LOOCV were further evaluated using box plots to compare the residuals of the stand component models’ sum (Figure 6). Both GAM and NSUR indicated larger prediction errors for Korean larch. However, GAM performed better in mean prediction errors, and NSUR showed lower variance in Model 1. The residual variance was almost similar in GAM and NSUR for Mongolian pine. Regarding Korean pine, GAM residuals were closer to 0 in Model 2, and its overall error distribution was similar. As a whole, GAM and NSUR reflected similar results for stand total and aboveground biomass predictions.

3.4. Parameter Estimates

Table 4. Parameter estimates and standard errors (SE) of the NSUR-based stand biomass models for three conifer species.

Model Type	Components	Korean Pine			Korean Larch			Mongolian Pine
		ci1	ki1	mi1	ci2	ki2	mi2	ci3	ki3	mi3
Model 1	Belowground	0.5676 ** (0.0122)	1.0160 ** (0.0068)	0.1898 ** (0.0060)	0.1378 ** (0.0121)	1.0948 ** (0.0284)	0.6488 ** (0.0347)	0.8164 ** (0.0223)	1.0083 ** (0.0082)	−0.0180 * (0.0067)
	Stem	1.6851 ** (0.1178)	0.9574 ** (0.0227)	0.2259 ** (0.0210)	0.5606 ** (0.0437)	1.0805 ** (0.0260)	0.6020 ** (0.0311)	1.4037 ** (0.1161)	0.8019 ** (0.0195)	0.5426 ** (0.0240)
	Branch	0.0425 ** (0.0052)	1.0562 ** (0.0378)	1.0357 ** (0.0359)	0.3374 ** (0.0121)	1.0075 ** (0.0104)	0.0791 ** (0.0112)	0.2616 ** (0.0189)	0.8657 ** (0.0174)	0.3889 ** (0.0196)
	Needle	0.0923 ** (0.0058)	1.0199 ** (0.0204)	0.5448 ** (0.0206)	0.0033 * (0.0010)	1.2843 ** (0.0734)	0.7405 ** (0.1026)	0.2676 ** (0.0085)	1.0555 ** (0.0095)	0.0078 ns(0.0076)
Model 2	Belowground	0.2053 ** (0.0095)	0.9607 ** (0.0088)		0.0438 ** (0.0028)	1.2478 ** (0.0125)		0.2266 ** (0.0269)	0.9029 ** (0.0223)
	Stem	0.6402 ** (0.0218)	0.9322 ** (0.0065)		0.1815 ** (0.0094)	1.2096 ** (0.0101)		0.3665 ** (0.0260)	1.0668 ** (0.0135)
	Branch	0.0110 * (0.0032)	1.4739 ** (0.0560)		0.1553 ** (0.0082)	0.8121 ** (0.0101)		0.0692 ** (0.0038)	1.0266 ** (0.0104)
	Needle	0.0267 ** (0.0032)	1.1837 ** (0.0230)		0.0014 * (0.0004)	1.3603 ** (0.0546)		0.0753 ** (0.0087)	0.9106 ** (0.0220)

Note: (**): parameter estimates significant at p < 0.0001; (*), parameter estimates significant at p < 0.05; (ns), non-significant parameter estimates.

and

Table 5. Parameter estimates, and standard errors (SE) of the GAM-based stand biomass models for three conifer species.

Model Type	Components	Korean Pine	Korean Larch	Mongolian Pine
		$\mathit{\beta }$	$\mathit{\beta }$	$\mathit{\beta }$
Model-1	Belowground	28.9231 ** (0.0816)	30.0867 ** (0.3025)	25.0433 ** (0.0877)
	Stem	77.5495 ** (0.3035)	101.9671 ** (0.9081)	95.3928 ** (0.6195)
	Branch	22.4097 ** (0.3237)	10.7153 ** (0.0249)	14.6117 ** (0.0705)
	Needle	11.9541 ** (0.0966)	1.7628 ** (0.0456)	8.6684 ** (0.0277)
Model-2	Belowground	28.9231 ** (0.0753)	30.0867 ** (0.2004)	25.0430 ** (0.1530)
	Stem	77.5495 ** (0.2135)	101.9671 ** (0.5549)	95.3928 ** (0.3911)
	Branch	22.4100 ** (0.3360)	10.7153 ** (0.0529)	14.6117 ** (0.0414)
	Needle	11.9541 ** (0.0816)	1.7628 ** (0.0408)	8.6684 ** (0.0594)

Note: (**): parameter estimates significant at p < 0.0001.

present the estimated parameters and corresponding standard errors (SE) for Model 1 and Model 2 based on NSUR and GAM for the species analyzed. The majority of estimated parameters were statistically significant at p < 0.05. After testing various smooth functions, the final selection for the GAM-based models included Gaussian process smooths (GP) for the components of Korean pine, thin-plate regression splines (TP) for the components of Korean larch, and B-splines (BS) for the components of Mongolian pine.

Model 1 showed positive powers of G for all stand biomass components across the species, and positive powers of Hm for all components except for belowground biomass in Mongolian pine. These results suggest that an increase in G for the same Hm led to an increase in stand biomass. Likewise, the increase in Hm for the same G increased the biomass for all components of the species, excluding the belowground biomass of Mongolian pine, which decreased. In Model 2, all parameters for stand biomass components were positive. This indicates that there is a positive relationship between V and both stand total and component biomass, meaning that an increase in V leads to an increase in biomass (Table 4).

4. Discussion

Forest management and assessment studies require quantitative information on trees and stands. In forestry, a forest stand refers to a contiguous cluster of trees that share similar attributes such as age distribution, species composition, site quality, and structure. It is considered a fundamental management unit in forestry [19,75]. Therefore, precise estimation of stand biomass is of utmost importance for carrying out scientific research on ecosystem productivity, carbon cycles, and nutrient and energy transfer at the local, regional, or global scale.

This study developed two alternate models (Model 1 and Model 2) for predicting the stand biomass of Korean pine, Korean larch, and Mongolian pine. Model 1 contained commonly used stand variables, while Model 2 was based on stand volume. Both models were evaluated using GAM and NSUR methods for predicting total biomass and its components. We employed the LOOCV to verify the performance of models in estimating the stand’s total and component biomass. Overall, the modeling technique of GAM was better than NSUR, and Model 2 was superior to Model 1 for the species analyzed.

4.1. The Development of Stand Biomass Model Systems

Stand biomass models link different biomass components with commonly measured variables in forest inventories and, thus, serve as an essential tool for forest biomass estimation. The variables documented in such inventories generally include basal area, stand volume, and stand density. This study found basal area (G) as an essential variable in the model system of stand variables for the three conifer species, which supports previous studies [12,76]. However, to improve the accuracy of stand biomass models, Xin et al. [77] have recommended incorporating stand mean height (Hm) as a secondary variable. In agreement with their findings, stand mean height was significant statistically for the majority of the stand components models in this study (Table 4).

Introducing additional stand variables may improve the prediction accuracy of the two model systems in this study, but a simplified model system may be more practical and applicable for forest managers, even though introducing additional stand variables may enhance the prediction accuracy of the two model systems in this study [78]. Thus, the biomass prediction selects the appropriate number of variables with reasonable accuracy to fulfill management requirements and minimize the workload and cost of field surveys [79,80]. Moreover, introducing multiple independent variables in a model might create complexity in the modeling process. Therefore, two parsimonious model systems were developed in this study.

Researchers have used different linear and non-linear methods to develop additive biomass model systems in recent decades. The additive biomass model includes simultaneous equations and separate modeling. Some methods used in this model involve predicting based on component proportions, such as proportion adjustment, fractional multinomial logit regression, Dirichlet, three-step proportional weighting system, and log-ratio regression approach. [33,37,81,82]. The parameter estimation methods used include ordinary least squares, seemingly unrelated regression, Dirichlet regression, and two error-in-variable models [25,38,83]. When using these methods, it is important to consider using the same variables for fractional multinomial logit regression [81].Using different variables may prevent the simplification of the model structure in the three-step proportional weighting system. The NSUR-based simultaneous equations were chosen for stand biomass modeling and compared with the GAM-based system due to its applicability, flexibility, and popularity. We excluded the linear model that logarithmically transforms stand biomass and stand variables due to the potential systematic bias introduced by anti-log transformations [84,85]. Additionally, the use of correction factors can make it challenging for the system to maintain additivity in predicted values.

Based on the statistical criteria, Model 2 exhibited superior predictive power compared to Model 1 in estimating all biomass components across three conifer species (Table 3). This outcome can primarily be attributed to the strong correlation between volume and stand biomass (Figure 2). Another factor contributing to this superiority is the higher proportion of stand stem biomass to total biomass, as observed by Bi [86]. Additionally, He et al. [56] discovered that the stand biomass model utilizing stand volume outperformed models based on other stand variables for Larix plantations in China. Regarding stand volume as the only predictor. Similarly, researchers have reported findings for various tree species in Northeast China, with stand volume as the sole predictor [57].

It is crucial to note that the stand volume was estimated using the individual tree volume model, which means that propagation errors were expected in the estimation of stand biomass [87]. In this study, it was observed that the needle component of Korean larch had the lowest R² value (0.65–0.73) among the species. This occurred because the leaf biomass was smaller compared to the other two species at the individual tree level (Figure 2). The relatively smaller leaf biomass could lead to more discrete observed values if there were variations in the number of trees within the sample plot (Figure 4 and Figure 5).

4.2. GAM vs. NSUR

The comparison of GAM and NSUR methods in this study showed that the former had better performance for three conifer species and can be widely applied for stand biomass estimation. The goodness-of-fit and validation statistics showed that, in terms of RMSE and RRMSE values, the GAM-based Model 1 and Model 2 outperformed NSUR in most biomass models for all species (Table 2 and Table 3). In particular, the GAM-based models delivered the lowest ME values for stand total and component biomass (Table 3). The minimizing loss function in NSUR simultaneously minimizes the loss of nonlinear biomass models of the system and accounts for the inherent correlations between total, subtotal, and component biomass. However, the minimizing loss function might sacrifice the accuracy of the component models, although it ensures additivity [37,88]. Even if the GAM-based component model was developed separately, the predicted values of the total or subtotal were acquired through aggregating the predicted values of the components. However, the prediction accuracy of GAM’s stand total and subtotal biomass was better than NSUR, and most GAM-based component models had better statistical indicators in Model 1 and Model 2 for the species analyzed (Table 3). The data-driven GAM-based stand biomass models indicate that it can fit complex trends in the data than model-driven, which is supported by tree biomass prediction and taper modeling [9,41,48].

The scatter plots reveal slight variations between the NSUR and GAM techniques, as evidenced by the values of the slope and intercept. In general, the GAM-based models demonstrated greater accuracy in estimating the total and component biomass of the analyzed species (Figure 4 and Figure 5). At the same time, the variance of the stand total residuals was similar for the NSUR and GAM methods (Figure 6), although the latter performed more consistently in the overall evaluation.

Stand biomass modeling is convenient with the GAM technique due to fewer restrictions on statistical assumptions. The NSUR may face problems of model selection and heteroskedasticity. In such cases, the GAM can serve as an alternative that is flexible and robust for modeling non-linear and non-constant variance structures [89]. The data-driven GAM deals with the relationship between the response and predictor variables a suitable relationship through smoothing functions. One concern about GAMs is over-fitting, which happens when a model predicts noise as a pattern. GAM automatically sets a suitable degree of smoothness for splines during the modeling process [9]. Furthermore, LOOCV was applied to minimize the approximate out-of-sample prediction error and prevent over-fitting [90]. The GAM predictions should be biologically realistic to avoid flagrant errors such as negative values for stand biomass predictions. We did not find such a deviation in this study (Figure 4 and Figure 5).

The NSUR method, unlike GAM, considers the inherent correlation between biomass components in the additivity system (Figure 3). Therefore, NSUR has better statistical efficiency in the component model estimation of the system [34,35]. Moreover, the error term of total and aboveground biomass was acquired by summing up the error terms of the components (Equations (3) and (4)). It should be noted that the error terms may result in a singular variance and covariance matrix across equations. This situation may occur when either no weight function was applied or the same weight function was utilized for the entire system [91,92]. Due to the use of an optimization iterative algorithm (Gauss-Newton) in NSUR, this algorithm heavily relies on the initial parameter values provided. If the initial iterated values for the parameters differ significantly from the estimated model parameters, it is possible that the parameter estimation may not converge, especially when dealing with a large number of model parameters. Therefore, when using NSUR, it is necessary to carefully set the initial iterated values based on experience or literature to ensure parameter convergence [35,63].

Finally, this study demonstrated that GAM-based prediction systems could reduce the slight uncertainties in estimating stand total and component biomass. However, the slight uncertainties might be amplified at the stand level when stand biomass is extrapolated to forest or national levels. Therefore, with the main aim of stand biomass prediction, GAM could serve as a flexible and promising alternative modeling approach.

5. Conclusions

Two stand biomass models (Model 1 and Model 2) were developed for three conifer species in Northeast China using the GAM and NSUR methods in this study. The developed stand biomass models included estimates of the total, aboveground, and component biomass. The fitting and validation statistics and the graphical analysis showed that the GAM-based models were more appropriate in predicting stand total and component biomass and provided a reliable alternative to the traditional NSUR method. In addition to better predictions, the GAM requires fewer statistical assumptions that induce efficient forestry data modeling and convenient predictions. The volume-derived model (Model 2) was superior to the stand biomass model, including observed stand variables (Model 1). Therefore, the GAM-based Model 2 was recommended for the biomass estimation of three conifer species. However, the GAM-based Model 1 also provides an alternative where stand volume is unavailable. Future work focusing on different species in the region would broaden the scope of this research to a larger landscape.