Above-ground biomass equations for Pinus radiata D . Don in Asturias

Aim of the study: The aim of this study was to develop a model for above-ground biomass estimation for Pinus radiata D. Don in Asturias. Area of study: Asturias (NE of Spain). Material and methods: Different models were f itted for the different above-ground components and weighted regression was used to correct heteroscedasticity. Finally, all the models were refitted simultaneously by use of Nonlinear Seemingly Unrelated Regressions (NSUR) to ensure the additivity of biomass equations. Research highlights: A system of four biomass equations (wood, bark, crown and total biomass) was develop, such that the sum of the estimations of the three biomass components is equal to the estimate of total biomass. Total and stem biomass equations explained more than 92% of observed variability, while crown and bark biomass equations explained 77% and 89% respectively.


Introduction
Tree biomass estimation is needed for the sustainable planning of forest resources, so the development of biomass equations has been, and indeed remains, one of the main lines of work of many researchers (Cunia, 1986(Cunia, , 1988;;Waring and Running, 1998).In addition, the importance of forests as carbon sinks is closely linked to the amount of biomass available in them, which makes it essential to have equations to quantify the contribution of forest areas to the global carbon cycle (IPCC, 2003).
The first studies arose from the need to estimate the biomass productivity of different species, the studies on Larix decidua and Picea abies in Switzerland (Burger, 1945(Burger, , 1953) ) being the most notable.Subsequently research began to focus on the determination of the dry weight of tree components (e.g.wood, bark, branches), especially those components of greater importance for forest companies.The fields of ecology and physiology also showed an interest in this research line and con-tributed to the development of sampling methods, especially those to quantify foliage (Kittredge, 1944;Ovington, 1957).Many forest studies develop equations using regression techniques for specific geographic areas and tree species.Biomass equations relate tree biomass (kg) or stand biomass (t/ha), as well as their components, with easily measurable variables.Tree variables commonly used are diameter at breast height (d) and total height (h), while in stand equations, stand variables like stocking (trees/ha), basal area (G) and dominant height (H 0 ) are used.
Studies of the estimation of forest biomass have increased in recent decades, taking on importance and incorporating a large number of species and stand structures (Pardé, 1980;Zeide, 1987;Waring and Running, 1998).Zianis et al. (2005) made a review of biomass and volume equations for tree species in Europe and found 607 biomass equations, although only a minority relevant to Southern Europe.In many cases, these studies go beyond quantification or characterization, and cover more global aspects, such as those concerning the harvesting and further use of the wood (for pulp, firewood, biomass, etc.) and their application in research (e.g.study of carbon cycle, nutrient balance of forest systems).In recent years, remote sensing has proved to be an effective tool in quantifying forest biomass (Montes et al., 2000;Drake et al., 2002), although this method needs real data relating to tree size, so it is always necessary to quantify individual trees to later determine forest biomass and express it per unit area (Diéguez-Aranda et al., 2003).
Several studies on radiata pine (Pinus radiata D. Don) biomass have been developed, most of them conducted in New Zealand, Chile and South Africa (e.g.Madgwick, 1983, Guerra et al., 2005;Moore, 2010).In Spain, tree biomass equations for this species have been fitted in Galicia (Balboa-Murias et al., 2006), along with general models for the country as a whole (Montero et al., 2005).In addition, other studies have been conducted that relate biomass production of radiata pine with soil properties or nutrient accumulation (Rey et al., 2001;Merino et al., 2003).Recently, Balboa-Murias et al. (2006) carried out a study of the temporal variations and distribution of carbon storage in radiata pine and maritime pine pure stands under different silvicultural management.Given these findings highlighting the variability of growth depending on a number of different factors, it is evident that it is of interest to develop models which are specific to different geographical areas, with their own specific soil, meteorological and geographical conditions, in this case Asturias.
This work is part of a total quantif ication of the potential of radiata pine for carbon storage in Asturias and for the use of residues (thinnings and final fellings) for energy.The aim of this study was to develop equations to predict the biomass of above-ground components (wood, bark and crown) for radiata pine in Asturias.The analysis required the use of appropiate statistical techniques to correct the heterogeneity of variance of the residuals in the adjustment of different biomass components equations, as well as a simultaneous adjustment to ensure additivity between the estimates of biomass equations of the different components and the total tree biomass equation.

Data
The data used to develop above-ground biomass equations was based on a sample of 27 radiata pine trees in five stands managed by Forest Services in Asturias (sample size was dictated by the availability of resources).Trees from all diameter classes were selected (minimum diameter equal to 5 cm), avoiding trees with severe defects or imbalances due to their proximity to stand edges, forest tracks, etc.
Before felling, diameter at breast height (d) and total height (h) were measured.After felling, total height measurements were verified and trees were destructively sampled, separating above-ground biomass into the following components: stem (logs with bark with a thin-end diameter of 7 cm), branches of a diameter larger than 7 cm, thick branches (2-7 cm) and branches of a diameter less than 2 cm.
The fresh weight of each component was measured in the field with a portable balance (accurate to 0.1 kg).Disks, with bark, were cut from each stem at regular intervals.These disks together with representative samples of all tree components were then taken to the laboratory where disks were separated into wood and bark, and branches with a diameter of less than 2 cm were further separated into needles, twigs (diameter less than 0.5 cm at the insertion) and thin branches (0.5 to 2 cm diameter).
Subsamples of the different biomass components were oven-dried at 105°C to determine field moisture content, which was then used to convert fresh weight to dry weight.
Table 1 gives the descriptive statistics of diameter at breast height, total height, total dry weight and dry weight of the different components (incorporating the crown as the sum of branches and needles).

Adjustment of biomass equations
Once the dry weight of the different components had been determined, different models were developed to relate dry weight with measurable tree variables.Diameter is the most common variable used in aboveground biomass equations [w = f(d)], although the inclusion of tree height as the second predictor variable [w = f(d, h)] significantly improves the adjustments (Wang, 2006), hence both were used for this study.
A variety of models have been used in the estimation of tree biomass (Cunia and Briggs, 1984;Reed and Green, 1985;Reed et al., 1996), although all are derived from three mathematical models: linear, non-linear with additive error and non-linear with multiplicative error (Pardé, 1980;Snodow, 1985;Parresol, 1999): In this work, the first step was to fit an individual model for every component, using as independent varia-Tree biomass of Pinus radiata in Asturias bles the diameter and total height of individual trees.In the second step, in order to ensure additivity between the sum of the estimates of the biomass equations of the different components and the total biomass equation, and considering it unrealistic to assume that the errors of each equation are not correlated (Borders, 1989;Parresol, 1999Parresol, , 2001)), the equations of all the components of above-ground biomass were fitted simultaneously.
Initially, models for every sampled component were fitted, producing unsatisfactory results for branches, twigs and needles.In addition there were problems of convergence and high multicollinearity in the simultaneous adjustment, probably due to the sample size.Branches, twigs and needles were therefore combined in a single component named "crown", which resulted in increased accuracy of the estimations.This corresponds with the three components (wood, bark and crown) into which above-ground biomass is usually divided (Evert, 1985;Parresol, 1999).

Individual adjustment
Linear and non-linear models were fitted for each component considered in this study (wood, bark and crown).Linear models were fitted using the method of least squares with the REG procedure of SAS/STAT ® software (SAS Institute Inc., 2004a).Non-linear models were fitted by least squares with the SAS/STAT ® NLIN procedure, using the iterative method of Gauss-Newton (Hartley, 1961).
Normally, biomass data exhibit heteroscedasticity, that is, the error variance is not constant (Parresol, 2001).Although it is true that in the presence of heteroscedasticity the least squares estimates of the parameters remain unbiased and consistent, they are inefficient, and therefore their standard errors are not correct and the usual signif icance test cannot be applied.To confirm the presence of heteroscedasticity, the graphics of studentized residuals vs. estimated values were analyzed and White (1980) and Breusch and Pagan (1979) tests were applied.Heteroscedasticity was corrected by weighted regression (Schlaegel, 1982;Cunia, 1986;Parresol, 1999), applying to every observation a weighting equal to the inverse of the variance of the residuals, using the potential function (Neter et al., 1996).The value of the exponent k was calculated using the optimization methodology proposed by Harvey (1976).Different weighted factors were tested (d -k , d 2 h) -k depending on the independent variables of each model.
The other usual problem in the adjustment of biomass equations is multicollinearity, which arises when the independent variables in the model are themselves correlated.The presence of severe multicollinearity makes parameter estimates unstable, i.e. they depend too much on the particular dataset used in the adjustment (Myers, 1990).To analyze the presence of multicollinearity in each evaluated model, the condition number (CN) was used.According to Belsey (1991, pp: 139-141), a CN < 10 indicates that collinearity is not a major problem, if 30 ≤ CN ≤ 100 there are problems of multicollinearity and if 1,000 < CN < 3,000 the problems are severe.

Simultaneous adjustment to ensure additivity
One of the most important properties of biomass equations for the different components is that the sum of their estimations should be equal to the estimation provided by the total biomass equation; a property known as additivity.Studies by Kozak (1970), Chiyenda and Kozak (1984) and Cunia andBriggs (1984, 1985) resulted in three different methods for forcing additivity in a lineal system of biomass equations.Parresol (1999) concluded that the most flexible and general is the method based on simultaneous adjustment using Zellner's (1962) seemingly unrelated regression (SUR).
It allows the use of different mathematical models, with different independent variables and weighting factors, for each fraction.Furthermore, this method corrects the problem of inherent errors between the estimations of each model.In this case, the independent variables of the total biomass equation are all the independent variables of the rest of the models and additivity is ensured by setting restrictions on the regression coefficients.SUR takes into account the existing correlations between biomass components and results in lower variance.The greatest benefit in using SUR methodologies is that confidence and prediction intervals for biomass estimates are reduced (Parresol, 2001).In the case of systems of non linear equations, the simultaneous adjustment requires the use of NSUR (Nonlinear Seemingly Unrelated Regressions) (Parresol, 2001).In this study, once individual models were selected for each component the parameters were estimated again by NSUR methodology, using the MODEL procedure of SAS/ETS ® (SAS Institute Inc., 2004b).In the simultaneous adjustment, the same weighting factors for every component were used, ensuring the additivity of biomass estimations, thus allowing the calculation of the biomass of two or more components, or total biomass, simply by adding together the dry weight of each component.

Model evaluation
Comparison of the performance of the models for each component was based on numerical and graphical analyses of the residuals.For this purpose, the following statistics were calculated: coefficient of determination (R 2 ), root mean square error (RMSE) and Akaike's information criterion differences (AICd).
Besides the statistics above, one of the most efficient ways of evaluating model performance is by visual inspection.Therefore, scatter plots of studentized residuals against the estimated values were analyzed to check for the presence of heteroscedasticity.

Individual adjustment
Once the different models for each component were adjusted, the best model was selected in each case based on goodness-of-fit statistics and graphical analyses.After the graphical analysis and the White (1980) and Breusch and Pagan (1979) tests, heteroscedasticity was detected in all the selected models for every component, so these were f itted again using weighted regression, selecting the weighting factor according to the optimization methodology proposed by Harvey (1976).Results of the heteroscedasticity tests and weighting factors used in the individual adjustments are shown in Table 2 and the correction can be observed (p-value greater than 0.05).Scatter plots of studentized residual against predicted values using the weighting factors also demonstrated heteroscedasticity to be corrected through weighted regression.
Table 3 shows the parameter estimates and their corresponding approximated standard errors for the selected models for each component using weighted regression.All the parameters were found to be significant at the 5% level.Additionally, this table shows the condition numbers and goodness-of-fit statistics calculated from the residuals obtained with the weighted w i : dry weight of component i (kg).β j : parameter of the model.d: diameter at breast height (cm).h: total height (m).
regression.The condition numbers calculated allow the conclusion to be drawn that the models have no problems of multicollinearity.

Simultaneous adjustment
Once each biomass component was individually adjusted, the best models selected were simultaneously fitted using NSUR methodology in order to ensure the additivity of the biomass components.
All the independent variables that appear in the models of wood, bark and crown biomass are included in the model of total biomass; the additivity is guaranteed by imposing restrictions on the parameters.The models of different components were adjusted using the same weighting factors as in the individual adjust-ments (Table 2).Table 4 shows the parameter estimates, significance tests and goodness-of-fit statistics for the models obtained by simultaneous adjustment of the system of the four biomass equations.All parameters were significant at the 5% level.The homogeneity of error variance was confirmed by graphical analyses of the studentized residuals and White and Breusch-Pagan tests.
The equations finally proposed for the estimation of above-ground biomass of radiata pine trees in Asturias are then:  [4] where w i is the dry weight of each biomass component i (kg), d the diameter at breast height (cm) and h the total height (m).Fig. 1 shows the graphics of the observed values of the different biomass components and of the total biomass compared to predicted values using the corresponding equation.

Discussion
The simultaneous adjustment provided a sensitive improvement in most of the goodness-of-fit statistics.Other authors have also verified such an increase in some of the statistics, as well as a slight decrease in others (Paulo et al., 2002;Carvalho and Parresol, 2003).
Although the changes in the goodness-of-fit statistics were small and the parameters estimated through simultaneous adjustment were similar to those obtained by individual adjustment, the parameter standard errors decreased noticeably (a direct result of taking into account the contemporaneous correlation between the components), and hence the parameter estimation was more efficient.Similar results have been obtained in other studies which verif ied that NSUR methodology improves the accuracy of biomass predictions (Parresol, 1999(Parresol, , 2001;;Carvalho and Parresol, 2003).
The best adjustments were obtained in the equations of wood and total biomass.For crown biomass, variance in its estimation was greater, in relative terms, than that obtained in the estimations of wood biomass, which is due to the variability of the crown structure, the number of branches and variation in wood density along branches (Pardé, 1980).Maybe the inclusion of other independent variables (e.g.live crown length) could improve the quality of the adjustment, although it involves a great sampling effort to be able to apply the model.The condition number in the simultaneous adjustment was 51.7, indicating that the model does not have serious problems of multicollinearity.

Conclusions
In this study a system of equations of above ground biomass for radiata pine in Asturias was developed through simultaneous adjustment with the generalized least squares method known as NSUR (Nonlinear Seemingly Unrelated Regressions), which ensured the additivity of the equations.A system of four biomass equations (wood, bark and crown biomass and total biomass) was developed, such that the sum of the estimations of the three biomass components is equal to the estimate of total biomass.The wood and total biomass equations explained more than 92% of the variability of the observed data, while the bark and crown biomass equations accounted for 77% and 89% respectively.In further studies it would be advisable to increase the sample size in order to obtain a system of equations which consider all the components initially sampled individually (stem, bark, branches, twigs and needles) and thus improve the quality of the fit.However, the results obtained in the current work allow the estimation of the dry weight of different components and total dry weight with sufficient accuracy for the aim of this study, and, furthermore, these equations have better behavior than the models that have existed until now for this species in Spain.

Figure 1 .
Figure 1.Scatter plots of observed values against predicted values of the different biomass components and total biomass.

Table 1 .
Descriptive statistics of sample of 27 trees used in the adjustment

Table 2 .
Tree biomass of Pinus radiata in Asturias 411 Weighting factors used in the adjustment by weighted regression of different component biomass equations and results of heteroscedasticity tests

Table 3 .
Selected equations for each component, parameter estimates, associated approximate standard errors, condition number and goodness-of-fit statistics of the individual adjustments using weighted regression
wi: dry weight of component i (kg).βj:parameter j of component i (wood, bark or crown).d:diameter at breast height (cm).h:total height (m).R 2 : coefficient of determination.RMSE: root mean square error.wi