Comparison of Tree Biomass Modeling Approaches for Larch ( Larix olgensis Henry) Trees in Northeast China

: Accurate quantiﬁcation of tree biomass is critical and essential for calculating carbon storage, as well as for studying climate change, forest health, forest productivity, nutrient cycling, etc. Tree biomass is typically estimated using statistical models. Although various biomass models have been developed thus far, most of them lack a detailed investigation of the additivity properties of biomass components and inherent correlations among the components and aboveground biomass. This study compared the nonadditive and additive biomass models for larch ( Larix olgensis Henry) trees in Northeast China. For the nonadditive models, the base model (BM) and mixed e ﬀ ects model (MEM) separately ﬁt the aboveground and component biomass, and they ignore the inherent correlation between the aboveground and component biomass of the same tree sample. For the additive models, two aggregated model systems with one (AMS1) and no constraints (AMS2) and two disaggregated model systems without (DMS1) and with an aboveground biomass model (DMS2) were ﬁtted simultaneously by weighted nonlinear seemingly unrelated regression (NSUR) and applied to ensure additivity properties. Following this, the six biomass modeling approaches were compared to improve the prediction accuracy of these models. The results showed that the MEM with random e ﬀ ects had better model ﬁtting and performance than the BM, AMS1, AMS2, DMS1, and DMS2; however, when no subsample was available to calculate random e ﬀ ects, AMS1, AMS2, DMS1, and DMS2 could be recommended. There was no single biomass modeling approach to predict biomass that was best for all aboveground and component biomass except for MEM. The overall ranking of models based on the ﬁt and validation statistics obeyed the following order: MEM > DMS1 > AMS2 > AMS1 > DMS2 > BM. This article emphasized more on the methodologies and it was expected that the methods could be applied by other researchers to develop similar systems of the biomass models for other species, and to verify the di ﬀ erences between the aggregated and disaggregated model systems. Overall, all biomass models in this study have the beneﬁt of being able to predict aboveground and component biomass for larch trees and to be used to predict biomass of larch plantations in Northeast China.


Introduction
Plantation forests typically have high growth rates and thereby absorb large amounts of carbon dioxide and help mitigate global climate change. Larch (Larix olgensis Henry) is an important tree species for afforestation and acquiring commercial timber in Northeast China. This species is the We aimed to explore the difference in biomass predictions of the BM, MEM, AMS1, AMS2, DMS1, and DMS2. Therefore, the objectives of our study were: (1) to construct two nonadditive biomass models (i.e., BM and MEM) of each component based on DBH only or both DBH and H; (2) to construct four additive biomass models (i.e., AMS1, AMS2, DMS1, and DMS2) based on D only or both D and H with weighted NSUR; (3) to verify the performance of the biomass models with jackknife resampling; and (4) to compare the accuracy of model fitting and performance of the different nonadditive and additive biomass models.

Data
Data used in this study are obtained from a large dataset of trees sampled from larch plantation stands in Heilongjiang Province, one of the largest forestry provinces in Northeast China. The sample plots were established in July/August of 2007-2016. A total of 40 plots were selected from 10 sites ( Figure 1). Each plot was 20 × 30 m or 30 × 30 m in size. For each of the 40 plots, 1-2 sample trees for each of the dominant, codominant, intermediate, and suppressed types were analyzed, and thus, 5-7 trees were sampled from each plot. In total, the dataset consisted of measurements on 229 larch trees using destructive biomass sampling and included the following parameters: diameter at breast height (DBH), total tree height (H), green weights of stems, branches and foliage, green weights of the sampled disks, branches, and foliage, and the moisture content of disks, branches, and foliage. The specific biomass determination method including field and laboratory measurements followed similar protocols as Dong et al. [5]. For each sampled tree, the sum of branches, foliage, and stem dry biomass was calculated as the aboveground biomass. A summary of the descriptive statistics for DBH and H, as well as the aboveground and component biomass of trees is listed in Table 1. The stem, branch, and foliage biomass of all sample trees as well as their relationships with DBH and H are shown in Figure 2.
Forests 2020, 11,202 3 of 21 (2) to construct four additive biomass models (i.e., AMS1, AMS2, DMS1, and DMS2) based on only or both and with weighted NSUR; (3) to verify the performance of the biomass models with jackknife resampling; and (4) to compare the accuracy of model fitting and performance of the different nonadditive and additive biomass models.

Data
Data used in this study are obtained from a large dataset of trees sampled from larch plantation stands in Heilongjiang Province, one of the largest forestry provinces in Northeast China. The sample plots were established in July/August of 2007-2016. A total of 40 plots were selected from 10 sites ( Figure 1). Each plot was 20 × 30 m or 30 × 30 m in size. For each of the 40 plots, 1-2 sample trees for each of the dominant, codominant, intermediate, and suppressed types were analyzed, and thus, 5-7 trees were sampled from each plot. In total, the dataset consisted of measurements on 229 larch trees using destructive biomass sampling and included the following parameters: diameter at breast height ( ), total tree height ( ), green weights of stems, branches and foliage, green weights of the sampled disks, branches, and foliage, and the moisture content of disks, branches, and foliage. The specific biomass determination method including field and laboratory measurements followed similar protocols as Dong et al. [5]. For each sampled tree, the sum of branches, foliage, and stem dry biomass was calculated as the aboveground biomass. A summary of the descriptive statistics for and , as well as the aboveground and component biomass of trees is listed in Table 1. The stem, branch, and foliage biomass of all sample trees as well as their relationships with and are shown in Figure 2.

Base Model
In general, is used as the first variable in tree biomass models because it is simple and easy to apply. To improve the performance of biomass models, is the second variable incorporated into the model. According to the data obtained by the visual inspection of stems, branches, and foliage, biomass data can be modeled by a power function with and ( Figure 1). Thus, two BMs were used to establish independent component biomass models. These BMs are given by: where represents the stem, branch, and foliage biomass in kilograms (k = s, b, and f, for stem, branch, and foliage, respectively); and represent the diameter at breast height and total tree height, respectively; , , and represents the model parameters; and is the model error term.

Mixed Effects Model
Considering the structure of biomass data, two MEMs with the sample plots as a random factor can be used to establish independent component biomass models. These MEMs are given by: where , , and are the random effect parameters in ith sample plot ( 1,2,3 … , ), and is the error term. In the equation, obeyed a normal distribution with zero mean and within sample plot variance and covariance matrix , and of the sample plot random effects ( ,

Base Model
In general, DBH is used as the first variable in tree biomass models because it is simple and easy to apply. To improve the performance of biomass models, H is the second variable incorporated into the model. According to the data obtained by the visual inspection of stems, branches, and foliage, biomass data can be modeled by a power function with DBH and H ( Figure 1). Thus, two BMs were used to establish independent component biomass models. These BMs are given by: where W k represents the stem, branch, and foliage biomass in kilograms (k = s, b, and f, for stem, branch, and foliage, respectively); DBH and H represent the diameter at breast height and total tree height, respectively; β k0 , β k1 , and β k2 represents the model parameters; and ε k is the model error term.

Mixed Effects Model
Considering the structure of biomass data, two MEMs with the sample plots as a random factor can be used to establish independent component biomass models. These MEMs are given by: where u ik0 , u ik1 , and u ik2 are the random effect parameters in ith sample plot (i = 1, 2, 3 . . . , m), and ε ik is the error term. In the equation, ε ik obeyed a normal distribution with zero mean and within sample plot variance and covariance matrix R, and u ik of the sample plot random effects (u ik0 , u ik1 and u ik2 ) obeyed a normal distribution with zero mean and sample plot variance and covariance matrix D.

Aggregated Model Systems
Aggregated model systems fit the component biomass data simultaneously, which explicates the instinctive correlations among the component biomass of the same sample tree. There are two structures for an aggregated model system, i.e., those with one constraint (AMS1) and those with no constraint (AMS2). •

Aggregated model systems with one constraint
Following the model structure specified in Parresol [30], AMS1 ensures additivity between the tree aboveground biomass and component biomass with one constraint: tree aboveground biomass is equal to the sum of the tree component biomass. The expressions of AMS1 with DBH only and combination of DBH, H are given as: where a, s, b, and f denote the aboveground, stem, branch, and foliage, respectively.

Aggregated model systems with no constraint
Following the model structure specified in Affleck and Diéguez-Aranda [31], the expressions of AMS2 with DBH only and a combination of DBH and H are:

Disaggregated Model Systems
Disaggregated model systems can ensure biomass additivity by directly partitioning tree aboveground biomass into three components (W s , W b , and W f ) in this study. There are also two structures for the disaggregated model systems, i.e., those without an aboveground biomass model (DMS1) and those with an aboveground biomass model (DMS2).

Disaggregated model systems without an aboveground biomass model
For DMS1, the given tree aboveground biomass (Ŵ a−DBH andŴ a−DBH,H ) are from the fitted model W a−DBH = eβ a0 DBHβ a1 andŴ a−DBH,H = eβ a0 DBHβ a1 Hβ a2 , respectively. The expressions of DMS1 with DBH only and a combination of DBH and H are as follows: where β b2 β s2 = α 2 , and β f 2 β s2 = γ 2 are model parameters to be estimated.

Disaggregated model systems with an aboveground biomass model
For DMS2, the tree aboveground biomass model (W a−DBH and W a−DBH,H ) and component biomass models were fitted simultaneously. The expressions of DMS2 with DBH only and a combination of DBH and H are as follows:

Parameter Estimation Methods
In our study, the model parameters in the BM and MEM were estimated using nonlinear ordinary least squares (OLS) and restricted maximum likelihood (REML) methods, respectively. For AMS1, AMS2, DMS1, and DMS2, the model parameters were estimated using nonlinear seemingly unrelated regression (NSUR). More detailed information on the theory and algorithm of the three methods can be found in the literature [7,[21][22][23].
All the above models were fitted using different procedures, such as PROC NLIN, PROC NLMIXED, and PROC MODEL in SAS 9.3 software [33].

Weighting Function for Heteroscedasticity
Because of the heteroscedasticity in the model residuals shown by the tree biomass data, to correct the heteroscedasticity of variance, a power variance function of the following equation var(ε ik ) = σ 2Ŵ 2φ ik was used in the tree biomass models, where ε ik is the mixed model residual, φ is a parameter to be estimated,Ŵ ik is the estimated biomass of k component of each tree on the ith sample plot using fixed part of the mixed effects model, and σ 2 is a scaling factor for error dispersion [23,35].

Model Evaluation and Validation
As described above, the six biomass modeling approaches evaluated in this study were: (1) BM, (2) MEM, (3) AMS1, (4) AMS2, (5) DMS1, and (6) DMS2. It is well known that the quality of model fitting does not entirely reflect the quality of future prediction, so that model validation is necessary to assess and evaluate the predictive quality of the different biomass models. Several authors suggest that the most applicable models should be validated using a Jackknifing technique, also known as the "leave-one-out" method or Predicted Sum of Squares (PRESS) [5,7,9]. Thus, in this study, the biomass equations were fitted to the entire dataset (sample trees N of all plots), while model validation was accomplished by Jackknifing technique in which the biomass models were constructed using all but some sample trees of one sampled plot data (fitting data included the sample trees of m − 1 sample plots), and then the fitted models were used to predict the value of the dependent variable for some sample trees of the excluded sample plot. The model fitting was assessed by two goodness-of-fit statistics (Equations (13) and (14)), and the model performance was evaluated by four model validation statistics of Jackknifing (Equations (15)-(18)) as follows: The coefficient of determination (R 2 ): The percent root mean square error (RMSE%): and the predictive performance (mean prediction error (MPE), mean prediction error percent (MPE%), mean absolute error (MAE), and mean absolute percent error (MAE%) of the tree biomass prediction equation: where W ijk is the observed k (k = s, b, f, and a, represented stem, branch, foliage and aboveground, respectively) component biomass value of jth tree in ith sample plot;Ŵ ijk is the estimated biomass value of jth tree in ith sample plot of the model to which were fitted all the observations of m sample plots; W k is the average value of the biomass;Ŵ ijk,−i is the predicted value of the model to which was fitted all the observations of m−1 sample plots; n j was the number of observations in ith sample plot; and m the number of sample plots. It is important that the parameters of random effects were predicted using the best linear unbiased predictions (BLUPs) method for MEM [36,37]. A vector of random effects parameters of sample plot i was calculated using: where the element of the Z ik matrix was the partial derivative of the nonlinear function with respect to its random parameters;ε ik was obtained by the difference between the observed and predicted biomass using the model, including only fixed effects;D was the variance-covariance matrix estimated in the modelling process; andR was the corresponding variance-covariance matrix of within sample plot.

Comparison of Different Biomass Modeling Approaches
In this study, once the biomass estimate was calculated for each component, the component estimates were summed up to generate the aboveground biomass estimate. We used the six modeling approaches to calculate the biomass of each tree component, and we compared the difference among these modeling approaches using the values of R 2 , RMSE, MPE, MPE%, MAE, and MAE%.

BM and MEM Based on DBH and H
The stem, branch, and foliage biomass equations were fitted independently by weighted OLS and REML. The parameter estimates and their standard errors for the BMs and the MEMs are listed in Table 2. The parameters estimates were significantly different from 0 in each biomass equation (Table 2). For both BM and MEM, the parameter β i1 of DBH was positive for each biomass component, while the parameter β i2 of H was positive for the stem component and negative for the branch and foliage components. These data show that stem biomass increased with an increase in H; however, the branch and foliage decreased with an increase of H for the same DBH. For BM and MEM, the aboveground biomass should be estimated by summing the estimated stem, branch, and foliage biomass.  (3) and (4)]. β k0 , β k1 and β k2 : fixed parameters; u ik0 , u ik1 and u ik2 : random effect parameters; σ 2 u ik0 : variance of u ik0 ; σ 2 u ik1 variance of u ik1 ; σ 2 u ik2 : variance of u ik2 σ u ik0 σ u ik1 : covariance of u ik0 and u ik1 ; σ 2 : residual variance.

Model Variables
Components

AMS1 and AMS2 Based on DBH and H
The parameter estimates of the stem, branch, and foliage biomass models in AMS1 and AMS2 were fitted simultaneously by weighted NSUR and listed in Table 3. The parameters β i0 , β i1 , and β i2 were highly significant, and their negative or positive prescriptions were consistent with those values in BM and MEM, and they had biological significance. Compared to AMS1, there was no constraint of W a = W r + W s + W b + W f and three equations (i.e., stem, branch, and foliage equations) in AMS2. Thus, the aboveground biomass was also estimated by summing the estimated stem, branch, and foliage biomass. For AMS2, a constant 3 × 3 matrix was assumed for cross-correlations among all three equations, and a 4 × 4 matrix was assumed for AMS1.
For example, in this study, there were strong correlations between aboveground and stem biomass, and between branch and foliage components in AMS1 and AMS2 based on DBH, as shown in the following two correlation matrices:

DMS1 and DMS2 Based on DBH and H
The only difference between DMS1 and DMS2 lies between the tree aboveground biomass model being fitted independently or simultaneously. The parameter estimates for DMS1 and DMS2 were fitted by weighted NSUR and listed in Table 4. All the other parameter estimates except γ 1 were highly significant in each biomass equation. Similarly, a 3 × 3 and 4 × 4 matrixes for cross-correlations in DMS1 and DMS2 based on DBH are given by:

DMS1
Stem As with AMS1 and AMS2, there were high correlations between the aboveground and stem biomass, between the branch and foliage components, and weak correlations between the aboveground and branch components and between the aboveground and foliage components.  Table 5 shows the R 2 , RMSE%, and weight functions for each biomass equation developed using the six different biomass modeling approaches. The results indicated that a majority of the biomass equations fit the biomass data well, with R 2 > 0.70. All models fit the aboveground and stem biomass data best, while small R 2 values were observed in the branch and foliage biomass equations. Figure 3 shows the scatterplots of the observed aboveground and component biomass as well as the model predictions by the six biomass modeling approaches. Table 5 and Figure 3 indicated that the MEM with DBH only and combination of DBH and H yielded higher R 2 and smaller RMSE% than did the other five biomass modeling approaches, and the addition of H improved the goodness-of-fit for the majority of the aboveground and component biomass equations. The DMS1 and AMS2 produced a slightly larger R 2 value for most of the aboveground and component biomass models, compared with AMS1 and DMS2, and had a very close values of R 2 and RMSE% for predicting the aboveground and component biomass (Table 5). Table 5. Goodness-of-fit statistics and the parameter of weight functions (φ) for the biomass models developed using the BM, MEM, AMS1, AMS2, DMS1, and DMS2.

Approaches
Model Furthermore, we used the jackknifing technique to assess the validity of these biomass models. The model validation statistics were computed and presented in Table 6, in which MPE and MPE% represented the average prediction error, and MAE and MAE% represented the magnitude of prediction error. The model prediction error varied across the different biomass modeling approaches. Model jackknife statistics indicated the stem biomass models of MEM based DBH and DMS1 based combination of DBH and H slightly overestimated the tree stem biomass (MPE < −0.5 kg and MPE% < −0.3%), and the other biomass models underestimated stem, branch, foliage, and aboveground biomass (MPE > 0.3 kg and MPE% > 0.3%) ( Table 6). For all biomass modeling approaches, the average prediction error of the branch biomass was the largest (MAE% was between 19%~39%) and the average prediction error of the aboveground biomass was the smallest (MAE% was between 6%~16%). Biomass models based on a combination of DBH and H seemed preferable to biomass models based on DBH only.
The validation statistics suggested that MEM was better than BM, AMS1, AMS2, DMS1, and DMS2 at predicting the aboveground and component biomass; no significant differences among between BM, AMS1, AMS2, DMS1, and DMS2 were observed (Table 6). To better compare the different biomass modeling approaches, prediction across the diameter classes would be a good method to validate tree biomass models. In this study, the MPE, MPE%, MAE, and MAE% in tree components of different biomass modeling approaches are based on DBH only, and the combination of DBH and H are shown in Figure 4. The figure indicates that the six biomass modeling approaches with DBH only and combination of DBH and H fitted well for the aboveground and component biomass for larch. Furthermore, the MEM did better than the other five biomass modeling approaches for most diameter classes and adding H into biomass models reduced the prediction error in all classes, especially in the largest diameter class (D > 25 cm).
Overall, based on R 2 , RMSE, MPE, MPE%, MAE, and MAE% for the aboveground and component biomass, the MEM was the best for the aboveground and each component biomass, and the DMS1 was slightly better than BM, AMS1, AMS2, and DMS2 for most component biomass. Compared with the AMS1, AMS2 decreased the most of MPE, MPE%, MAE, and MAE% for aboveground, stem, branch, and foliage biomass, and it would be better for estimating the aggregated models. Similarly, DMS1 was better for estimating the disaggregated models.

Discussions
For the selection of biomass model variables, DBH is an indispensable predictor of biomass models. In practice, tree biomass models constructed with only DBH require basic forest inventory data in their application [10,20,38]. Our results showed that DBH was the primary explanatory variable in the component biomass models. This may originate from the intimate correlations between components and tree diameter [7,10,20,39]. However, within the given DBH, there is usually some variation among the aboveground and component biomass values, which highlight that it is insufficient to predict the aboveground and component biomass by the biomass model constructed with only DBH. Thus, to improve the prediction accuracy of the aboveground and component biomass models, another variable should be added into the biomass model [20]. An increasing number of scholars have often considered H as another commonly used and vital predictor variable to reduce biased estimates of biomass models because tree height usually reflects the site factors [15,18]. In our study, to simulate the biomass allometric relationships for larch trees, W k = e β k0 DBH β k1 and W k = e β k0 DBH β k1 H β k2 were selected to construct basic equations. Six biomass modeling approaches were constructed and validated with the jackknifing technique. Our findings were consistent with previous studies, i.e., that a combination of DBH and H significantly improved the prediction accuracy of aboveground and component biomass [6,7,9,40].
Tree biomass models are categorized as nonadditive or additive models. Nonadditive models cannot synchronously consider the aboveground and component biomass data, leading to unequal aboveground biomass. The additive biomass models comprise a desirable characteristic for a system of equations used for the tree biomass prediction, which explicate the instinctive correlations among component biomass of the same sample, and thus, they have a great statistical efficiency [23,30,41]. In this study, we applied six biomass modeling approaches to develop the biomass models. The BM and MEM have separately fitted aboveground and component biomass models, and they did not hold the additivity property for aboveground biomass. Therefore, the sums of the predictions of tree component models were usually larger or smaller than the predictions of the aboveground biomass models, although the differences were small. In contrast, AMS1, AMS2, DMS1, and DMS2 successfully accounted for the correlations among component biomass by a covariance matrix, in which the aboveground biomass prediction was aggregated from the predictions of the tree component models or disaggregated into tree component biomass. Thus, the AMS1, AMS2, DMS1, and DMS2 held the additivity property for the aboveground biomass.
The AMS1 and AMS2 was fitted with independent nonlinear biomass models, in which there is no random effect in each model. The AMS2 contains no constraint, while AMS1 contains one constraint that guarantees that aboveground predictions will be exactly equal to the sum of the biomass prediction of the stem, branch, and foliage component. Our results indicated that both AMS1 and AMS2 fitted the data and performed well in terms of the average prediction errors for aboveground and component biomass predictions. In this study, we demonstrated the differences between AMS1 and AMS2 because of the constraints imposed on the model system. Furthermore, the AMS1 and AMS2 in our study had smaller standard errors of parameters compared to BM, although AMS1, AMS2, and BM possess similar R 2 , RMSE%, MPE, MPE%, MAE, and MAE%, which was consistent with results from Parresol [30]. In comparison to the BM, the AMS1 and AMS2 accounted for correlations between component biomass and focus on additivity. Therefore, we recommend AMS1 and AMS2 over the BM. In addition, AMS2 was indeed better for predicting the aboveground and component biomass than AMS1, even though AMS1 actually uses aboveground biomass model as a dependent equation. These data are in accordance with those from Zhao et al. [23]. Thus, the AMS2 is more suitable to construct the biomass models for aggregated model systems.
Both DMS1 and DMS2 maintained the properties of additivity for the aboveground and component biomass. In our study, the prediction accuracy of the aboveground and each component biomass model using DMS1 were higher than those from using DMS2. This is likely because disaggregated model systems depend on the aboveground biomass model, which is commonly thought to be the most accurate among the aboveground and component biomass models [7]. Although DMS1 was slightly superior to AMS2, the advantage of AMS2 over DMS1 lies in the fact that it has been successfully implemented for individual biomass estimation, and that it is more maneuverable in practical applications.
The results in Table 5, Table 6, and Figure 4 show that the biomass models were obviously improved on the fit and validation statistics after including the sample plot as a random effect into MEM. Consequently, the MEM in the above six approaches selected is probably more suitable for this study. Many studies have shown that biotic factors (e.g., DBH, H) and abiotic factors (e.g., origin, site, and climate) affect the biomass prediction accuracy [7,9,21,23,27]. The random effect "plot" added into the MEM could relate to climatic factors or site factors. In other words, BM, AMS1, AMS2, DMS1, and DMS2 consider the influence of biotic factors only, the MEM also takes abiotic factors into account, making it the most efficient among the six biomass modeling approaches. In fact, the fixed effects parameters in MEM has larger standard errors among all six approaches (Table 2). Thus, when a subsample of biomass is available to predict the random effects, the MEM is more efficient than the other five biomass modeling approaches. However, if the subsample is available, the MEM without the random effects would obtain less efficient estimates. At this point, we would recommend the use of the AMS2 or DMS1.