1. Introduction

Radiata pine (Pinus radiata D. Don) is a common plantation species in a number of Southern Hemisphere countries, including New Zealand, Chile, Australia, and South Africa, with more than 4 million hectares of plantations having been established (Lavery and Mead, 1998). In New Zealand, the total area of radiata pine plantations is approximately 1.6 million ha, which is 89 per cent of the total area of exotic plantations (Ministry of Agriculture and Forestry, 2007). Since the 1960s a considerable number of biomass studies have been conducted in this species in New Zealand, Australia and South Africa. These were generally undertaken to estimate forest productivity (as measured by dry matter production), nutrient budgets (e.g. Baker and Attiwill, 1985; van Laar and van Lill, 1978; Will, 1964) and how these are affected by a range of site and management factors (e.g. Barton, 1984; Cromer et al., 1985; Mead et al., 1985), as well as the partitioning of biomass among foliage, branches, stems and roots (e.g. Beets and Pollock, 1987). Much of this earlier research into radiata pine biomass has been summarised by Madgwick (1994).

Recently, there has been a renewed interest globally in biomass research in most tree species, including radiata pine, due to the need to predict forest carbon stocks and the potential amount of biomass available as a source of energy. One of the main drivers for this has been the Kyoto Protocol of the United Nations Framework Convention on Climate Change (UNFCCC), which stipulates mechanisms whereby post-1990 carbon storage in forests can be used as an allowable carbon sink to offset some greenhouse gas emissions (IPCC, 2003). Under Article 3.3 of the Kyoto Protocol, net carbon stock changes on land subject to afforestation, reforestation, and deforestation must be estimated each year over the defined commitment period from information on land area and its corresponding carbon density. To help meet its obligations under the Kyoto Protocol, the New Zealand Government established the Emissions Trading Scheme (ETS), which introduced a price for greenhouse gases in order to provide an incentive to reduce emissions and increase forest sinks (New Zealand Government, 2008). Forestry is the first sector to enter the ETS and therefore a robust method is needed for determining the carbon stocks in forests registered under this scheme.

A number of approaches exist for estimating tree biomass. These include physiologically-based models, such as 3PG (Landsberg and Waring, 1997) that use information on temperature, solar radiation and rainfall to predict total biomass; biomass expansion factors (e.g. Lehtonen et al., 2004) that convert stem volume into whole-tree biomass; and allometric equations that predict total biomass from easily measured variables, such as diameter at breast height and height (e.g. Parresol, 1999; Ter-Mikaelian and Korzukhin, 1997; Zianus et al., 2005). Allometric equations have been developed for predicting the biomass of individual radiata pine trees and stands within a number of previous biomass studies (see Keith et al. 2000 for a catalogue of equations). However, many of these equations are stand- or site-specific (e.g. Baker et al., 1984; Dargavel, 1970) and so potentially introduce high uncertainty when applied at a national level where a wide range of site types, stand structures and age classes exist. A number of individual-tree equations were developed by Snowdon et al. (2000) using a dataset consisting of 351 trees compiled from various Australian and New Zealand biomass studies. However, additional data are known to exist from studies performed in New Zealand, which could be used to develop nationally-applicable functions for predicting biomass. Many of these data have already been used to develop functions for different biomass components (e.g. Madgwick, 1983; 1985), but not total above-ground biomass.

Therefore, the objective of this paper was to develop functions that can be used to predict the total above-ground biomass of individual radiata pine trees growing in New Zealand. Mixed effects models are used to account for any peculiarities of the stands investigated or different measurement methodologies used in the different studies, following a similar approach taken by Wirth et al. (2004).

2. Materials and methods

2.1. Radiata pine biomass data

A database was compiled consisting of biomass measurements for 637 individual radiata pine trees from 13 different studies (Tab. I). Only those trees that had both diameter at breast height (D) and height (H) measurements were considered for subsequent analysis. These trees ranged in age from 2 to 42 y, in D from 7 to 806 mm and in height from 1.47 to 42.1 m. However, there was a relatively large number of smaller trees in the dataset (Fig. 1). Of the 637 trees in the dataset, 620 had information available on the stem, bark, branch and foliage components in addition to total above-ground biomass. The methods used to measure biomass differed between studies. In general, the biomass of smaller trees was determined by sub-dividing the whole tree into its components (i.e., stem wood, bark, branches and foliage) and weighing each of these components. In larger trees, the basal area of every branch in the tree was typically measured and a sub-sample of these branches were weighed. The relationship between branch component mass and branch basal area was determined and used to estimate these biomass components for the whole tree.

Figure 1 Distribution of (a) diameter at breast height and (b) total height of the 637 radiata pine trees used to develop the biomass equations.

The majority of records in the database come from studies that have been carried out in the North Island of New Zealand where rainfall is in excess of 1000 mm per annum. Approximately half of the observations are from the studies conducted by Beets and Pollock (1987) and Smith et al. (1994). Most of the studies in eastern areas of the South Island, where rainfall is often less than 700 mm, have measured the biomass of trees less than 10 years old (e.g. Bandara, 1997; Clinton, 1990; Mason, unpublished data). In addition, a number of other studies have reported stand level estimates of biomass (see Madgwick, 1994 for details) or refer to individual tree measurements (e.g. Will, 1965); however, it was not possible to obtain these individual tree measurements in the course of assembling the database.

2.2. Development of biomass equations

A preliminary visual inspection of the data was undertaken to check for outliers and to investigate potential predictive relationships. Regression equations were developed to relate total above-ground tree mass (AGB) to different measures of tree size, including diameter at breast height (D [cm]), height (H [m]) and D²H. The first of these equations was based on the traditional non-linear allometric relationship which has the following form: where: AGB = total biomass or one of its components, X_l= tree dimension variable (l = 1,..., p), and β_l are the model parameters (l = 0,1,..., p). In equation (1) either an additive error term or a multiplicative error term can be assumed (Parresol, 1999). Because biomass data generally exhibit non-constant variance, a common approach taken when developing predictive equations is to convert equation (1) into a linear equation with an additive error term by taking the natural logarithm of both sides, i.e., where ε is the residual error term. This transformation linearises the allometric relationships and stabilises the error variance (Baskerville, 1972). In this paper, the first set of equations fitted to the radiata pine data followed this approach, which also enabled parameter estimates and goodness of fit statistics to be compared with those from previous biomass equations developed for radiata pine (e.g. Snowdon et al., 2000) and other species (e.g. Ter-Mikaelian and Korzukhin, 1997; Parresol, 1999; Zianus et al., 2005). The following equations with this form were fitted to the data: To account for methodological and site differences between studies, random fluctuations in the model parameters were included, such that or in vector notation: where AGB_ij is the above ground biomass at the ith measurement in the jth study, b_lj is the random deviation from the mean regression coefficient β_l due to study j, and ε_ij is the residual error ( ∼ N(0, σ²). Models were fitted using the linear mixed effects model function within the R software program (R Project Core Team, 2009). Following the approach taken by Wirth et al. (2004), random effects representing differences between studies were only included for the intercept and coefficients of the first order terms in Equations (3)–(9). Random effects were assumed to be independent of the residual error and normally distributed with covariance matrix Ψ ( b _j ∼ N( 0 , Ψ )). While random effects for the different studies were assumed to be independent, random effects for an individual study were not assumed to be independent, and therefore the variance-covariance matrix for the random effects was positive definite symmetric (Pinheiro and Bates, 2000). Equations (3), (5) and (8) were also re-fitted to data from trees older than five years to allow direct comparison with the equations developed by Snowdon et al. (2000). Predictions from the two sets of equations were compared using Lin’s concordance correlation coefficient (Lin, 1989; 2000), which considers both the correlation and agreement between the sets of predictions.

To correct for the bias that occurs when estimates from logarithmic equations are back transformed to arithmetic units, the modified smearing estimate developed for linear mixed models by Wirth et al. (2004) was used. In matrix notation, the modified smearing estimate (i.e. the estimated expected value) of the total above ground biomass for a new tree (AGB_new) is given by: where X _new  = (1,X_1new,X_2new,.....,X_pnew), β̂ is the vector of estimated fixed effects from the model and cf₁ and cf₂ are correction factors that are given by: where b̂ _j is the vector of estimated random effects associated with study j, Z _new is the vector of predictor variables for which random fluctuations were calculated, and ε̂_ij is the vector of estimated residuals from the model. From equation (12b) it can be seen that the value of cf₂ depends on the vector Z _new . For this study the mean value of cf₂ was calculated for the underlying dataset used to develop the models.

As an alternative to correcting for the bias that occurs when converting estimates from logarithmic models back to arithmetic units, a series of linear and non-linear models were fitted to data on the original arithmetic scale. Specifically, the following models were fitted to the data: Weighted analyses were used when fitting the models given by equations (13)–(20) due to the non-constant error variance. Weightings for each observation were determined by assuming that in each model the variance (Var(ε_ij)) was a power of the absolute value of the variance covariate (v_ij), i.e. . When fitting equations (13) and (18), D was chosen as the variance covariate, while in the remaining equations it was D²H. The value of the variance function coefficient δ was determined using varPower constructor within the R software program (R Development Core Team, 2009) and a common value was assumed across all studies. The linear models (Eqs. (18)–(20)) were fitted following the same approach taken for the models fitted to the transformed data, with random effects included for both parameters in each model. For the non-linear models (Eqs. (13)–(17)), random effects were included for each parameter following the same approach taken for the linear models; however, in order for the algorithm to converge, the variance-covariance matrices for the random effects had to be simplified by assuming a block diagonal structure, rather than a symmetric structure (Pinheiro and Bates, 2000).

Figure 2 Relationship between total above ground biomass (AGB) for an individual tree and tree height (H), diameter at breast height squared (D2) and D2H.

2.3. Model evaluation

Comparisons between models were made using several goodness-of-fit statistics. Models with the same response variable (i.e., those based on the transformed allometric equation or those fitted on the arithmetic scale using weighted analysis) were compared using the Akaike Information Criterion (AIC). The AIC was calculated from the results of model fitting as  − 2 log-likelihood + 2 P, where P is the number of parameters in the model (including elements of the variance-covariance matrix). In order to calculate AIC, models were re-fitted and the parameter estimates obtained using the maximum likelihood method (see Pinheiro and Bates, 2000). Using the AIC statistic, the ’best’ model was taken to be the one with the smallest value of AIC. Because AIC depends on the scale of the response variable, it cannot be used to compare models with different response variables. Therefore, Furnival’s index (FI) (Furnival, 1961) was used to compare transformed models with weighted linear and non-linear models. This index was calculated as: where f^′(Y) is the derivative of the dependent variable with respect to total above-ground biomass (the square brackets signify the geometric mean) and RMSE is the root mean square error of the fitted equation. For untransformed models, FI is equivalent to the RMSE and a rescaled RMSE for transformed or weighted models. A fit index (R²) was also calculated for each model using the following equation (Parresol, 1999): where TSS is the total sum of squares and RSS is the residual sum of squares. These were calculated as: where is the mean value of above ground biomass across all trees in the dataset and AĜB_i is the predicted value for the ith tree. In calculating these fit indices, model predictions were based on the fixed effects and predictions from the logarithmic models were back transformed and corrected for bias.

Models were also evaluated by examination of residuals, with an emphasis on investigating the existence of any possible trends with age, stand density and wood density as these factors could contribute to site level differences in total above ground carbon for a given tree diameter and height. Mean wood density for the whole stem was calculated by dividing total oven-dry stem mass by the total stem volume under bark as determined from sectional measurements.

3. Results

Visual examination of the data (Figs. 2 and 3) revealed the existence of strong relationships between AGB and both D² and D²H on the original scale of measurement. The relationship between AGB and tree height was clearly non-linear and there was also evidence of increasing variance in AGB with increasing values of all three tree size variables. On the logarithmic scale strong relationships were observed between log(AGB) and log(H), log(D) and log(D²H) with a reduced heterogeneity in variance compared with the arithmetic scale (Fig. 3).

Figure 3 Relationship between total above ground biomass (AGB) for an individual tree and tree height (H), diameter at breast height squared (D2) and D2H on a logarithmic scale for x and y axes.

3.1. Models fitted to transformed data

On the transformed scale, the models predicting AGB as a linear function of D (Eq. (3)) or the combined variable D²H (Eq. (8)) were the poorest performed as measured by AIC (Tab. II). Comparison of the residuals from Equation (3) against predicted values (Fig. 4) (and similarly for Eq. (8)–not shown) revealed evidence of a systematic trend which indicated that the relationships were possibly not linear and that an additional higher order term could be required. The corresponding models containing additional quadratic terms (i.e., Eqs. (4) and (9)) had lower values of AIC than their linear counterparts. Further improvement to the model containing linear and quadratic diameter terms was obtained by the inclusion of tree height terms (Eqs. (6) and (7)). Overall, the ‘best’ model on the logarithmic scale in terms of AIC (i.e., the one with the lowest value) was that given by Equation (7) (Tab. II). Values for the first correction factor for the smearing estimate (cf₁) ranged from 1.0254 to 1.0474, while the mean value of cf₂ for the underlying dataset ranged from 1.0527 to 1.0970. Values of the fit indices (R²) based on the corrected back transformed predictions from the models ranged from 0.92 for the model based on a linear function of D (Eq. (3)) up to 0.95 for the models given by Equations (6) and (7), indicating that the latter are able to explain more of the variation in total above ground biomass.

Figure 4 Comparison of residuals against predicted values of above ground biomass (on the logarithmic scale) for Equation (3) which contains diameter at breast height as the sole predictor. The solid line represents a locally-weighted smoothing function (lowess) fitted to the data.

3.2. Weighted arithmetic models

On the arithmetic scale, the linear model (Eq. (20)) containing the combined variable (D²H^0.5) had the lowest value of AIC (5237). The linear models containing D² and D²H (Eqs. (18) and (19)) had the highest values of AIC (5721 and 5571, respectively). However, the higher weighting (δ = 1.00) given to results from smaller trees in the model containing D² (Eq. (18)), coupled with some apparent non-linearity in the relationship with AGB, meant that this model considerably under-predicted AGB in larger trees (Fig. 5a). It was possible to rectify this problem by adding a higher order term to the model (either D³ or D⁴). The resulting models exhibited more satisfactory behaviour for larger trees (Fig. 5b). The three and four parameter nonlinear models containing D and H (Eqs. (16) and (17)) also had relatively low values of AIC (5347 and 5349, respectively). For these models the value of the variance function coefficient δ was approximately 0.5 indicating the that variance was linearly proportional to D²H. However, with this weighting there was some indication that Equations (16) and (17) under-predicted when AGB was greater than 1500 kg (data not shown), although there were fewer than 30 observations for trees with a total biomass greater than this value.

Figure 5 Comparison of measured above-ground biomass with that predicted by (a) Equation (18) and (b) the same Equation but with a higher order diameter term (D3) added.

3.3. Model evaluation

Comparison of the values of Furnival’s index showed that overall the transformed models performed better than the weighted arithmetic models (Tab. II), with Equations (6) and (7) having the lowest value of FI (10.67 and 10.66, respectively) of all the models developed. These models were able to explain approximately 95% of the variation in AGB, while the best weighted model (Eq. (20)) was able to explain approximately 97% of the variation. All three of these models had both tree height and diameter as inputs. Of the three models (Eqs. (4), (13) and (18)) that only had diameter as an input, Equation (4) was the best performing and was able to account for 95% of the variation in AGB.

Figure 6 Comparison of residuals from the model given by Equation (7) against tree age, stand density and wood density.

For the best overall model in terms of Furnival’s index (Eq. (7)), there was some evidence that it over-predicted biomass in stands younger than five years old (Fig. 6). Many of these trees were sampled by Bandara (1997) and this study had the largest random deviations from the model parameters of all the studies considered. An attempt was made to include a stand age term in this model, but this term was not significant (p > 0.05). This was possibly due to the high correlation between stand age and tree height (r = 0.92) which meant that stand age was already implicitly included in the model. There was no clear evidence of a trend in the residuals with stand density (Fig. 6). Mean wood density was only able to be calculated for 258 of the 637 trees in the biomass dataset (Tab. IV), with most of these trees from the studies conducted by Beets and Pollock (1987) and Smith et al. (1984). For these 258 trees, mean stem density ranged from 297 kg m^-3 up to 534 kg m^-3 with an average value of 344 kg m^-3. There was some indication that the model under-predicted AGB for the nine trees with the highest values of wood density, but there were few observations at higher values of wood density upon which to assess the presence of any systematic trend.

3.4. Comparison with previous models

Re-fitting Equations (3), (5) and (8) to data from trees greater than five years of age resulted in the following parameter estimates: These were very similar to the parameter estimates for the equations developed by Snowdon et al. (2000). Comparison of the predictions from the two sets of equations for the range of tree sizes in the database assembled for this study showed that there was almost perfect agreement between them, as values of the concordance correlation coefficient all exceeded 0.99.

4. Discussion and conclusions

This meta-analysis built on the work of Magwick (1983) and Snowdon et al. (2000) and has produced a series of models that can be readily applied to predict the total above-ground biomass of individual radiata pine trees. These equations are intended to be applied at a national level, primarily for carbon accounting purposes. At a national level, unbiased estimates of carbon stocks are required, but biomass estimates obtained by applying these national level equations to individual stands may be biased as the relationship between tree dimensions and biomass may differ between stands. In developing these equations a mixed modelling approach was used to account for random perturbations caused by methodological and site differences, which had not been done in the studies by Magwick (1983) and Snowdon et al. (2000) that used ordinary least squares regression. However, comparisons showed that the functions developed in this study were very similar to those developed by Snowdon et al. (2000) using a dataset containing measurements from 351 radiata pine trees greater than five years in age.

The mixed modelling approach has been used by Wirth et al. (2004) to develop functions for Norway spruce (Picea abies (L.) Karst.) using data obtained from studies conducted by different authors. Wirth et al. (2004) assumed that differences between studies reflected methodological and site differences. In the radiata pine dataset, methodological differences were often related to tree size, with the biomass of small trees determined from weighing the entire tree, whereas sub-sampling was normally undertaken to determine the biomass of older larger trees. However, such methodological differences may not always be captured by the random effects structure assumed in the models as several of the studies sampled trees from across an age series (e.g. Beets and Pollock, 1987; Madgwick et al., 1977). Site differences may reflect different silvicultural treatments such as initial spacing and thinning intensity as well as fertiliser treatments; many of the data available were from studies investigating the effects of such treatments on biomass accumulation and nutrient uptake (e.g. Barton, 1984; Beets and Pollock, 1987). In some instances different treatments were applied within a study, but this level of intra-study variation was not included as a random effect as records linking individual trees to a particular experimental treatment were incomplete.

Site differences could also be due to broad scale differences in stem form and wood density. In New Zealand, wood density of radiata pine can vary by up to 25–30% across the latitudinal range over which it is grown (Cown et al., 1991), and is related to variation in air temperature and the ratio of carbon to nitrogen in the soil (Beets et al., 2007a). While there was no strong evidence of bias across the range of wood density values for the best model developed in this study, there was some indication of an underestimate of AGB at high wood density and biomass data were only available for a sub-set of the latitudinal range on which radiata pine is grown in New Zealand. In particular, data were lacking from older trees grown in the southern parts of the South Island, i.e., Canterbury, Otago and Southland, where wood density is lower. Therefore, if the equations developed in this paper were applied to these stands, they could potentially over-predict above-ground biomass. For example, when the equations developed by Madgwick (1983) were applied to an independent dataset from Victoria, Australia they were found to underestimate stem biomass by 19% (Baker et al., 1984). One reason given for this under-prediction was the high wood density of the Australian trees (Baker et al., 1983). It may be possible to capture some of the variation in the factors affecting wood density through the inclusion of appropriate site-level variables such as latitude, elevation and site index in biomass equations. Elevation and site index were included in the equations developed for Norway spruce branch and needle biomass by Wirth et al. (2004) and a similar approach could be taken for radiata pine; however, more information is required about the characteristics of those sites for which data are available and additional biomass data are required from a broader range of sites. In addition, explicitly representing these factors in national level biomass functions could reduce the potential bias that may occur when applying these functions to individual stands or regions.

The functions developed in this paper do not explicitly account for the effect of live crown pruning on total above-ground biomass. Approximately, 62 percent of the radiata pine forest estate in New Zealand is, or is expected to be, pruned to a height of at least four metres (Ministry of Agriculture and Forestry, 2007). Therefore, ignoring the effect of green crown pruning when estimating the total above-ground biomass of young radiata pine trees could result in above-ground biomass being significantly overestimated as the ratio of crown biomass to total above-ground biomass is relatively high (Cromer et al., 1985; Madgwick et al., 1977). In order to explicitly account for pruning, functions need to be developed to predict foliage and branch mass that include some measure of crown size, e.g. crown length. Using this approach, separate functions would be developed for each component of total above-ground biomass and the predicted values from these functions combined to give the predicted total for the entire tree. Ideally, these equations would be additive so that the sum of predicted biomass for tree components would equal the predicted total for the whole tree (Parresol, 1999).

Although estimates of below ground biomass are often required for estimating total live biomass, only above-ground biomass was considered in the analysis reported here. One of the main reasons for this was that data on root biomass were only available for 31 of the 637 trees in the dataset. However, it is possible to estimate below ground biomass using previously developed allometric functions (Jackson and Chittenden, 1981) or through the ratio of above-ground to below-ground biomass (Beets et al., 2007b).

The models developed in this paper require validation against independent datasets. As is often the case in studies in which biomass functions have been developed for different species, full use was made of all available data. While it is possible to evaluate the accuracy of biomass equations by various cross-validation techniques, this was not done here. The priority for further work is to incorporate data from other studies that have been conducted into the biomass database. Many of the papers that present stand-level biomass estimates describe the measurement of biomass for individual trees that was carried out in order to develop these stand-level estimates. It would also be worthwhile to combine the dataset of measurements from Australian studies that was used in the analysis performed by Snowdon et al. (2000). This would enable an even broader meta-analysis to be performed that would allow the possible significance of site and stand level factors in biomass equations to be tested.

The predictions from the biomass equations presented in this paper should also be compared with those obtained from the C_change model developed by Beets et al. (1999). This model predicts the amount of carbon in the stem, crown, roots and forest floor by using functions for partitioning tree growth into different components (i.e., stem, bark, branches and foliage) coupled with functions that model the mortality and decay of these components. This approach is being used to predict the carbon stocks and stock changes in New Zealand’s planted forests as part of the Land Use and Carbon Analysis System (LUCAS), which is used to fulfil New Zealand’s reporting requirements under the Kyoto Protocol (Stephens et al., 2008). It is therefore important that the two approaches give consistent results when they are used to calculate carbon storage in forests.

Acknowledgments

I would like to thank those researchers who willingly contributed their data which made it possible to develop the functions contained in this report. Drs Herb Madgwick and Craig Trotter, Richard Woollons and Mark Kimberley provided numerous useful advice in the preparation of this report. Funding for this project was provided by the New Zealand Ministry for the Environment under contract 03/04-0271-L and the New Zealand Ministry of Agriculture and Forestry under contract MAF POL 0809-19.

References

All Tables

All Figures

	Figure 1 Distribution of (a) diameter at breast height and (b) total height of the 637 radiata pine trees used to develop the biomass equations.
In the text

	Figure 2 Relationship between total above ground biomass (AGB) for an individual tree and tree height (H), diameter at breast height squared (D2) and D2H.
In the text

	Figure 3 Relationship between total above ground biomass (AGB) for an individual tree and tree height (H), diameter at breast height squared (D2) and D2H on a logarithmic scale for x and y axes.
In the text

	Figure 4 Comparison of residuals against predicted values of above ground biomass (on the logarithmic scale) for Equation (3) which contains diameter at breast height as the sole predictor. The solid line represents a locally-weighted smoothing function (lowess) fitted to the data.
In the text

	Figure 5 Comparison of measured above-ground biomass with that predicted by (a) Equation (18) and (b) the same Equation but with a higher order diameter term (D3) added.
In the text

	Figure 6 Comparison of residuals from the model given by Equation (7) against tree age, stand density and wood density.
In the text