Modeling biomass allocation strategy of young planted Zelkova serrata trees in Taiwan with component ratio method and seemingly unrelated regressions

Trees accumulate biomass by sequestrating atmospheric carbon and allocate it to different tree components. A biomass component ratio is the ratio of biomass in a tree component to total tree biomass. Modeling the ratios for Zelkova serrata, an important native reforestation tree species in Taiwan, helps in understanding its biomass allocation strategy to design effective silvicultural treatments. In this study, we applied Component Ratio Method (CRM) to relate biomass component ratios of main stem, large branch, twig, and foliage to tree attributes of Z. serrata from a 9-year-old plantation. Nonlinear and linear CRM models were fitted with Seemingly Unrelated Regression to account for model correlations. Linear CRM models with dbh as the predictor had the best fit with model correlations as high as 80%. About 46% and 40% of total tree biomass was allocated to main stem and large branch, respectively. However, main stem biomass decreased by 1.9% with every 1-cm increase in dbh, but large branch biomass increased by 2.2% instead. Results suggest that dominant Z. serrata trees tend to branch and fork, while smaller trees invest in larger main stem. An early pruning treatment should focus on dominant trees to maintain crown ratio and ensure wood quality.

A unique contribution of trees to ecosystem services beneficial to human society is the accumulation of carbon in the form of biomass. As such, the Paris Agreement formally recognizes that forests play an important role in addressing the impact of climate change by sequestrating carbon from atmosphere 1 . As the world is moving towards decarbonization, many mitigation methods have been developed such as negative emissions technologies (NETs) 2 and radiative forcing geoengineering 3 . Among the different NETs, afforestation and reforestation approach 2 is immediately relevant to forest management. By planting reclaimed lands or degraded forests, standing trees accumulate biomass throughout their life cycles albeit at different rates. Forests could also contribute to other climate change mitigation methods such as biochar production 4 . Matovic suggested that about 4.8 Gt of carbon could be sequestrated if 10% of the world biomass was converted to biochar, and part of the biomass could be sourced from forest management activities 4 . Osman et al. comprehensively reviewed several decarbonization technologies and found that using plant biomass as fuel in the oxyfuel combustion route could promote bioenergy and carbon capture and storage (BECCS) system as an effective way to achieve decarbonization 5 . Thus, there are many pathways available for forests to assist climate change mitigation.
Tree biomass is not directly measurable and is usually estimated by different methods. The most common approach is first estimating tree volume from forest inventory data and converting it to tree biomass by biomass expansion factors 6,7 . An alternative approach is using allometric equations to estimate tree total biomass or biomass of each tree component from stem diameter at breast height (dbh) 6,8 . Modeling how a tree partitions its total biomass into different tree components is needed. It is because understanding distribution of wood biomass www.nature.com/scientificreports/ within a tree is important to estimate its wood utilization potential and economic value of tree components 9 . This also helps with designing appropriate silviculture treatments that promote biomass accumulation in certain tree components. Jenkins et al. proposed a Component Ratio Method (CRM) in which ratio of biomass of a tree component was nonlinearly associated with dbh 10 . A biomass component ratio is calculated as the proportion of biomass of a tree component (e.g., main stem) to total tree aboveground biomass 10 . As a result, ratios of all biomass components should sum up to one. CRM has the advantage of additivity such that predicted biomasses of all the components in a tree sum to the predicted total tree biomass 11,12 . A hybrid approach could be used to estimate tree component biomass: tree total biomass is first predicted from tree volume by a biomass expansion factor, and CRM is then used to consistently distribute this total to each tree component 6 . Thus, CRM could assist in designing silvicultural treatments that maximize wood utilization potential of a whole tree or certain tree components. While CRM is very useful, there are very few studies that model relationship between biomass component ratios and dbh 13,14 .
In Taiwan, forests occupy about 60% of the total land area of about 36,200 km 2 , which could be grouped into five general types: broadleaf forest, coniferous forest, mixed broadleaf-coniferous forest, mixed bamboo broadleaf forest, and mixed bamboo coniferous forest 15,16 . The 4th Taiwan National Forest Resource Inventory conducted between 2008 and 2012 reported that total carbon storage in Taiwan forests was 754.3 Mt CO 2 or 347.9 t CO 2 ha −115,17 . Of the total carbon storage, the largest proportion of carbon is stored in broadleaf forests (468.9 Mt CO 2 ) followed by coniferous forests (156.3 Mt CO 2 ) and mixed broadleaf-coniferous forests (103.5 Mt CO 2 ) 15 . However, on a per unit area basis, mixed broadleaf-coniferous forests store the most amount of carbon (604.7 t CO 2 ha −1 ) followed by coniferous forests (522.3 t CO 2 ha −1 ) and broadleaf forests (319.1 t CO 2 ha −1 ) 15 . Zelkova serrata (Thunb.) Makino is an important reforestation broadleaf tree species native to Taiwan. It distributes from 300 to 1000 m a.s.l. and has high economic value due to its desirable wood properties for construction, furniture, flooring, and sculpture 18 . As such, Z. serrata is one of the ten major tree species widely planted since the start of the Taiwan National Reforestation Program in 1996 19 . A total of 1785 ha of land has been reforested with the tree species between 1997 and 1999 19 . Because of its importance as a plantation tree species, there is a continuing interest in understanding its role in carbon sequestration. The 4th Taiwan National Forest Resource Inventory reported that Z. serrata plantation had an average stand volume of 193.2 m 3 ha −1 , which stored about 452.2 t CO 2 ha −1 -higher than the average storage rate of broadleaf forests 15,17 . A 9-year-old Z. serrata plantation could hold up to 267.9 t CO 2 ha −118 . Depending on stand age, annual stand-level carbon sequestration rate could be between 1.81 to 4.11 t CO 2 ha −1 year −118,20 . The biomass expansion factors for a 25-year-old and a 46-year-old Z. serrata plantations were estimated to be 1.328 Mg m -3 and 1.528 Mg m -3 , respectively 20 . Lastly, CO 2 fixation rates of the upper-leaf and lower-leaf of Z. serrata species were estimated to be 5.52 g m -2 s −1 and 2.38 g m -2 s −1 , respectively 21 .
While past studies have assessed carbon sequestration potential and biomass of Z. serrata on a stand-level, they have not explored the strategy adopted by the tree species in distributing its total biomass among its various components. As mentioned above, understanding this allocation strategy has ecological, economic, and management implications. Thus, to fill in the knowledge gap, the goals of this study were to apply CRM to model treelevel relationship between proportion of biomass in each tree component and dbh of Z. serrata, and to suggest potential silvicultural treatments that improve wood utilization potential of the tree species.

Materials and methods
Study site. The study site was established at a Z. serrata plantation in the Neimaopu Forest District of the National Taiwan University Experimental Forest. The plantation was established in 1997 with a planting density of 1500 trees/ha and an area of 1.3 ha. It was located at 23° 40′ N and 120° 50′ E at 800 m a.s.l. Mean annual precipitation in the area was 1853 mm between 1997 and 2004. Mean annual temperature was 21.5 °C with mean relative humidity of 81.4%. This study was carried out in 2005 when the planted Z. serrata trees were 9 years old. All trees in the 1.3 ha plantation were censused for dbh and tree height (ht). A total of 921 trees were measured. The sampled trees were grouped into five diameter classes of 5-cm width (Table 1). A total of 12 trees were randomly selected from the first four diameter classes for biomass study with two trees from the diameter class of ≤ 5 cm, four trees from the diameter class of 5-10 cm, and three trees from each of the diameter classes of 10-15 and 15-20 cm (Table 1). Since the last diameter class has only 3 sampled trees, no tree was selected from it. Thus, the 12 selected trees represented the range of tree attributes in the study site. Table 1. Number of census trees and sampled trees for biomass study in each diameter class. Mean diameter at breast height (dbh) and mean tree height of the trees sampled for biomass study with their respective coefficient of variation (%) in brackets are reported for each diameter class. www.nature.com/scientificreports/ Biomass sampling protocol. For a sampled Z. serrata tree, the tree was felled at the base and separated into four components in the field: main stem, large branches, twigs, and foliage. Fresh weight (kg) of each component was measured. Stem analysis was carried out for the main stem. The main stem was separated into 1-m sections. A stem disc was collected at the top of each section, and its fresh weight (kg) was measured. Large branches, twigs, and foliage were subsampled, and the samples of the three components were measured for their fresh weights (kg). All stem discs and samples were brought back to laboratory and dried at 65 °C until constant weight. The oven-dried stem discs and samples of large branches, twigs, and foliage were measured for their dried weights (kg). A ratio of dried weight to fresh weight for each component (i.e., stem, large branch, twig, and foliage) was calculated from the samples. For each component, the ratio was applied to convert the fresh weight of the component recorded in the field to its dried weight biomass (kg).
Statistical analysis. The above ground biomass (AGB, kg) of a sampled Z. serrata tree was defined as the sum of its four component dried weight biomasses: stem biomass (B stem , kg), branch biomass (B branch , kg), twig biomass (B twig , kg), and foliage biomass (B foliage , kg). For each sampled tree, stem biomass ratio (R stem = B stem / AGB), branch biomass ratio (R branch = B branch /AGB), twig biomass ratio (R twig = B twig /AGB), and foliage biomass ratio (R foliage = B foliage /AGB) were calculated with the four ratios summed to one. The four ratios were used to build the CRM for each biomass component. Two biomass ratio models were applied based on a nonlinear model 10 (Eq. 1) and a linear model (Eq. 2) relating a biomass component ratio to a tree attribute, where, c denoted a biomass component, and X was a predictor. Three predictors were considered: dbh, dbh 2 , and dbh 2 ·ht. The first predictor was the tree diameter, the second predictor represented tree basal area, and the third predictor represented tree volume. As a result, there were a total of six combinations of model and predictor for developing the Z. serrata CRM. For R foliage , preliminary data analysis and model fitting showed that the parameter β 1 was not statistically significantly different from zero for the six combinations of model and predictor. This suggested that foliage biomass was not significantly associated with dbh, tree basal area, and tree volume. Hence, following the suggestions by Jenkins et al. and Radtke et al., the six combinations of model and predictor were only fitted to R stem , R branch , and R twig 10,13 . As the four ratios should sum to one, R foliage was calculated by subtracting the sum of the other three component ratios from one, i.e., R foliage = 1 -(R stem + R branch + R twig ). As a result, the three component ratio equations (i.e., R stem , R branch , and R twig ) was integrated as a system for each combination of model and predictor. To properly develop such a system, one should consider that the component ratios were dependent and the residuals were correlated because the same tree gave the values to the three component ratios 22 . To account for potentially correlated residuals, Seemingly Unrelated Regression (SUR 12,23 ) was used to fit a system of the three component ratio equations for each of the six combinations of model and predictor. In particular, Nonlinear Seemingly Unrelated Regression (NSUR) was applied to Eq. (1), and Linear Seemingly Unrelated Regression (LSUR) was applied to Eq. (2). Comparison between the combinations was made by examining residual standard error (RSE) and residual plots. The best system for the three component ratio equations was chosen. All analyses were carried out in R using systemfit package 24,25 . Ethics declarations/protocol compliance. The experimental and field protocols of collecting plant materials in this study were performed in accordance with relevant institutional and national guidelines and regulations.

Model selection.
There was a total of six combinations of two SUR models and three predictors fitted to the Z. serrata dataset to build the CRM for the three biomass components. Selecting the final Z. serrata CRM model was based on: (1) residual plots, (2) comparisons of RSE, and (3) levels of significance of the estimated parameter β 1 in the fitted CRMs. Residual plots of the three NSUR models for the three predictors (i.e., dbh, dbh 2 , and dbh 2 ·ht) depicted clustering of residuals over a small range of predicted values (Fig. 1). The clustering of residuals was especially prominent for the NSUR models with dbh 2 and dbh 2 ·ht (Fig. 1d-i). For example, for the NSUR model with dbh 2 ·ht, the residuals ranged from − 0.15 to 0.15% for a predicted value of about 0.43 for R stem (Fig. 1g). This implied that predicted values of a component ratio for a tree attribute were very similar even though the actual observed values were different. This could be an issue when predicting a component ratio for a new tree.
On the contrary, residuals of three LSUR models for the three predictors were more dispersed over the range of predicted values (Fig. 2). Moreover, the range in the residuals of the three LSUR models was smaller than their NSUR counterparts. The LSUR model with dbh generally produced more homogeneously dispersed residuals without an obvious trend across the range of predicted values consistently for the three biomass component ratios (Fig. 2a-c) compared to the residuals from the two LSUR models with dbh 2 and dbh 2 ·ht (Fig. 2d-i).
Agreeing with the residual plots, the RSEs of the three NSUR models were consistently larger than their LSUR counterparts, which could be 20-60% larger depending of the component ratio ( Table 2). The discrepancy was particularly large for R stem . Moreover, for the two NSUR models with dbh 2 and dbh 2 ·ht, the estimated β 1 for the R twig were not significantly different from zero ( Table 2). Among the three LSUR models, the LSUR model with www.nature.com/scientificreports/ dbh generally had the lowest or comparable RSE than the two LSUR models with dbh 2 and dbh 2 ·ht (Table 2). Furthermore, its estimated β 1 for the three biomass component ratios were more highly significant than the estimated β 1 of the two LSUR models with dbh 2 and dbh 2 ·ht, i.e., smaller p-values (Table 2). Considering the consistency across the three assessment criteria, the LSUR model with dbh performed the best and was chosen to build the Z. serrata CRM system.

Component ratio model (CRM). The largest biomass component ratio of the sampled Z. serrata trees was
R stem with an average and standard deviation of 0.463 ± 0.118 (range = 0.264 to 0.604). The second largest biomass component ratio was R branch with an average and standard deviation of 0.399 ± 0.136 (range = 0.237 to 0.624). R twig had an average and standard deviation of 0.092 ± 0.033 (range = 0.052 to 0.144). R foliage was the smallest with an average and standard deviation of 0.046 ± 0.023 (range = 0.018 to 0.096). Despite obvious difference in the average values, the range of the four tree component ratios was fairly wide. Especially, the range showed large overlapping between R stem and R branch , and between R twig and R foliage . The fitted LSUR model had an overall R 2 of 0.73 suggesting that the model overall goodness of fit was good with about 73% of the total variance in biomass component ratios linearly explained by dbh. Fitted LSUR suggested a very high negative correlation between R stem and R branch linear models (− 0.834; Table 3), but a more moderate negative correlation between R branch and R twig linear models (− 0.341; Table 3). Thus, the moderate to high correlation between two component ratio models highlighted the need to apply SUR in model fitting.
In general, dbh explained the variance of each biomass component ratio relatively well with multiple R 2 ranging from 0.52 to 0.75 (Table 3). However, the linear relationship between dbh and R stem , R branch , and R twig was different. For R stem , the relationship was negative with an increase of 1 cm in dbh correlated with a decrease of 0.019 in R stem (p-value = 0.0005; Fig. 3a, Table 3). On the other hand, for R branch , the relationship was positive with an increase of 1 cm in dbh correlated with an increase of 0.022 in R branch (p-value = 0.0003; Fig. 3b, Table 3). Lastly, for R twig , the relationship was negative with an increase of 1 cm in dbh correlated with a decrease of 0.0045 in R twig (p-value = 0.0082; Fig. 3c, Table 3).

Discussion
Many past studies focused on modeling the relationship between dry weight biomass of tree components and tree attributes 26,27 . However, very few studies have modeled the relationship with ratio of biomass in each tree component to AGB 10,13 . Woodall et al. applied the CRM models 10 to estimate biomass and carbon content of trees in USA using the national forest inventory data. Our study contributes to the continuing modeling efforts to understand relationship between biomass component ratios and tree attributes. Our study is also unique in that it is the first to model CRM under SUR framework anticipating that there would be correlation between models. This has been supported in the results with correlation as high as 80%. Carvalho and Parresol suggested that it would be more realistic to consider component biomasses being dependent and residuals being correlated 28 . SUR should lower estimated variances of regression parameters, which means higher efficiency in estimating parameters and producing reliable prediction intervals 11,12 . Thus, results would be more reliably interpreted when applying SUR, and would in turns lead to more confidence in decision making such as designing effective silvicultural treatments. For that matter, CRM should be analyzed under the SUR framework as would other studies on dried weight biomass 22 .
For studies on biomass component ratio 10,13 and on dry weight biomass component 6,29 , nonlinear relationship in the form of exponential distribution has been used to relate biomass ratio or biomass to dbh. In contrast, the nonlinear model (Eq. 1) in our study had poor predictability. The residual plots suggested that the fitted nonlinear Table 3. Estimated parameters and properties of the final fitted linear seemingly unrelated regression models (Eq. 2). The estimated parameters are reported for the main stem, large branch, and twig biomass component ratios. RSE is residual standard error.  www.nature.com/scientificreports/ models predicted similar values of biomass component ratios for trees of different dbhs. On the contrary, the linear model (Eq. 2) had better predictability. A possible explanation that the linear model performed better could be due to small sample size from the young Z. serrata stand. That being said, the sampled trees covered a wide range in dbh from 3.9 cm to 19.7 cm suggesting the scope of inference for the final LSUR CRM models should be adequate. Hence, multiple model forms should be compared when building a CRM system for a tree species.
Past studies have assessed carbon sequestration potential of Z. serrata tree species on a per area basis, e.g., with forest inventory data for Taiwan 30 , with remote sensing data in an urban forest in USA 31 , and for urban forests in South Korea 32,33 . None of these studies has assessed how Z. serrata distributing biomass among its components on tree-level. While the sampled trees in our study are of the same age, the Z. serrata stand exhibited strong horizontal and vertical stratification with a wide range in dbh and ht. The fitted CRM showed clear tendency of dominant Z. serrata trees to allocate biomass into developing larger branches at the expense of main stem biomass. One would expect that investing in larger branches is for crown development. However, the final fitted CRM suggested otherwise as there was no statistically significant increase in foliage biomass in dominant trees. Moreover, twig biomass in dominant trees were less than that in smaller diameter trees. On the contrary, for intermediate or suppressed Z. serrata trees, majority of sequestrated carbon is allocated to developing main stem according to the final fitted CRM system. While Z. serrata is highly adaptive to grow in a range of environments, thus preferred as a reforestation species, they are prone to branching and forking 34 . In a study 35 , 82% of Z. serrata trees in a five-year-old plantation developed forks with 44% of the trees forked at 1.3 m and below while 39% of the trees above 1.3 m. This corresponds with our study in that Z. serrata tends to branch when local growing conditions are favorable for it to become dominant. Thus, pruning was necessary to increase wood utilization of Z. serrata trees unless seedlings were planted in high density or with genetic selection 34 . This is supported by our study especially for dominant Z. serrata trees, which should be pruned early to avoid undesirable wounds.
From an economical perspective, early pruning of dominant trees and planting seedlings in high density both incur additional operational costs. However, economic gain from planting in high density could potentially offset the additional costs of planting materials and labor. It is generally observed that planting seedlings in high density tends to limit individual tree diameter growth due to increase competition 36 . From our modeling results, we speculate that high planting density of Z. serrata would lead to greater allocation of biomass to main stem instead of forming large branches, i.e., less forking and branching. This would reduce the cost of pruning at early stand development and increase extraction ratio when the stand is mature for thinning operation. Extraction ratio is defined as the ratio of harvested wood transported out of a forest to total wood harvested 37 . Higher extraction ratio implies greater economic returns from production of wood products, which could also serve as long-term carbon storage or fuel for BECCS through methods such as oxyfuel combustion 5 . However, an in-depth economic study such as net present value analysis is necessary to fully understand the implications on stand-and landscape-level 38,39 . Nevertheless, our study shows that modeling biomass allocation strategy of Z. serrata would have economical implication for Taiwan forestry as the tree species will continue to be important in reforestation effort.
It would be fairly easy to apply the developed CRM system in our study to assess biomass allocation in a Z. serrata plantation with a hybrid approach 6 . For a Z. serrata tree in a sample plot, its volume is first estimated and converted to total tree dry weight biomass with the biomass expansion factors from Lin et al 20 . The biomass component ratios of the four tree components (R stem , R branch , R twig , and R foliage ) are predicted from its dbh accordingly with our fitted LSUR models (Table 3). From which, one could then estimate dried weight biomass of the four tree components of the sample tree, which in turns could be expanded to per area basis with appropriate expansion factors associated with the sample plot 40 .
Despite that our study was carried out in a single even-aged stand, it is the first to suggest that there are significant differences in biomass allocation strategy for Z. serrata trees of different sizes at the early stage of stand development. Most of the Z. serrata plantations established during the Taiwan National Reforestation Program should be currently about 20 years old. We conjecture that dominant trees that are already in the main canopy in early stand development stage will likely continue the same growth trajectory and biomass allocation strategy, and so would the suppressed trees. However, future study should resample in the same study site, which is now 23 years old, to test our hypothesis.

Conclusion
This study is the first to model biomass allocation strategy of planted Z. serrata trees. It is one of the few studies to model biomass allocation with the CRM approach, and is also the first to model CRM under SUR framework to properly account for correlations between models. Our developed CRM could also be used to approximately predict tree component biomasses of a Z. serrata plantation when only carbon estimate per unit area is available. For example, a Z. serrata plantation stores about 452.2 t CO 2 ha −115 . Of this amount, our CRM models suggest that 209.4, 180.4, 41.6, and 20.8 t CO 2 ha −1 are stored in main stem, large branch, twig, and foliage, respectively. Contrary to other studies, our results supported a linear relationship between biomass component ratios and dbh instead of a nonlinear relationship. The fitted linear relationship suggests that Z. serrata trees in the main canopy have larger sized crowns because of the tendency in forking and branching, which could more effectively compete for resources and suppress development of the other trees. Pruning is necessary not only to improve wood utilization potential of dominant Z. serrata trees by reducing knots but also to allow other trees in a stand to develop their utilization potential. Future work could sample Z. serrata trees across stand development stages and elevation to examine whether there is any change in biomass allocation strategy of the tree species under different stand age and growing conditions. Moreover, studying potential effects of various silvicultural treatments on biomass allocation strategy of the tree species and economic tradeoff could lead to better planning of wood utilization in long-term carbon storage or bioenergy production.