Additive tree biomass equations for Betula platyphylla Suk. plantations in Northeast China

Abstract

Key message

A new system of additive tree biomass equations was developed for juvenile white birch plantations based on tree diameter at breast height (DBH) and tree height (HT). Compared with previous equations developed for natural white birch forests, the new system included one more biomass component and provided more accurate predictions.

Context

Accurate estimates of tree component and total biomass are necessary for evaluating alternative forest management strategies for biomass feedstock, carbon sequestration, and products. Previous biomass equations developed for white birch trees in natural stands provided substantially biased predictions for white birch plantations.

Aims

A new system of additive tree biomass equations was developed for juvenile white birch plantations in the northeastern China.

Methods

With destructive biomass sampling data from 501 trees sampled from white birch provenance and family trails at ages 7, 9, 10, and 13 in three provinces, a system of nonlinear additive tree biomass equations based on DBH and tree height was developed using the nonlinear seemingly unrelated regressions (NSUR) approach.

Results

Compared with previously published equations developed for natural white birch forests, the new system provided more accurate predictions of white birch tree component and aboveground and total biomass, especially of branch, foliage, and root biomass.

Conclusion

The new system extended the applicability of biomass equations to white birch plantations in the northeastern China.

1 Introduction

Accurate estimates of total tree biomass and tree component biomass are necessary for evaluating alternative forest management strategies for biomass feedstock, carbon sequestration, and products (Castedo-Dorado et al. 2012; Zhao et al. 2015). Individual tree biomass equations are needed to estimate biomass in foliage, branches, stems, and other tree components (Affleck and Dieguez-Aranda 2016). A desirable feature of a system of equations is that the predictions for tree components sum to the prediction for total tree biomass (Parresol 1999). Additive biomass equations have been used to estimate tree component biomass for some conifer tree species (Bi et al. 2010; Zhao et al. 2015) and broad-leaved tree species (Bi et al. 2004, 2015; Zheng et al. 2015). When estimating a system of additive biomass equations, taking into account the inherent correlation among the biomass components has greater statistical efficiency (Parresol 1999, 2001).

White birch (Betula platyphylla Suk.) is widely distributed in the temperate broad-leaved forests of East Asia (Kuang et al. 1979). It can grow in pure stands or in mixtures with Acer, Larix, Picea, Tilia, and another species of Betula. White birch is an important reforestation tree species and plays an important ecological role as a pioneer species of secondary broad-leaved forests of northern China (Zhang et al. 2012; Wang et al. 2015). It is also one of the most important commercial tree species for paper, furniture, and plywood production in China (Li et al. 1995). White birch is the most extensively distributed broad-leaved tree in northern and southwestern forest areas of China because of its rapid growth and adaptation to a wide variety of sites (Li et al. 1995). White birch represents one of the largest standing volumes among the hardwood species in the area (Fang et al. 2011). In Daxing′an Mountains, nearly one third of the forests are dominated by white birch (Zeng et al. 2003). The natural birch forest area in China has been gradually dwindling since the 1980s because of pests, disease, and environmental damage. In response, region-wide provenance and family trials were established starting in 1991with the main goal of identifying fast-growing seed sources to enhance white birch plantation productivity (Zhu et al. 2001).

Additive biomass equation systems have been developed for white birch trees in natural stands (Dong et al. 2013, 2015). However, these equations may not be appropriate for estimating tree biomass grown in white birch plantations, which are extensively distributed in the region. Tree biomass allocation among different components may also differ between plantation forests and natural forests, because management treatments such as thinning, pruning, fertilization, and planting specific genotypes selected for high yield can affect the accumulation and partitioning of biomass among tree components (Albaugh et al. 2009).

In this study, we have done extensive destructive sampling to collect biomass data on individual trees and developed a new system of additive biomass equations for white birch trees growing in plantations. The new system fulfills the additivity property, takes into account inherent correlations among tree biomass equations, and overcomes the heteroscedasticity problem. The predictive performance of the new system was compared with the biomass equations developed by Dong et al. (2013, 2015) for white birch trees in natural forests.

2 Materials and methods

2.1 Study sites and experiment design

A total of 501 trees were sampled from white birth plantations on three sites (HLJ, JL, and LN) in Northeastern China. Detail information about geographic coordinates, climate variables, and experimental sites can be found in Table 1 and Fig. 1. Three plantations are on the HLJ site located in Maoer mountain, Heilongjiang province. The first one was a 13-year-old provenance trial consisting of 16 provenances: 12, 3, and 1 from northeastern-, central-, and northwestern China, respectively. The trees from the different provenances were planted in a randomized complete block design of 4 blocks and 6 tree-row plots. The planting density was 1.5 m × 1.5 m. The second plantation was a 9-year-old half-sib family trial with 60 half-sib families that were planted in a randomized complete block design of 4 blocks and 10 tree- 2 row plots, respectively.Their female parents were from 10 provenances (DFH, FL, HR, HN, LS, MES, QY, WYL, XBH, and XJ). The planting density was 2 m × 2 m. The third plantation was a 7-year-old full-sib family trial with 20 full-sib families that were planted in a randomized complete block design of 4 blocks and 10 tree- 2 row plots, respectively. Their parents were from MES, XBH, and Finland (not shown in Fig. 1).The planting density was 2 m × 2 m. The JL site, in Jilin province, consisted of a single 10-year-old half-sib family trial plantation with 9 half-sib families that were planted in a randomized complete block design of 4 blocks and 10 tree- 2 rows plots. Their female parents were from 5 provenances (DFH, HR, MERS, XBH, LS). The planting density was 2 m × 2 m. The LN site was located in Caohekou, Liaoning province. The plantation was a single 13-year-old provenance trial plantation with 16 provenances (12 from northeastern, 3 from central-, and 1 from northwestern China) planted in a randomized complete block design of 4 blocks and 6 tree-row plots. The planting density was 1.5 m × 0.75 m.

Table 1 Geographic coordinates, temperature, and precipitation at the sites

Fig. 1

Geographic locations of white birch provenances and experimental sites

2.2 Sampling and biomass measurements

On the HLJ site, 106 trees were destructively sampled from stands of the provenance trials in August 2011, 175 trees were sampled from half-sib families, and 60 trees sampled from full-sib families in August 2012. On the JL site, 26 trees were sampled from stands of half-sib families in August 2012; on the LN site, 134 trees were sampled from stands of the provenance trials in August 2011. Two sample trees per provenance per block, one sample tree per half-sib family per block, and one sample tree per full-sib family per block were randomly selected. Destructive biomass sampling was not conducted on some stands due to poor survivals.

Biomass measurements were conducted following the procedures of Wang et al. (2015). Aboveground biomass was divided into stem wood, stem bark, branch, and leaf components. The stem (including stem wood and stem bark) was divided into 1 m long sections. The fresh mass of each stem section was recorded. A 5-cm disc subsample was cut from the middle of each section. The fresh mass of bark and wood from the 5 cm disk were recorded for determining the ratio of fresh mass of bark and wood of the stem. The total fresh mass of all tree branches (including leaves) of each tree were measured and recorded. The branches of each tree were divided into three crown positions (lower, middle, upper) and one third of the branches of each tree sub-sampled for determining the ratio of leaf to branch. Then, subsamples were randomly taken from mixed branches (about 500 g/tree) and leaves (about 500 g/tree) of canopy crown positions of each tree. The entire roots were dug out manually and carefully with shovels. The total fresh mass of roots were recorded. Subsamples of fresh roots were randomly sampled (about 700 g/tree). All samples were oven-dried at 70 °C until constant mass and the ratio of dry to fresh mass were calculated. Dry biomass of each tree component was calculated by multiplying its fresh mass by the respective dry/fresh mass ratio.

Summary statistics for DBH, total height, and component biomass of all sample trees are shown in Table 2. The relationships of stem wood, stem bark, branch, leaf, and root biomass with tree DBH and height are shown in Fig. 2.

Table 2 Summary statistics of tree height (HT), diameter at breast height (DBH), and biomass for the sampled trees (n = 501)

Fig. 2

Relationship between stem wood, stem bark, branch, leaf, aboveground, root, and total biomass of the sampled trees and tree diameter at breast height (DBH) and height

2.3 Model description

A system of seven equations with additive error terms, cross-equation constraints on the structure parameters, and cross-equation correlation for five three tree biomass components (stem wood, bark, branch, leaf, and root), tree aboveground, total tree biomass with additivity was used:

$$ {\displaystyle \begin{array}{l}\ \mathrm{Stem}:\kern1.5em {Y}_1={f}_1\left({\mathbf{X}}_1,{\boldsymbol{\upbeta}}_1\right)+{\boldsymbol{\upvarepsilon}}_1\\ {}\ \mathrm{Bark}:\kern1.5em {Y}_2={f}_2\left({\mathbf{X}}_2,{\boldsymbol{\upbeta}}_2\right)+{\boldsymbol{\upvarepsilon}}_2\\ {}\ \mathrm{Branch}:\kern0.5em {Y}_3={f}_3\left({\mathbf{X}}_3,{\boldsymbol{\upbeta}}_3\right)+{\boldsymbol{\upvarepsilon}}_3\\ {}\ \mathrm{Leaf}:\kern1.75em {Y}_4={f}_4\left({\mathbf{X}}_4,{\boldsymbol{\upbeta}}_4\right)+{\boldsymbol{\upvarepsilon}}_4\\ {}\ \mathrm{Root}:\kern1.75em {Y}_5={f}_5\left({\mathbf{X}}_5,{\boldsymbol{\upbeta}}_5\right)+{\boldsymbol{\upvarepsilon}}_5\\ {}\ \mathrm{Above}:\kern1em {Y}_6={f}_1\left({\mathbf{X}}_1,{\boldsymbol{\upbeta}}_1\right)+{f}_2\left({\mathbf{X}}_2,{\boldsymbol{\upbeta}}_2\right)+{f}_3\left({\mathbf{X}}_3,{\boldsymbol{\upbeta}}_3\right)+{f}_4\left({\mathbf{X}}_4,{\boldsymbol{\upbeta}}_4\right)+{\boldsymbol{\upvarepsilon}}_6\\ {}\ \mathrm{Total}:\kern1.5em {Y}_7={f}_1\left({\mathbf{X}}_1,{\boldsymbol{\upbeta}}_1\right)+{f}_2\left({\mathbf{X}}_2,{\boldsymbol{\upbeta}}_2\right)+{f}_3\left({\mathbf{X}}_3,{\boldsymbol{\upbeta}}_3\right)+{f}_4\left({\mathbf{X}}_4,{\boldsymbol{\upbeta}}_4\right)+{f}_5\left({\mathbf{X}}_5,{\boldsymbol{\upbeta}}_5\right)+{\boldsymbol{\upvarepsilon}}_7\end{array}} $$

(1)

where, Y_i represents the vector of stem wood, stem bark, branch, foliage, root, above ground biomass, and total tree biomass, respectively; f_l(X_l, β_l) is a nonlinear function for tree biomass component (l = 1, …, 5) for stem wood, bark, branch, leaf, and roots, respectively); ε_i is the n × 1 vector of residuals for the i^th equation (i = 1, …, 7), and n is the number of observations (trees).

Observations from different trees are generally taken to be independent. Assume ${\boldsymbol{\upvarepsilon}}_i\sim N\left(\mathbf{0},{\sigma}_i^2{\boldsymbol{\uppsi}}_i\right)$, where ψ_i is a (n × n) diagonal matrix. Heteroscedasticity in the i^th equation is decribed by ψ_i of which the diagonal elements are not all identical.

Let${\mathbf{V}}_i={\sigma}_i\sqrt{{\boldsymbol{\uppsi}}_i^{-1}}$ and $V={\left[\begin{array}{cccc}{\mathbf{V}}_1& \mathbf{0}& \dots & \mathbf{0}\\ {}\mathbf{0}& {\mathbf{V}}_2& \dots & \mathbf{0}\\ {}\dots & \dots & \dots & \dots \\ {}\mathbf{0}& \mathbf{0}& \dots & {\mathbf{V}}_7\end{array}\right]}_{\left(7n\times 7n\right)}$.

Inherent correlations among biomass components measured on the same tree are described by the correlations matrix among biomass equations:

$$ \boldsymbol{\uprho} =\left[\begin{array}{cccc}1& {\rho}_{12}& \dots & {\rho}_{17}\\ {}{\rho}_{21}& 1& \dots & {\rho}_{27}\\ {}\dots & \dots & \dots & \dots \\ {}{\rho}_{71}& {\rho}_{72}& \dots & 1\end{array}\right]. $$

Let C = ρ ⊗ I_n, where ⊗ is the Kronecker product; R = VCV. Now, the distribution for the residual terms ε = (ε₁, ε₂, …, ε₇)^′ is assumed to follow ε ∼ N(0, R).

Tree component biomass can be modeled as a power function of tree dimensions as:

$$ {f}_l\left({\mathbf{X}}_l,{\boldsymbol{\upbeta}}_l\right)={\beta}_{l0}{DBH}^{\beta_{l1}}{HT}^{\beta_{l2}} $$

(2)

where X_j is tree DBH, total tree height (HT), and β_l = (β_l0, β_l1, β_l2). Each component equation can contain its own independent variables.

The system of equations was fitted using the four-step fitting method (Zhao et al. 2015) with nonlinear seemingly unrelated regression (NSUR) and using the SAS/ETS® MODEL Procedure (SAS Institute Inc. 2011). This approach guarantees additivity in biomass equations, accounts for the inherent correlation among the biomass equations, and addresses heteroscedasticity by having a unique weighting function for each equation.

2.4 Model assessment and evaluation

Four fit statistics were obtained for each equation and used to evaluate the goodness of fit for the biomass prediction system: mean residual (E), mean of the absolute value of residuals (MABE), root mean square error (RMSE), and the coefficient of determination (R²). Mathematical expressions of these criteria are:

$$ {E}_i=\frac{\sum \limits_{j=1}^n\left({Y}_{ij}-{\hat{Y}}_{ij}\right)}{n} $$

(3)

$$ {\mathrm{MABE}}_i=\frac{\sum \limits_{j=1}^n\mid {Y}_{ij}-{\hat{Y}}_{ij}\mid }{n} $$

(4)

$$ {\mathrm{RMSE}}_i=\sqrt{\frac{\sum \limits_{j=1}^n{\left({Y}_{ij}-{\hat{Y}}_{ij}\right)}^2}{n}} $$

(5)

$$ {R}_i^2=1-\frac{\sum \limits_{j=1}^n{\left({Y}_{ij}-{\widehat{Y}}_{ij}\right)}^2}{\sum \limits_{j=1}^n{\left({Y}_{ij}-{\overline{Y}}_i\right)}^2} $$

(6)

where ${Y}_{ij}\ \mathrm{and}\ {\hat{Y}}_{ij}$ are the jth observed and predicted values of biomass for the ith component, and ${\overline{Y}}_i$ is the mean of n observed values for the same component.

In this study, the biomass equation system was fitted to the entire data set (N = 501 trees). Model validation was accomplished by the leave-one-out (LOO) cross-validation technique, in which the model system was fitted using all-but-one tree (leaving one tree out), and then, the fitted model system was used to predict the values of all component and total tree biomass for that left-out tree. The summary statistics were calculated using the same formulas (3, 4, 5, and 6).

The predictive performance of the new model system was compared graphically with the previously published equations for white birch trees in natural forests (Dong et al. 2013, 2015). The biomass estimates from the equations of Dong et al. (2013, 2015) were also assessed with the same criteria as the new biomass equations, using formulas (3, 4, 5, and 6).

Data availability

The datasets generated and/or analyzed during the current study are available from the corresponding author.

3 Results

3.1 Biomass equations based on DBH and HT

The fitted biomass equation system is shown below. The parameters acting as powers of DBH and HT were highly significant in each biomass equation (Table 3).

$$ {\displaystyle \begin{array}{l}\mathrm{Wood}:\kern0.75em {\hat{Y}}_1={\widehat{\beta}}_{10}{DBH}^{{\widehat{\beta}}_{11}}{HT}^{{\widehat{\beta}}_{12}}\\ {}\mathrm{Bark}:\kern1.25em {\hat{Y}}_2={\hat{\beta}}_{20}{DBH}^{{\hat{\beta}}_{21}}{HT}^{{\hat{\beta}}_{22}}\\ {}\mathrm{Branch}:{\hat{Y}}_3={\hat{\beta}}_{30}{DBH}^{{\hat{\beta}}_{31}}{HT}^{{\hat{\beta}}_{32}}\\ {}\mathrm{Leaf}:\kern1.25em {\hat{Y}}_4={\hat{\beta}}_{40}{DBH}^{{\hat{\beta}}_{41}}{HT}^{{\hat{\beta}}_{42}}\\ {}\mathrm{Roots}:\kern0.75em {\hat{Y}}_5={\hat{\beta}}_{50}{DBH}^{{\hat{\beta}}_{51}}{HT}^{{\hat{\beta}}_{52}}\\ {}\mathrm{Above}:\kern0.5em {\hat{Y}}_6={\hat{Y}}_1+{\hat{Y}}_2+{\hat{Y}}_3+{\hat{Y}}_4\\ {}\mathrm{Total}:\kern1em {\hat{Y}}_7={\hat{Y}}_1+{\hat{Y}}_2+{\hat{Y}}_3+{\hat{Y}}_4+{\hat{Y}}_5\end{array}} $$

(7)

Table 3 Parameter estimates and their asymptotic standard error and p values for the additive biomass equation system (7)

The coefficients of DBH were positive in each biomass component, while the coefficients associated with HT were positive in stem wood and stem bark, and negative in branch, leaf, and root components (Table 3). The positive coefficients of DBH and HT suggested their positive relationship with stem wood and bark biomass. This implies that, for the same DBH, tree stem wood and bark biomass increased with increasing tree height. The positive coefficient of DBH and the negative coefficient of HT in branch, leaf, and root components implied that their biomass increased with increasing DBH, but for the same DBH, branch, leaf, and root biomass decreased with increasing tree height.

The system of biomass equations fitted well for stem wood, bark, aboveground biomass and total tree biomass, but fitted poorly for branch, leaf, and root biomass (Table 4). The weighting functions involved either DBH or HT or both, depending on the component biomass equation, suggesting the necessity of different weighting functions for each biomass equation (Table 4). The residual variances were stabilized with weight functions (Fig. 3). There were high correlations among most of the biomass equations, as shown in the following correlation matrix:

$$ {\displaystyle \begin{array}{l}\kern2.75em \mathrm{Stem}\kern0.5em \mathrm{Bark}\kern0.5em \mathrm{Branch}\kern0.5em \mathrm{Leaf}\kern0.5em \mathrm{Above}\kern0.5em \mathrm{Root}\kern0.5em \mathrm{Total}\kern0.75em \\ {}\begin{array}{c}\mathrm{Stem}\\ {}\mathrm{Bark}\\ {}\mathrm{Branch}\\ {}\mathrm{Leaf}\\ {}\mathrm{Above}\\ {}\mathrm{Root}\\ {}\mathrm{Total}\end{array}\ \left(\begin{array}{ccccccc}1& 0.378& 0.198& 0.066& 0.816& 0.096& 0.607\\ {}& 1& 0.025& 0.128& 0.474& -0.169& 0.221\\ {}& & 1& 0.404& 0.637& 0.370& 0.662\\ {}& & & 1& 0.459& 0.070& 0.361\\ {}& & & & 1& 0.207& 0.802\\ {}& & & & & 1& 0.739\\ {}& & & & & & 1\end{array}\right)\end{array}} $$

Table 4 Weight functions and fit statistics for each biomass component in the additive biomass equation system (7)

Fig. 3

Pearson residual plots for each biomass component in the equation system fitted using NSUR method and different weight functions for each system equation with its own weight function

3.2 Biomass model validation and comparison

The leave-one-out cross-validation statistics indicated that the additive biomass equation system (7) slightly underestimated stem bark (0.3%), branch (1.1%), leaf (0.9%), and total biomass (1.4%), underestimated root (6.5%), and slightly overestimated stem wood (1.3%) and total aboveground biomass (0.5%) (Table 5).

Table 5 Leave-one-out (LOO) cross-validation results for each biomass component in the additive biomass equation system (7)

Comparing the actual component biomass with the biomass predicted by the newly developed equation system as well as with the equations developed by Dong et al. (2013, 2015) showed that the new equation system predicted all biomass components better than the previous equations (Table 6,Fig. 4). The new equations substantially improved the estimation of component and total biomass. The equations of Dong et al. (2013, 2015) substantially overestimated stem wood, aboveground, and total biomass.

Table 6 Statistics for predicting the component biomass of white birch trees from the newly developed equations (NM), the equations developed by Dong et al. (2013) (D13) and Dong et al. (2015) (D15)

Fig. 4

Comparison of stem wood, leaves, branch, aboveground, root, and total biomass predictions for the new equation system (NM) and the equation systems of Dong et al. (2013) (D13) and Dong et al. (2015) (D15)

4 Discussion and conclusions

We developed a set of biomass equations for trees growing in white birch plantations. Compared with previous equations developed for natural white birch forests, the new system of equations included one more biomass component—stem bark component—and provided more accurate predictions of branch, foliage, and root biomass. The new system of equations was particularly accurate for predicting stem wood, total aboveground, and total tree biomass in white birch plantations. Biomass equations developed by Dong et al. (2013, 2015) for trees in natural white birch forests are not suitable for estimating component biomass and total biomass of trees in white birch plantations (Table 6 and Fig. 4). Natural forests and plantations are different in many ways including stand density and the level of competition from adjacent trees so the pattern of biomass partitioning should also differ (Satoo and Madgwick 1982).

The new equations were developed with a large variety of genotypes across a wide geographical area and a wide range of growth rates. Planting materials came from 16 provinces from northeastern to northwestern China and included the progenies of 60 half-sib families and 20 full-sib families obtained by crossing, even a parent of full-sib families was from Finland. The new equations were developed with tree biomass data from juvenile white birch plantations (7, 9, 10, and 13 years old). The previous equations were developed with tree data from natural forests of unknown ages (Dong et al. 2013, 2015). The average DBH of sampled trees used in previous studies was in the DBH range used in the current study. Tree DBHs in this study ranged from 3.4 to 14.2 cm, while tree DBHs in previous studies ranged from 8.0 to 33.1 cm (Dong et al. 2013) or from 5.4 to 33.1 cm (Dong et al. 2015). Models built to predict growth for smaller trees and juvenile stands are not as common as those for older stands (Vaughn 2007). Due to the lack of biomass equations for small diameter trees, the carbon stored in this component of temperate forests has been often ignored (Daryaei and Sohrabi 2016). Tree size has a strong effect on the biomass partitioning patterns (Mensah et al. 2016). Modeling juvenile growth of trees is important for a better understanding of the whole process of stand development and helping to schedule appropriate silvicultural treatments for young stands (Zhao et al. 2015).

The model structure as in Eq. 1 was used to ensure the additivity property of nonlinear biomass models. Most importantly, the system was fitted with the weighted NSUR and using a large dataset. Gerbing and Anderson (1985) investigated the effects of several variables and found that sample size had the largest effect on the variance in parameter estimates. In the previous studies, Dong et al. (2013) used 66 white birch trees and Dong et al. (2015) used 98 white birch trees. Biomass dataset used in our current study consisted of 501 white birch trees. Larger sample size can reduce parameter estimation uncertainty. Heteroscedasticity almost always exists in biomass models, as our results show. It can inflate the standard error of the estimate of the parameters. In this study, we addressed this problem by having a unique weighting function for each biomass equation and would achieve minimum variance estimates and reliable prediction intervals. The correlation matrix among biomass equations showed that the equations do have correlated errors. So the NSUR used in the current study to take the cross-equation error correlation in account would result in more precise estimate of the parameters than separate estimation of equations in the system.

The new equations substantially improved the estimation of branch, foliage, and root biomass. However, both the new and previous system of tree biomass equations for white birch which are based on DBH and HT, did not estimate branch, foliage, and root biomass well. Zhao et al. (2015) reported that more detailed crown size measurements including crown length, crown width, and diameter at the base of live crown could improve crown biomass prediction. In future work, additional crown information will be collected and integrated into the system of biomass equations to improve the prediction of branch, foliage, and root biomass.

The total biomass (or subtotal biomass) was obtained by adding component biomass together. Its variance should be a function of the component biomass function and cross-correlation parameters (Affleck and Dieguez-Aranda 2016). In our model system, however, the total (or subtotal) biomass equation was treated as like component equations, specifying one additional variance function and additional cross-correlation parameters. It should be noted that the residual variance of the total biomass estimated in the way we did could likely be biased. Although this bias is not a big deal for most applications, it is worth while to study further.