The study area is in the Central Highlands of Vietnam (Fig. 1). The average annual temperature is 21.60C and the average annual rainfall is 2,213 mm (Hydro-meteorological Station in the Central Highlands). Located at an altitude of 400-800 m above sea level, the soil types in the study area include red-brown soil on basalt, gray soil, and red-yellow soil on granite.

The studied agroforestry model includes the native multi-purpose tree species L. glutinosa combined with cassava. L. glutinosa is combined in different proportions with cassava ranging from 50% to 80% of the agroforestry area and is managed in a 6-7 year rotation cycle. Planting density (Nplant) ranged from 500 to 1967 plants ha^-1. L. glutinosa was grown from seeds collected at local tropical evergreen broadleaf forests in the first rotation, and in the following cycles, regenerated stems from coppice were used. The average number of stems plant^-1 (Nstemplant) in the second and third cycles ranged from 1 to 3.4. The number of stems ha^-1 (Nstem) ranged from 500 to 5900. After the third cycle, L. glutinosa was replanted from seeds (Table 1).

Twenty-two 300-m² circular plots were sampled covering a full range of ages (A) (1-7 years), a range of planting density, and a range of rotation cycles 1-3 in the agroforestry model. For all trees in the sample plots D (cm), H (m), and the Nstemplant were recorded. In each sample plot, one selected tree having the same D as the quadratic mean diameter (Dg), was destructively sampled to obtain a tree dry biomass/carbon dataset and its components. Age distributions of all sampled plots and destructively sampled trees are shown in Fig. 2. Before felling the sample tree, D, H, and crown diameter (CD, m) in two cardinal directions, North-South and East-West, and the sampled tree’s age (A, year) were recorded. The crown area of the trees sampled (CA, m²) was calculated by the formula

, where CD (the average crown diameter) was obtained from two cardinal direction measurements. After felling, tree heights were remeasured and fresh weights of the stem, branches, leaves, and bark were recorded. The four sampled tree components including stem, bark, branches and leaves were separated and weighed for fresh biomass in the field using a scale with a precision of 0.01 kg. Stem wood and bark samples were obtained from the base, middle, and top sections of the sampled tree. Sampled branches included small and large branches, and the sampled leaves on these branches. For four sampled tree components, a total of 88 samples were collected for the analysis; 100-300 g of each sample was weighed for fresh biomass on site using an electronic scale with a precision of 0.01 g.

The sample materials were dried at 105°C until a constant weight was attained. This provided the average fresh-to-dry mass ratio for all tree components to calculate the respective dry biomass – stem (Bst, kg tree^-1), branches (Bbr, kg tree^-1), leaves (Ble, kg tree^-1), and bark (Bba, kg tree^-1). Total tree aboveground biomass was computed as the sum of component biomass i.e., AGB = Bst + Bbr + Ble + Bba (kg tree^-1). The samples were analyzed after drying and the percentage of carbon was estimated by the Walkley & Black method (1934) and the carbon content of tree stem (Cst, kg tree^-1), branches (Cbr, kg tree^-1), leaves (Cle, kg tree^-1), bark (Cba, kg tree^-1), and the total tree aboveground carbon (AGC = Cst + Cbr + Cle + Cba, kg tree^-1) was calculated. Table 2 shows the summary statistics for each tree predictor and the destructively sampled tree response variables.

In this study, we compared two methods to fit systems of equations for AGB/AGC and their components: weighted nonlinear least-squares (WNLS) and weighted nonlinear seemingly unrelated regression (WNSUR) fit by the generalized least squares method.

— Relationship among tree biomass-carbon components and selection of predictors: The factor analysis method was performed to examine the relationship among the tree biomass and carbon components and select their predictors from D, H, A, and CA (Kim & Mueller, 1978; DeCoster, 1998). This method supports both principal components and classical factor analysis, which produces a linear combination of multiple quantitative variables and explains the largest percentage variation among those variables. The values of the variables were standardized by subtracting their means and dividing by their standard deviations.

— Model calibration – Independent fit: Fig. 3 shows that the relationship between tree biomass-carbon components, AGB and AGC vs. D conforms to the power law that was used in the study. Preliminary analysis showed that power models fit by the nonlinear method produced higher reliability than the log-linear model. This is consistent with previous studies (Huy et al., 2016a, 2016b, 2016c). Therefore, we used the nonlinear method to fit biomass-carbon modeling systems, and the heteroscedasticity in residuals was accounted for by using appropriate weighting (Davidian & Giltinan, 1995; Picard et al., 2012; Huy et al., 2016a, 2016b, 2016c, 2019). The WNLS models were fitted using nls function in the statistical software R (R Core Team, 2021). The model forms used were as follows:

where Y_i is the Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba or AGB/AGC in kg for the i^th sampled tree; a and b are the parameters of the model; X_i is the predictor(s) selected by factor analysis such as a combination of D (cm), H (m), A (year), D²H for the i^th sampled tree; and ε_i is the random error associated with the i^th sampled tree. The weighting variable is 1/D^δ or 1/(D²H)^δ and ^δ is selected in a range of ±2 (Picard et al., 2012).

— Model calibration – Simultaneous estimation: Independently fitted component models do not ensure that the AGB/AGC calculated as the sum of predicted tree component models is the same as the AGB/AGC estimated from independently developed AGB/AGC models (Sanquetta et al., 2015; Affleck & Dieguez-Aranda, 2016; Poudel & Temesgen, 2016; Gonzalez-Benecke et al., 2018; Huy et al., 2019; Trautenmüller et al., 2021). The SUR can solve that limitation by allowing simultaneous estimation of the tree biomass or carbon component and AGB/AGC (Parresol, 2001). Additionally, different weighting factors can be used for each equation to account for heteroscedasticity. The WNSUR models were fitted using SAS procedure Proc Model (SAS Inst., 2014) with the generalized least squares method. The modeling systems for the component and total tree biomass and carbon had the following general forms:

where Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC are biomass/carbon of stem, branches, leaves, bark, and total in kg, respectively; a_i and b_i are parameters of the power model i(i = 1, 2, 3, 4 for the stem, branches, leaves, and bark respectively); Xij is the predictor variables (D, H, D²HWWD) for the i^th equation and the j^th predictor; and ε_i is the residuals for the i^th equation (i = 1, 2, 3, 4, 5). The weighting variable is 1/D^δ or 1/(D²H)^δ and δ is selected in a range of ±2 (Picard et al., 2012).

Leave-one-out cross validation (LOOCV) is a special case of K-fold cross-validation, where K is set to the number of samples in the dataset. LOOCV has the maximum computational cost. It requires that one model be created and evaluated for each sample in the dataset (Cheng et al., 2017). LOOCV was applied for the cross-validation of modeling systems. The dataset was split into two parts with n-1 samples for model development and one for cross-validation. The process was repeated for all the observations, and the error statistics were computed at each iteration and averaged over the number of iterations.

The Akaike Information Criterion (AIC) (Akaike, 1973) compares and selects the best model. The adjusted R² describes the variability in the dependent variable explained by the set of predictors accounting for the number of predictors used in the model, whereas the diagnostic plots are used to examine residuals for any departure from model assumptions. The lower the AIC and the closer the adj. R² is to 1, the better the model. At the same time, cross-validation errors based on each sample such as bias (%), root mean square error (RMSE, kg tree^-1), and mean absolute percent error (MAPE, %) were used and averaged over the n realizations of LOOCV. The models with the smallest values of cross-validation errors were preferred.

where n is the number of realizations (number of samples); and y_i and ŷ_i are observed and predicted Bst/Cst, Bbr/Cbr, Ble/Cle, Bba/Cba and AGB/AGC for the i^th realization, respectively.

Bias and MAPE statistics estimated using Eqs. (8) and (10) have been the most widely used by many authors (e.g., Chave et al., 2005; Basuki et al., 2009; Huy et al., 2019, 2022). In some cases, the errors based on Eqs. (8) and (10) tend to be biased on negative errors, y_i < ŷ_i than in positive errors; therefore, the average systematic error (ASE) and the mean percent standard error (MPSE) were used (Zeng et al., 2017). The only difference among these metrics is the denominator; that is, instead of the observed value, the predicted value of the response is used in the denominator (Zeng et al., 2017; Huy et al., 2022).

Factor analysis was carried out with two groups of eight variables each. Group 1 included tree biomass components (Bst, Bbr, Ble, and Bba) and tree-level predictors (D, H, A, and CA) and group 2 included tree carbon components (Cst, Cbr, Cle and Cba) and tree predictors (D, H, A and CA). Factor analysis provided a small number of factors that account for most of the variability in the eight variables and supported an examination of the relationship between biomass and carbon components of trees. In both groups, two factors were extracted since two factors had eigenvalues greater than or equal to 1.0. Together, they accounted for 89.24% (group 1) and 89.30% (group 2) of the variability in the original data.

According to the results of the factor analysis, there were strong relationships between tree biomass and carbon components and tree predictors (D, H, A). They were placed close together in Fig. 4. They had a high coefficient (over 0.83) while CA had the lowest coefficient (below 0.30), which was shown in the first common factor equations (Factor 1) below. CA was excluded from further analysis because of its negligible influence on the biomass and carbon of the tree components.

The first common factor equations for tree biomass/carbon components are:

The results for selecting independent models using the WNLS method for tree biomass and carbon components and AGB/AGC vs. different combinations of variables D, H, and D²H by LOOCV are presented in Tables 3 and 4. We also observed that tree age (A) was marginally significant for all components and total (p > 0.05). The selected models from the independent fitting of component and total AGB and AGC were as follows:

The plots of observed vs. fitted values and weighted residuals vs. the fitted values (Fig. 5) depict a good model fit for independent fitting of component and total biomass. The results also showed that the model selection retained identical predictors or a combination of predictors when independent tree biomass and carbon components models were fitted (Table 3 and Table 4).

The results of calculating the average carbon fraction (CF = Carbon/Biomass weight) for tree biomass components and total and the variation within the 95% confidence range are shown in the box plot in Fig. 6.

The tree age was excluded from further analysis based on the results of independently fit models. Sixteen WNSUR equation systems were developed using D, H, and their combinations as predictors to predict AGB, AGC, and their components (Tables S1 and S2 [suppl]). Unlike the independent model fit, the results showed that the model selection retained the different predictors or a combination of predictors when these models were fit simultaneously. Based on the errors statistics of LOOCV, we found two optimal combinations of tree predictors for simultaneous modeling systems for tree biomass and carbon components and the totals are as follows:

The parameters of the two selected simultaneous modeling systems that were fit by the WNSUR method using the entire dataset are presented in Table 5. In addition to the simplicity of application, since measuring tree H is costly, two modeling systems for simultaneously predicting tree biomass and carbon were created with only D (Table 6).

Both WNLS and WNSUR addressed the issue of heteroscedasticity of the residuals. The WNLS residuals fluctuated uniformly along with the fitted values of the models, whereas the WNSUR model reduced the variability of the residuals significantly compared to the WNLS model (Fig. 7). The result of comparison of LOOCV statistics between two methods of independent WNLS and simultaneous WNSUR to fit selected AGB and AGC models is shown in Table 7.

Tree age (A) was insignificant in tree parts and total biomass and carbon models (Tables 3 and 4). Therefore, variables D and H (Dutca et al., 2018) may be sufficient to describe the variability in biomass and carbon in Litsea trees planted. Wood density (WD) is commonly used in mixed-species biomass models (e.g., Chave et al., 2005; Basuki et al., 2009; Chave et al., 2014; Huy et al., 2019). However, as WD reflects the differences in biomass accumulation among species, it may not be a critical predictor for the species-specific biomass models.

The biomass and carbon in branches (Bbr/Cbr), leaves (Ble/Cle), and bark (Bba/Cba) of the tree have greater variation than the biomass and carbon in the stem (Bst/Cst) and total (AGB/AGC). When the models for estimating AGC and its carbon components were developed separately, the errors in the selected models for stem carbon (Cst) and total (AGC) were the smallest with MAPE = 19.6-20.1% (Table 4). Meanwhile, the errors of the selected branches and leaves carbon models were much higher, with MAPE = 29.8-33.9% (Table 4). On the other hand, with the chosen simultaneous modeling system, the errors in the models of stem carbon (Cst), and total (AGC) were the smallest, with MAPE = 15.0-16.2% (Table S2), and the errors in the models of branches, leaves, and bark carbon were much higher, MAPE = 25.1-34.8% (Table S2). This may be because of the high variability in the stand density (Table 1) of the agroforestry model studied here. As the density is low, branches and foliage grow stronger and wider, so there are large variations in Bbr/Cbr and Ble/Cle, making these fitted models have larger errors than other models like Bst/Cst and AGB/AGC.

Simultaneous modeling systems fit by WNSUR reduced MAPE by 1.3% - 3.9% of AGB and AGC estimates compared to independent models fit by WNLS (Table 7). This result is consistent with Huy et al. (2019) and Trautenmüller et al. (2021). As expected, the sum of the biomass of the components predicted with independently fit component models differed from the estimates obtained from the AGB/AGC models. Using non-additive independent biomass models of tree components and total AGB produces biologically inconsistent estimates that scale up to a large area, which affects the final total biomass estimated. Simultaneous fitting leads to higher efficiency than independent fitting because the variance and covariance information of the tree biomass- carbon components and total are included in the model (Parresol, 2001; Poudel & Temesgen, 2016; Trautenmüller et al., 2021).

Given the above observations along with Fig. 7, a procedure using a simultaneous and weighted modeling system fit by WNSUR is recommended to generate native tree biomass-carbon models and their components. This recommendation is consistent with the comments of Huy et al. (2019) and Trautenmüller et al. (2021).

The agroforestry model aims to create economic, social, and ecological efficiency. The design and implementation of an agroforestry model that balances economic and environmental factors is always the primary concern, but is not easy to study. Using the selected modeling systems fit by WNSUR to simultaneously predict the biomass-carbon sequestration of L. glutinosa including Bba, AGB and AGC associated with different densities of Nplant and Nstem to calculate total stem bark biomass ha^-1, total AGB ha^-1, total AGC ha^-1 and total CO² equivalent absorbed ha^-1 for each agroforestry model.

As a result, a popular agroforestry model with a density of 1633 Nplant and 1960 Nstem of L. glutinosa, uses 65% of the space and the remaining 35% of the space to grow cassava until the end of the 6-7 year rotation cycle of L. glutinosa, the most valuable component of L. glutinosa is the stem bark biomass, which reached 4.7 tons ha^-1 and the total AGB, total AGC, and total CO² equivalent accumulation reached 15.0 tons ha^-1 (2.5 tons ha^-1 year-1), 7.1 tons ha^-1 (1.2 tons ha^-1 year-1) and 26.0 tons ha^-1 (4.3 tons ha^-1 year-1), respectively. In Austrian mountain agroecosystems, carbon is sequestered in perennial biomass by up to 3.1 tons ha^-1 year-1 (Bertsch-Hoermann et al., 2021), the potential of agroforestry systems in tropical India to accumulate carbon is estimated at 0.3-15.2 tons ha^-1 year-1 (Dhyani et al., 2020). In comparison, the carbon accumulation in the agroforestry model studied here is lower than that of the agroforestry systems in Europe and averages in the Asian region.

Most of the published models have been used to predict tree biomass and apply the default CF of IPCC (2006) for plants is 0.47 to carbon conversion. The CF of the stem bark observed in this study was slightly lower (CF = 0.46), but the CF of the total and other components was higher (CF = 0.48 for branches and total; CF=0.49 for leaves) (Fig. 6). Comparison with this result showed differences in CF among tree parts, therefore, using the same default CF = 0.47 for all parts of the tree could result in a significant error in the conversion. Therefore, the selected modeling system in this study will provide greater confidence in predicting the tree component carbon sequestration than using the tree biomass models and the IPCC’s CF.