Modelling between tree and longitudinal variation in green density within Pinus radiata: implications for estimation of MOE by acoustic methods

Instruments based on resonance are widely used in the forest industry to predict modulus of elasticity (MOE) and segregate logs of varying quality for different end uses for fast growing softwoods such as Pinus radiata D. Don. Predictions of MOE, made using resonance instruments, often assume a constant green density, ρg, of 1,000 kg m–3. However, little research has been done to test the robustness of this assumption. The objective of this research was to describe changes in predictive precision of MOE as ρg is increasingly well characterised. Longitudinal measurements of velocity, V, and ρg taken from eighty 17-year old unthinned P. radiata trees growing at two sites in Chile were used to calculate MOE. Predictions of MOE were then made by substituting measurements of ρg for values predicted by the following models (i) Model 1 - assuming a constant ρg of 1,000 kg m−3, (ii) Model 2 - using the mean tree ρg of 914 kg m−3, (iii) Model 3 - using a model with fixed effects to account for the mean longitudinal variation in ρg, (iv) Model 4 - inclusion of previous terms and random effects to account for tree level variation and (v) Model 5 - inclusion of previous effects (in model 4) and a random quadratic term. Differences in MOE determined from measurements of ρg and the five predictions of ρg were expressed as both a percentage difference, (D) and an absolute percentage difference (Da) to assess precision and bias. At the tree level, values for mean D and Da (in brackets) were −9.9 (10.4)%, −0.459 (5.49)%, −0.262 (4.15)%, −0.045 (0.232)% and −0.0406 (0.189)%, for Models 1, 2, 3, 4 and 5, respectively. At the log level, considerable longitudinal bias in D was evident for Model 1 where over-prediction of MOE was greatest between relative heights of >0.1–0.4, with D reaching maximum values of −33.8% between relative heights of >0.1–0.2. Assuming constant ρg can result in substantial error in estimates of MOE using acoustic instruments particularly when predictions are made at the log level. The mixed effects modelling approach described here demonstrates a useful method for characterising variation in ρg allowing more accurate estimates of MOE to be made using acoustic methods.


Background
Modulus of elasticity (MOE) is a measure of the resistance of wood to deformation under an applied load, and it is a useful indicator of corewood quality in the widely grown plantation softwood Pinus radiata D. Don. Modulus of elasticity is used as a threshold criterion in machine stress-grading of structural timber and is also a key property for determining quality of laminated veneer timber. Plantation-grown Pinus radiata has a relatively poor MOE (Moore 2012) when compared with other internationally traded structural timber species.
Instruments that use sound or stress waves to estimate velocity and MOE have been widely adopted as they are simple, compact and easy to operate. Although the use of sound to measure MOE has been researched for over 50 years (Schultz 1969), it is only recently that the method has seen widespread application (Dowse and Wessels 2013;Wang et al. 2007). The two types of tools available are resonance tools for use on logs and time of flight tools (TOF) that are suitable for use on both standing trees and logs. Tools based on TOF are generally considered to be less accurate than resonance tools (Andrews 2000(Andrews , 2002. Resonance tools are routinely used within the forest industry. Considerable research has been undertaken with these instruments to characterise variation in MOE at a range of scales. These instruments are also widely used by the forest industry during harvesting operations to segregate logs. Resonance tools can be used to screen logs and divert those below the structural-grade threshold to other applications avoiding the expense of processing wood that will not meet final specifications (Tsehaye et al. 2000;Matheson et al. 2002).
Despite the widespread use of resonance methods, there are errors associated with this type of methodology. One of the potentially largest sources of error is the assumption of a constant green density, ρ g , when predicting MOE from velocity. Green density is used in combination with velocity (V), determined from resonance, to predict MOE using the following equation: Typically values for ρ g of 1,000 kg m −3 are assumed and variation in MOE using this assumption have been widely described among sites and among trees (Watt et al. 2009;Lasserre et al. 2004).
Although constant values for ρ g are often assumed, ρ g is highly dependent on the ratio of heartwood to sapwood within a tree. In P. radiata, sapwood is close to saturation and typically values for ρ g are 1,100 kg m −3 (Loe and Mackney 1953;Cown 1999) while heartwood has a moisture content just above the fibre saturation point and low ρ g that is typically about 600 kg m −3 (Cown 1999). The ratio of heartwood to sapwood has been found to be affected primarily by tree age and height and to a lesser extent by site and thinning (Cown et al. 1991;Kininmonth 1991;Wilkes 1991;Cown 1974). In mature P. radiata ρ g is often lower than 1,000 kg m −3 because the percentage of heartwood increases with tree age (Chan et al. 2012). Indeed, whole-stem values as low as 798 kg m −3 have been found in 36-year-old trees (Chan et al. 2012). Although variation in ρ g has been studied (Chan et al. 2012;Cown and McConchie 1982), little research has fully partitioned and modelled variation in ρ g among sites, among trees and/or within trees, or investigated the impact of this variation on estimates of MOE.
In this study longitudinal measurements of velocity and ρ g obtained from P. radiata growing at two contrasting sites were used to calculate MOE. Predictions of MOE were then made by substituting measurements of ρ g for values predicted by five different models. These models ranged in complexity from the simplest which assumed constant ρ g of 1,000 kg m −3 through to the most complex which was a mixed-effects model that included fixed and random effects to determine ρ g . Differences between measured MOE and that predicted by the five models was determined to quantify the magnitude of errors that occur as ρ g is increasingly well characterised. Finally, a methodology is described for accurately predicting longitudinal variation in ρ g for a single tree from measurements obtained from that tree at a given height.

Study sites
Data for this research were obtained from two sites spanning a climatic range located at Santa Lucia (latitude 37°12 41, longitude 71°48' 31") and Santa Isabel (latitude 37°24 45, longitude 72°15' 29") in the Bio region (region VIII) of Chile. Santa Lucia is located in the Andean foothills at an elevation of 727 m. It has a relatively low mean annual temperature (10.9°C) and high mean annual rainfall (2,022 mm year −1 ). In contrast, Santa Isabel is located at a lower elevation (170 m) in the Central Valley and has higher mean annual temperature (12.9°C) and lower mean annual rainfall (1,166 mm year −1 ) than Santa Lucia. The soil types underlying the sites were volcanic ash and sandy, respectively, for Santa Isabel and Santa Lucia.

Plant material, treatments and stand history
Plant material at both sites were drawn from a half-sib family named Colicheu produced through the breeding programme of the Chilean company Forestal Mininco. Seedlings and one-year-old cuttings were used as treatments at both sites.
For each treatment, a permanent sample plot (PSP) of 630 m 2 (30 m × 21 m) was established in which there were six rows of 15 plants. At both sites, the initial planting density was 1,428 trees ha −1 (2 × 3.5 m) and the PSPs were not thinned or pruned prior to the measurements.

Measurements
Measurements were made during September 2012 on 17-year-old trees that were planted in 1995. A total of 40 sample trees at each site (20 per plant material treatment) were selected for measurement. These represented all diameter classes. Malformed trees (including leaning trees or those with forked trunks or with large branches) were not selected.
Following the felling of each tree, measurements of height and diameter at breast height (DBH) were taken. Stem slenderness was determined from these measurements as height/DBH. Each tree was cut into 3-m long sections until a diameter of 10 cm was reached. All branches were removed. Discs that were 30-mm thick were taken at the top and base of each section. Green density, ρ g , (kg m −3 ) was determined on each sample disc following a standard protocol (ASTM D2395) from the following equation, where W g is green weight (kg), and V g is green volume (m 3 ) determined using the immersion method. All values of ρ g reported in this paper use values of ρ g averaged over the log from measurements made at the top and base. Velocity was determined for each log without removing bark using a resonance tool (Hitman HM200, Fibre-Gen). Velocity (V) was determined from sample length (l) and the resonant frequency (f) as V = 2 fl.

Description of the modelling approach
Longitudinal measurements of velocity and ρ g were used to calculate MOE from Equation (1). Although these calculated values were not direct measurements of MOE, as they were determined from V and ρ g , they will be referred to hereafter as measured MOE for clarity. Predictions of MOE were then made by substituting measurements of ρ g for values predicted by the following models (i) Model 1 -assuming a constant ρ g of 1,000 kg m −3 , (ii) Model 2 -using the mean ρ g across all trees of 914 kg m −3 , (iii) Model 3 -using a model with only fixed effects to account for the mean longitudinal variation in ρ g , (iv) Model 4 -using a model with the previous fixed effects and random effects to account for variation in intercepts for ρ g between trees and (v) Model 5 -using a model with the previous effects (in Model 4) and a random quadratic term to account for differences in the longitudinal pattern of variation in ρ g among trees. A summary of the terms included within each model is given in Table 1.

Model development
Using SAS (SAS Institute Inc. 2008) a random-coefficients mixed-effects model was used to characterise variation in ρ g at the log level for Models 3, 4 and 5. Mixed effects models can include both fixed and random effects. An effect is termed fixed if the levels in the study include all possible levels of the factor, or at least the levels about which inference is to be made. In contrast, a factor is considered random if its levels plausibly represent a larger population with a probability distribution. Only fixed effects (relative height) were used for predictions of ρ g made by Model 3 while fixed effects and the random intercept were included for Model 4 ( Table 1). The full random-coefficients model (Model 5) included fixed effects and random effects for both the intercept and quadratic terms (Table 1).
Relative height, which was defined as the ratio of the measurement or prediction height to the total tree height, was used as a predictive variable within the model rather than actual measurement height. This was because the Akaike Information criterion, which measures the relative quality of a statistical model (Littell et al. 2006), was lower (4673.4 vs. 4725.6) for the final model with the relative height term included. The fixed-effects part of the model included relative height as a third-order polynomial. Neither site nor treatment, nor the interaction of these two terms with any of the terms in the relative height equation, was significant at P = 0.05. Random effects included the intercept (P < 0.001), relative height (P < 0.001) and the square of relative height (P < 0.001). A full description of Model 5 is given in the next section of the methods.
Predictions of ρ g against relative height using only the fixed effects in the model are shown as the black line while predictions of ρ g using all random effects for individual trees (Model 5) are shown as the blue lines in Figure 1. Predictions of longitudinal variation in ρ g for individual trees were markedly offset from the mean value. These results supported the significance of different intercepts for different trees but followed a similar longitudinal pattern to that of the mean, for most trees, although there were exceptions ( Figure 1). Terms included in the model

Analyses
Differences between measured values (y m ) and predictions (y p ) of ρ g and MOE using the five models were expressed at the mean tree level (by averaging across logs) and log level. These differences were expressed as both a percentage difference, D (D = ((y m -y p )/y m ) · 100) and an absolute percentage difference, D a (D a = |D|). For a given model, both D and|D|were identical between ρ g and MOE . Consequently results showing these differences from the measured values can be interpreted to represent either wood property. The distribution and mean values for D and D a were examined for each model to respectively assess model bias and precision. The relationship between relative height and both D and D a was also examined to determine within tree bias and precision for the five models.
Calibration of Model 5 to predict longitudinal variation in ρ g from a single disc The methodology by which the most complete and accurate model (Model 5) can be used to predict longitudinal variation in ρ g from measurements obtained from a single disc is demonstrated below. Use of the linear mixed-effects model described in Model 5 permitted the estimation of a mean response (population-specific) or a calibrated response (cluster-specific) for a new tree (Verbeke and Molenberghs 2000). A mean response can be obtained using only the fixed-effect components of Model 5, where the vector of random effects u k for a new kth individual is assumed to have expected value E(u k ) = 0. In contrast, a calibrated response can be obtained when auxiliary information is available allowing for prediction of random parameter components. In this case, auxiliary information corresponds to a measure or prediction of ρ g at a given stem height for a new tree. Using this information, the vector of random parameters can be predicted using an approximate Bayes estimator of u k (Vonesh and Chinchilli 1997;Rencher 2000).
where W k ¼Ẑ k DẐ 0 k þ R k Â Ã −1 and the error vector The design matrixẐ k and the variance-covariance matrix for the random-effects D and residuals R k are defined for Model 5 as: : : and R k ¼ Var e k ½ ¼ σ 2 I n k , where I n k is an identity matrix of dimension (n k x n k ) and n k is the number of observations used for calibration. Under this structure, random errors are assumed to be uncorrelated and have constant variance (σ 2 ).
Longitudinal predictions of ρ g can be made from Model 5 using ρ g data obtained from a single disc. As an example, ρ g measured at breast height (1.3 m) for a randomly selected sample tree in the dataset was 985.08 kg m −3 . The relative height of the disc was 0.041 (disc height of 1.3 m divided by the total tree height of 31.6 m). Thus, the difference between the observed value and the estimated mean-response value from Model 5 is determined from the fixed effects component of the model as: The design matrixẐ k and the variance-covariance matrices for the random coefficients and error estimated for Model 5 are determined as, Now, replacing the matrices in Equation 3 gives the following predictions for the random parameters of this specific tree: u 0k =49.5454, u 1k = −49.9722 and u 2k = 18.8294  Figure 1 Predictions of green density using Model 5 plotted against relative height. Black line: mean relationship between relative height and green density, predicted from the fixed effects; blue lines: tree-level estimates (predicted using both fixed-and random effects). and the calibrated parameters for the model are given below in Equations 7, 8 and 9: The accuracy of this prediction is assessed by comparing the predicted longitudinal pattern to that obtained through measurements.

Variation in tree dimensions
Tree diameter ranged from 18.5 to 40.1 cm (average 27.6 cm) and was not significantly affected by either the main or interactive effects of treatment and site ( Table 2). The total tree height of cuttings at Santa Lucia was significantly greater than for the other three site × treatment combinations, but no other significant differences in height were noted among treatments (Table 2). Tree height ranged from 20.4 -35.5 m (average 28.0 m). Neither the main nor interactive effects of site and treatment influenced stem slenderness which ranged from 82.0-129 (average 103) ( Table 2).

Variation in green density and MOE
Neither site, treatment nor their interaction significantly affected tree level values of ρ g (Table 2). At the tree level, ρ g ranged from 792-1,018 kg m −3 with a mean value of 914 kg m −3 (Table 2). Within-tree variation in green density was also large and showed a complex pattern when plotted against relative height (Figure 2a). Green density declined rapidly from mean values of 922 kg m −3 at relative heights of 0 -0.1 reaching mean minima of 861 kg m −3 at relative heights of >0.1 -0.2 before increasing slowly to mean maxima of 1,006 kg m −3 at relative heights of >0.7-0.8 (Figure 2a, Table 3).
At the tree level, there were marginally significant (P = 0.04, Table 2) treatment differences in MOE. Although mean values for seedlings exceeded those of cuttings by 4% (least square means of 8.05 vs. 7.73 GPa), these differences were only significant at Santa Lucia (8.07 vs. 7.62 GPa). Neither site nor the interaction of site and treatment significantly affected MOE (Table 2).
At the tree level, MOE ranged from 6.29-9.83 GPa (average 7.89 GPa, Table 3). The within-tree pattern in MOE was the opposite to that of ρ g (Figure 2b). Modulus of elasticity increased from mean minimum values of 6.76 GPa at relative heights of 0 -0.1 reaching highest values at relative heights between >0.1 -0.4, with a mean maximum of 8.57 GPa reached at relative heights of >0.2 -0.3 (Table 3). Values of MOE declined at relative heights above 0.4 with a mean value of 6.83 GPa at relative heights ranging from >0.7 -0.8 (Figure 2b, Table 3).

Variation in green density and MOE from actual values for the five models tested Tree-level variation
Estimates of ρ g and MOE, at the tree level were considerably biased using Model 1. Using an assumed ρ g of 1,000 kg m −3 resulted in over-prediction of both ρ g and MOE for all but two of the 80 trees (Figure 3a). Mean D was −9.9% with range of −26.6 to 1.6% (Figure 3a) while mean D a was 10.4% with range of 3.78 to 26.6% (Figure 3b) for both properties.
Use of Model 2 to predict ρ g improved both bias and precision at the tree level for both wood properties. Mean D was −0.459% with range of −15.7 to 10.1% Values shown include the tree-level mean followed by the standard error in brackets 1 For variables where there was a significant effect, means followed by different letters significantly vary at P = 0.05. 2 F-values are shown followed by the P category, where asterisks **, *denote respective significance at P = 0.01 and 0.05; ns is not significant at P = 0.05.
( Figure 3a) while mean D a was 5.59%, with range of 2.00 to 15.7% (Figure 3b). Little reduction in precision was achieved at the tree level using Model 3, for which mean D was −0.262% with range of −15.2 to 9.74% (Figure 3a) while mean D a was 4.15%, with range of 0.968 to 15.2%. Accounting for tree-level variation in the intercept (Model 4) and variation in the longitudinal pattern (Model 5) for ρ g resulted in marked improvements over Models 1, 2 and 3 for both properties at the tree level. Precision and bias was very similar between Models 4 and 5, and Model 5 was the least biased and most precise of all models tested. For Model 5, the mean D was −0.0406%, with range of −1.27 to 0.477% (Figure 3a) while mean D a was 0.189% with range of 0.0088 to 1.27% (Figure 3b).

Within-tree variation
Using Model 1, over-prediction of MOE was most marked in the lower half of the stem reaching maximum over-prediction between relative heights of >0.1-0.2 (Figures 4a,f ). Maximum D a was 33.8% between these relative heights (Figure 4f ) while mean D a was 16.5% for both properties (Figure 5a).
In addition to over-predicting measured values, the assumption of 1,000 kg m −3 for ρ g also resulted in prediction errors in the longitudinal pattern in MOE. For Model 1, there was a sharp increase in predicted MOE with relative height that declined markedly after relative heights of >0.1-0.3 were reached (Figure 5b). In contrast, measured MOE increased more slowly with relative height, with highest values found across a broader range of relative heights from >0.1-0.4 (Figure 5b). The largest difference between measured values of MOE and that predicted by Model 1 occurred at the highest values of MOE at relative heights of between >0.1-0.4 (Figure 5b).
Estimates of MOE using predictions of ρ g from Model 3 were relatively unbiased and precise (Figures 4c,h). Predictions from Model 3 closely followed measured changes in ρ g (Figure 5c) and MOE (Figure 5b) across the relative height range. Consequently, across the relative height range, there was little bias in either D (Figure 4h) or D a (Figure 5a).
Inclusion of random effects to account for tree-level influences on ρ g resulted in substantial gains in precision over the other three models. Estimates using either Model 4 or Model 5 resulted in little apparent bias when (h) Figure 4 Relationships between relative height and the difference (left hand panels) and percentage difference (right hand panels) for MOE predicted using Model 1 (a, f), 2 (b, g) 3 (c, h), 4 (d, i) or 5 (e, j).
predictions were plotted against relative height for either ρ g (Figure 5c) or MOE (Figures 4d,e,i,j, 5b). Across the range in relative heights, mean D a was 2.11% for Model 4 and 1.44% for Model 5.

Calibration of predictions made using Model 5
Longitudinal variation in green density for a single sample tree was determined using either the mean response or the calibrated response for Model 5. Both these results are presented in Figure 6 along with measured values. Clearly, the calibrated response accurately represents longitudinal variation in ρ g ( Figure 6) and represents a substantial improvement in precision over mean-response predictions from Model 5 (which included fixed effects only).

Discussion
This study shows that assuming a constant value of ρ g can result in substantial errors in estimates of MOE using acoustic instruments. Differences from the assumed value of ρ g for the site, and significant variation in ρ g among and within trees were three sources of variation in ρ g found to substantially affect predictions of MOE. The modelling approach described here demonstrates a useful Modulus of elasticity (GPa) Relationship between relative height and measured green density (open black circles) for a tree randomly selected from the dataset. Longitudinal variation in green density generated from a single measured disc of green density (at a relative height of 0.014) using Model 5 (red squares) and predictions of green density using Model 5 containing only fixed effects (blue triangles) are also shown. Table 3 Variation in green density and modulus of elasticity (determined using values of green density) with relative height determined from 80 trees at the two sites method for characterising variation in ρ g so that more accurate estimates of MOE can be made using acoustic methods. Assuming a constant ρ g of 1,000 kg m −3 is likely to result not only in poor precision but also biased estimates of MOE. Site values of ρ g were considerably lower than 1,000 kg m −3 so assuming this value resulted in consistent over-prediction of both ρ g and MOE at the tree level. Previous research has shown that there is an increasing proportion of heartwood as P. radiata matures (Chan et al. 2012). Consequently, recorded values of whole stem ρ g in mature stands are often lower than 1,000 kg m −3 , as demonstrated in the present study and in previous research (Waghorn et al. 2007;Chan et al. 2012;Cown and McConchie 1982).
Assuming a constant ρ g of 1,000 kg m −3 also resulted in poor prediction of longitudinal variation in MOE. Measured ρ g exhibited significant longitudinal variation that was consistent with previously observed longitudinal patterns (Chan et al. 2012;Cown and McConchie 1982). As ρ g was lowest between relative heights of >0.1 − 0.4, results show MOE to be markedly over predicted for the lower half of the tree when a ρ g of 1,000 kg m −3 was assumed. In addition to over-predicting values of MOE, this assumption also resulted in erroneous prediction of the longitudinal pattern in MOE. The highest values for MOE were reached over a broader range of tree heights (relative heights of >0.1 − 0.3) than the more peaked pattern described using a constant ρ g of 1,000 kg m −3 .
At the tree level, assumption of a site mean value for ρ g (Model 2) resulted in a significant precision gain over the model that assumed ρ g of 1,000 kg m −3 (Model 1). This assumption effectively removed bias at the tree level and reduced the mean absolute error from 10.4% to 5.59%. Despite these improvements, there was still considerable error at the tree level that was attributable to among-tree variation in ρ g . Predictions using Model 2 also showed considerable within-tree bias that was attributable to the large longitudinal variation in ρ g not captured by site mean level values in ρ g .
Although the site mean value of ρ g was markedly lower than 1,000 kg m −3 , neither site nor treatment significantly affected values of ρ g . The lack of a significant site effect on ρ g in this study is likely to be attributable to the low number of sites sampled. Previous research has shown significant site variation in ρ g and also that ρ g is lowest in P. radiata stands growing in drier areas (Chan et al. 2012). For trees growing under drought conditions, whole stem values for 36-year-old trees as low as 798 kg m −3 have been recorded (Chan et al. 2012). Further research is needed to develop models that can characterise variation in site level ρ g across gradients in age and environment.
At the tree level, modest gains in predicting MOE were made by accounting for the longitudinal pattern in density. As the longitudinal pattern was relatively similar between trees, this did remove the within-tree bias that was evident in log-level predictions. The longitudinal variation in ρ g reported here is consistent with previous research within mature P. radiata stems (Chan et al. 2012;Cown and McConchie 1982). These authors also found a decline from the tree base to about one-third tree height, above which height values for ρ g increased reaching maxima within the top third of the tree. Although this pattern remains the same regardless of tree age, the values for the longitudinal minima decline as trees age. A study within New Zealand across a range of stand ages found minimum ρ g at ca. 30% of total tree height, with values for ρ g ranging from ca. 985 kg m 3 at stand age 12 to ca. 770 kg m −3 at stand age 52 (Cown and McConchie 1982).
The marked between tree variation in ρ g was well accounted for by addition of random terms to the model. Addition of random effects to Models 4 and 5 substantially improved predictions of ρ g compared with the model with fixed effects only (Model 3). The parameter D a declined by over an order of magnitude from 4.15% for Model 3 to 0.231% for Model 4 and 0.189% for Model 5. Most of the variation among trees was accounted for by addition of a random term allowing for different intercepts in the longitudinal pattern (Model 4). Although gains were significant, addition of a random term that allowed for variation in the longitudinal shape did not greatly improve D a . The relative importance of shape and intercept is visually evident through inspection of Figure 1 which shows a wide range in different intercepts between trees but somewhat less variation in the shape.
Clearly, use of a mixed-effects model that includes random terms to account for tree level variation, such as Model 5, provides the most accurate means of modelling longitudinal variation in ρ g and MOE between sites and trees. A methodology for parameterising such a model from measurements of ρ g at a single height has been demonstrated. This procedure could considerably reduce the sampling effort required to characterise longitudinal variation in ρ g . Further research should explore whether easier-to-obtain, non-destructive measurements (such as increment cores) could be used to parameterise this type of model.

Conclusion
Results from this study show that marked over-prediction of ρ g and MOE occurred when constant ρ g values of 1,000 kg m −3 were assumed. At the log level the most accurate predictions of ρ g were obtained using a mixedeffects model that included fixed effects to account for mean longitudinal variation in ρ g and random effects to account for deviations in this longitudinal variation between trees. Using this mixed effects model, log level D a