Acoustic Velocity — Wood Fiber Attribute Relationships for Jack Pine Logs and Their Potential Utility

This study presents an acoustic-based predictive modeling framework for estimating a suite of wood fiber attributes within jack pine (Pinus banksiana Lamb.) logs for informing in-forest log-segregation decision-making. Specifically, the relationships between acoustic velocity (longitudinal stress wave velocity; vl) and the dynamic modulus of elasticity (me), wood density (wd), microfibril angle (ma), tracheid wall thickness (wt), tracheid radial and tangential diameters (dr and dt, respectively), fiber coarseness (co), and specific surface area (sa), were parameterized deploying hierarchical mixed-effects model specifications and subsequently evaluated on their resultant goodness-of-fit, lack-of-fit, and predictive precision. Procedurally, the data acquisition phase involved: (1) randomly selecting 61 semi-mature sample trees within ten variable-sized plots established in unthinned and thinned compartments of four natural-origin stands situated in the central portion the Canadian Boreal Forest Region; (2) felling and sectioning each sample tree into four equal-length logs and obtaining twice-replicate vl measurements at the bottom and top cross-sectional faces of each log (n = 4) from which a log-specific mean vl value was calculated; and (3) sectioning each log at its midpoint and obtaining a cross-sectional sample disk from which a 2 × 2 cm bark-to-pith radial xylem sample was extracted and subsequently processed via SilviScan-3 to derive annual-ring-specific attribute values. The analytical phase involved: (1) stratifying the resultant attribute—acoustic velocity observational pairs for the 243 sample logs into approximately equal-sized calibration and validation data subsets; (2) parameterizing the attribute—acoustic relationships employing mixed-effects hierarchical linear regression specifications using the calibration data subset; and (3) evaluating the resultant models using the validation data subset via the deployment of suite of statistical-based metrics pertinent to the evaluation of the underlying assumptions and predictive performance. The results indicated that apart from tracheid diameters (dr and dt), the regression models were significant (p ≤ 0.05) and unbiased predictors which adhered to the underlying parameterization assumptions. However, the relationships varied widely in terms of explanatory power (index-of-fit ranking: wt (0.53) > me > sa > co > wd >> ma (0.08)) and predictive ability (sa > wt > wd > co >> me >>> ma). Likewise, based on simulations where an acoustic-based wd estimate is used as a surrogate measure for a Silviscan-equivalent value for a newly sampled log, predictive ability also varied by attribute: 95% of all future predictions for sa, wt, co, me, and ma would be within ±12%, ±14%, ±15%, ±27%, and ±55% and of the true values, respectively. Both the limitations and potential utility of these predictive relationships for use in log-segregation decision-making, are discussed. Future research initiatives, consisting of identifying and controlling extraneous sources of variation on acoustic velocity and establishing attribute-specific end-product-based design specifications, would be conducive to advancing the acoustic approach in boreal forest management. Forests 2018, 9, 749; doi:10.3390/f9120749 www.mdpi.com/journal/forests Forests 2018, 9, 749 2 of 28


Introduction
The Canadian forest sector has been increasingly embracing an aspirational trivariate management proposition for which the goal is to maximize end-product value while simultaneously realizing volumetric yield and ecological sustainability objectives (sensu [1,2]).Tangential to this transformative shift is the increased apprehension regarding the quality and associated end-product potential of the increasing wood volumes harvested from more intensely managed second-growth forests [3].Among the many prerequisites required for realizing this goal and addressing knowledge gaps in second-growth fiber quality, is the provision of enhanced in-forest intelligence in regards to the forecasting of end-product potential at the time of harvest or immediately afterwards (n., in-forest generically refers to all forest operations conducted or decisions made before mill-processing commences inclusive of operational felling areas, landings and log sorting yards).
End-product forecasts can be used to inform and optimize log-segregation, allocation and merchandizing decision-making thus contributing to the potential realization of the value maximization objective [4][5][6].More precisely, the end-product potential of the logs extracted from the main stem of individual coniferous trees is a function of (1) external morphological stem features (e.g., diameter, height, taper, sweep, crook, and eccentricity), and (2) internal anatomical characteristics (e.g., microfibril angle, tracheid wall thickness, tracheid radial and tangential diameters) and associated physical properties (e.g., modulus of elasticity, wood density, fiber coarseness and specific surface area) of the xylem tissue.Various performance measures which are used to define and classify the overall type and grade of derived products have been shown to be explicitly linked to these internal fiber-based attributes.Although only stiffness and wood density measures are currently used in machine stress grading systems for solid wood products (e.g., dimensional lumber [7]), microfibril angle, tracheid cell wall thickness, radial and tangential tracheid diameters, fiber coarseness and specific surface area, are also important determinates of end-product quantity and quality [8] (sensu Table 1).Thus, as the fiber supply increases in diversity with the arrival of wood from density manipulated forests (e.g., pre-commercially and/or commercially thinned natural-origin stands, and plantations) combined with increasing pressures to find economic efficiencies within the upper portion of the forest products supply chain, provision of explicit information on these additional attributes may become consequential.
Forests 2018, 9, 749 3 of 28 Attribute determination for use in forecasting end-product potential deploying destructive sampling approaches can be logistically challenging, time consuming and expensive (e.g., in-forest extracting of 12 mm diameter xylem cores followed by Silviscan-based laboratory processing and testing).Conversely, non-destructive sampling methods in which attributes can be estimated indirectly through their relationship with acoustic velocity has been shown to be viable alternative [9][10][11][12].More specifically, in-forest acoustic-based sampling tools provide forest practitioners with an ability to estimate the dynamic modulus of elasticity within harvested logs and standing trees.The dynamic modulus of elasticity is a wood stiffness measure that is used as a surrogate measure of its static counterpart (static modulus of elasticity).This latter static measure is commonly employed along with maximum knot size, sweep and physical dimensions, to categorize solid wood products into specific grade classes [7,9,13].Thus, the ability to forecast one of the key underlying determinates of end-product potential within the forest at the time of harvest or soon after, has important utility in terms of informing segregation, allocation and merchandizing decision-making [6].
Mathematically, as derived from engineering principles, the dynamic modulus of elasticity (MOE dyn denoted m e in this study; GPa) of a log can be expressed as a function of the density-weighted velocity of a longitudinal stress wave (v l ; km/s) that propagates through the xylem tissue upon its generation arising from a mechanically induced impact on one of the open cross-sectional log faces (Equation (1) and denoted the primary relationship herewith; [12]).
where w d is the green (fresh) wood density (kg/m 3 ) of the log's xylem tissue.The dynamic modulus of elasticity has also been shown to be correlated with other attributes underlying end-product potential [14][15][16].
More specifically, based on extracted clear xylem samples obtained from black spruce (Picea mariana (Mill) B.S.P.), red pine (Pinus resinosa Ait.) and jack pine (Pinus banksiana Lamb.) trees growing in the central North America, Silviscan-derived estimates of wood density (oven-dried wood density denoted w d in this study; kg/m 3 ), microfibril angle (MFA denoted m a in this study; • ), tracheid cell wall thickness (w t ; µm), radial (d r ; µm) and tangential (d t ; µm) tracheid diameters, fiber coarseness (c o ; µg/m) and specific surface area (s a ; m 2 /kg), have been shown to be correlated with the dynamic modulus of elasticity (Table 2).Statistically significant (p ≤ 0.05) (1) directly proportional relationships between m e and w d , w t and c o (m e ∝ w d , w t , c o ), and (2) indirectly proportional relationships between m e and m a , d r , d t and s a m e ∝ m −1 a , d −1 r , d −1 t , s −1 a , have been established for these species (Table 2).Supporting correlative evidence has also been reported for non-boreal conifers.For example, a directly proportional correlative relationship between m e and wood density (r = 0.75) and an inversely proportional relationship between m e and m a (r = −0.83)has been reported for Silviscan-derived attribute estimates obtained from clear wood samples extracted from 76 radiata pine (Pinus radiata D. Don.) boards [13].Based on this collective correlative evidence, an empirical-based expansion of the acoustic-based inferential framework was proposed for boreal conifer logs [15], as summarized in Table 2. Since wood density is also one of the principal end-product determinates (sensu Table 1) and analytically required for acoustic-based attribute predictions (e.g., Equation (1); Table 2), a simplified directly proportional relationship between acoustic velocity and wood density was included within this framework (i.e., w d ∝ v 2 l ).Collectively, this expanded acoustic inferential framework represents a more encompassing analytical platform for forecasting commercially relevant wood quality attributes than previous univariate approaches solely based on the primary relationship (Equation (1)).a Species-specific bivariate linear associations between Silviscan-determined accumulative area-weighted mean fiber attributes as previously presented in Table 2 of Newton [16]: Pearson product moment correlation coefficient for the relationship between the dynamic modulus of elasticity and each attribute at breast-height (1.3 m) for 50 semi-mature black spruce trees [14], 54 mature red pine trees [16], and 61 semi-mature jack pine trees (this study) where * denotes a significant correlation at the 0.05 probability level.b Based on the primary m e ∝ w d v 2 l relationship for all attributes except wood density.
Forests 2018, 9, 749 5 of 28 Although these statistical-based secondary relationships lack the causal underpinnings of the primary relationship (sensu [17]), the empirical results to date have been supportive.Specifically, applying this analytical framework to red pine logs, revealed that statistically viable acoustic-based relationships could be established for five of the eight attributes examined (m e , w d , w t , c o and s a ) [15].However, the applicability of this expanded framework of acoustic-based wood attribute relationships to other intensely managed softwood species is largely unknown.Furthermore, even the most fundamental acoustic velocity-attribute relationship (Equation (1)) has yet to be parameterized for most of the commercially important Canadian boreal conifers (e.g., an "acoustic velocity" keyword search for articles published within Canada's primary applied journal (Forestry Chronicle) was unsuccessful in identifying previous research in this field).Conversely, other forest regions, such as New Zealand, United States of America, Australia, Europe, and Asia, have adapted this acoustic-based innovation pathway in order improve segregation decision-making (e.g., [18,19]).Hence, the evaluation of the acoustic approach and the provision of prediction equations for the estimating a suite of key attributes underlying end-product potential for jack pine, could contribute to addressing this innovation void.Principally, by providing the forest practitioner with the prerequisite prediction equations, in-forest segregation decision-making could be improved and potentially yield overall efficiency gains within the upper portion of the jack pine forest products supply chain.Consequently, the objective of this study was to investigate the empirical applicability of the expanded acoustic-based inferential framework by examining the nature, strength and predictability of the relationships between acoustic velocity and the dynamic modulus of elasticity, wood density, microfibril angle, tracheid cell wall thickness, radial and tangential tracheid diameters, fiber coarseness and specific surface area, specifically for jack pine.

Sample Stands, Plot Establishment and Sample Tree Selection
Two geographically separated jack pine thinning experiments that were previously established in northeastern (denoted the Sewell site) and northcentral (denoted the Tyrol site) regions of Ontario, were selected for sampling.The experimental treatments were representative of the crop planning strategy historically followed in Ontario for this species.Specifically, allowing jack pine to naturally regenerate following a stand-replacing disturbance and if sufficiently successful in terms of stocking levels and free-to-grow status, (1) implementing precommercial thinning in order to moderate overall site occupancy and encourage individual tree growth thereby reducing the time-to-operability status in order to mitigate the effects of projected wood supply deficits, or (2) implementing a temporal sequence of thinning treatments (precommercial and commercial thinning) in order to optimize site occupancy and maximize end-product potentials throughout the rotation [20].
At the Sewell site, sample trees were selected within 6 sample plots that were established in 3 fire-origin jack pine stands (2 plots/stand).These approximately 53-year-old stands were situated on sites of medium-to-good quality (mean site index of 18 m at a breast-height age of 50 [21]) and geographically located within Forest Section B.7-Missinaibi-Cabonga of the Canadian Boreal Forest Region [22].The soils and topography were characterized as deep (>1 m) coarse-to-medium textured sandy types situated on gently undulating terrain.Silviculturally, the stands were subjected to 1 of 3 density manipulation treatments yielding 3 unique crop plans: (1) natural stand development with no density management (unthinned); (2) precommercial thinned (PCT) at age 11 (1971); and (3) PCT at age 11 followed by a light pseudo-commercial thinning (CT) at age 43 (2003).Procedurally, for the first stand, sample trees were selected within 2 previously established 0.06 ha circular monitoring plots.For these 2 plots, a stratified random sample protocol was used to select sample trees: the observed diameter frequency distribution at the time of sampling was stratified into 3 size classes from which 1 or 2 trees per class were selected.However, to preserve the original sample design within the remaining stands, temporary prism-based variable-size plots where established adjacent to 2 of the existing long-term monitoring plots within the second and third stands.Deploying the observed diameter frequency distribution derived from the prism sweeps, 3 size classes were similarly formulated from which 1 or 2 trees per class were selected.In total, 31 jack pine sample trees were selected at the Sewell site: 3 stands × 2 plots/stand × 3 size classes/plot × 1-2 trees/size class.
At the Tyrol site, sample trees were selected within 4 0.07 ha circular sample plots that were established in 2 fire-origin jack pine stands.These stands were situated on sites of good-to-excellent quality (mean site index of 21 m at a breast-height age of 50 [21]), geographically located within Forest Section B9 (Superior) of the Canadian Boreal Forest Region [22], and approximately 73 years of age when sampled.The soils and topography were characterized as deep (>1 m) finely textured sandy types on gently undulating terrain.Silviculturally, the stands were subjected to 2 density manipulation treatments: PCT at an age of approximately 20 (1962) followed by a light pseudo-CT treatment at an age of approximately 56 (1998).A size-based stratified random sample protocol was used to select trees resulting in approximately 2-3 trees being chosen from each of the 3 size (diameter) classes.In total 30 jack pine sample trees were selected from the Tyrol site: 2 stands × 2 plots/stand × 3 size classes/plot × 2-3 sample trees/size class.In summary, a total of 61 sample trees were selected from the 2 experimental areas which were grown under a range of density management regimes reflective of past and current forest management practices for this species and region (e.g., unthinned, PCT, PCT+CT).

Sample Tree Measurements, Stem Analysis Procedures, Log-Based Acoustic Measurements and Disk Sampling
Prior to felling and sectioning each sample tree using destructive stem analysis, diameter at breast-height (1.3 m) outside-bark diameter (D b ; cm), total height (H t ; m) and height-to-live crown (H c ; m) measurements were obtained.Note, all measurements and destructive sampling were completed on the trees at Sewell and Tyrol sites at the conclusion of the 2013 and 2015 growing seasons, respectively.The destructive stem analysis protocol consisted of felling each sample tree at stump height (0.3 m), delimbing the stem and then topping the stem at an 80% relative height position.The stem was then sectioned into 0%-20%, 20%-40%, 40%-60% and 60%-80% percentile-based log-length intervals employing a percent height sampling protocol.Immediately thereafter, the velocity of a mechanically induced longitudinal stress wave propagating throughout each log was measured using a Director HM200 acoustic velocity resonance tool (Fibre-gen Inc., Christchurch, New Zealand; www.fibre-gen.com).Specifically, a twice-replicated set of v l (km/s) measurements were obtained at the bottom and top of each log from which an arithmetic mean value was calculated.Log lengths, ambient air, and xylem temperatures were also taken at the time of each acoustic velocity measurement.
For each tree, cross-sectional samples were then extracted from the center of each log, corresponding to the relative height positions of 10%, 30%, 50% and 70%.The samples were placed in temporary cold storage within 8 h of sectioning and then transported to a permanent cold storage (<0 • C) facility until further processing.Overall, the deployed sampling design yielded 124 cross-sectional disks from the Sewell site (1 disk/log × 4 logs/tree × 31 trees) and 120 cross-sectional disks from the Tyrol site (1 disk/log × 4 logs/tree × 30 trees).Due to logistical challenges during field sampling at the Sewell site, however, one of the cross-sectional samples had to be disregarded.Hence a final total of 243 cross-sectional disk samples obtained from the midpoint position on each of the 243 sample logs were available for analysis.Table 3 summarizes the mensuration characteristics of the standing sample trees prior to felling.The trees from both sites exhibited a similar degree of variation for each measured characteristic as quantified by the coefficient of variation.Diameter and live crown ratio ((H t −H c )/H t ) exhibited the greatest amount of variation given that a size-based stratified random sampling protocol was used to select the sample trees (e.g., trees from all 3 diameter size-groups (small, medium, and large) where included).The minimal degree of age variation is largely due to the even-aged nature of the sampled stands.Table 4 provides a descriptive summary of the derived logs in terms of their dimensions and associated acoustic velocity measurements.

Silviscan-3 Estimation of Fiber Attributes and Preliminary Computations
A transverse 2 cm × 2 cm bark-to-pith-to-bark sample was extracted along the geometric mean diameter of each of the frozen cross-sectional disks.Annual area-weighted ring measures of m e , w d , m a , w t , d r , d t , c o and s a were determined along one of the transverse pith-to-bark radii using the SilviScan-3 system.SilviScan-3 analyses yields a set of commercially relevant fiber attribute estimates through the deployment of an integrated hardware-software processing system involving automatic image acquisition and analysis (cell scanner), X-ray densitometry, and X-ray diffractometry technologies.In this study, the following attribute estimates were used from the Silviscan-3 analysis: (1) wood density as measured at a 25 µm resolution and 8% moisture content (dry basis) and determined following the removal of resins using X-ray densitometry (i.e., w d = f (intensity of incident and transmitted X-ray beams; sample thickness; and X-ray travel distance); [23,24]); (2) microfibril angle as determined via X-ray diffraction patterns (i.e., m a = f (variance of the cellulose-I 002 azimuthal X-ray diffraction pattern; and dispersion of the microfibril orientation distribution); [25]); (3) dynamic modulus of elasticity as determined from a combination of X-ray densitometry and diffraction measurements (i.e., m e = f (X-ray densitometry density estimate; and coefficient of variation of the normalized intensity profile); [26]); and (4) fiber dimensions and derived metrics consisting of tracheid wall thickness, radial and tangential tracheid diameters, fiber coarseness and specific surface area as determined using the results from the image analyses (radial and tangential diameters) in combination with the wood density estimate (sensu [27]).Computationally, deploying the annual area-weighted ring estimates for m e , w d , m a , w t , d r , d t , c o and s a and proceeding in a pith-to-bark direction, an attribute-specific cumulative area-weighted moving average value was calculated for each of the 243 radial sequences (Equation ( 2)).
where F is the cumulative moving area-weighted fiber attribute average for each sequence, a i is the Silviscan-3-determined area of the ith annual ring (mm 2 ; i = 1, . . .,I; I = total cambial age), and f i is the mean Silviscan-3-determined fiber attribute value within the ith annual ring.Table 5 provides an attribute-specific descriptive statistical summary by log position, inclusive of measures of central tendency (arithmetic means and medians), ranges (minimums and maximums) and relative variation measures (coefficients of variation).

Preliminary Data Stratification for Cross-Validation Assessment
To minimize the potential effects of spatial correlation arising from the selection of the multiple logs from the same sample tree while at the same time providing a validation data subset, the full data set consisting of 243 m e , w d , m a , w t , d r , d t , c o , s a -v l observational pairs was systematically subdivided into calibration and validation data subsets of approximately equal size.This involved randomly selecting from each sample tree a sequential set of measurements consisting of either those associated with the 1st and 3th positioned logs or those associated with the 2nd and 4th positioned logs.The resultant pairs were then allocated to either the calibration or validation subsets.At the conclusion of this stratification process, the calibration and validation data subsets were comprised of 122 and 121 m e , w d , m a , w t , d r , d t , c o , s a -v l observational pairs, respectively.As presented in Table 6, the acoustic velocity measurements and attribute values within the subsets were very similar in terms of measures of central tendency, magnitudes, and relative variation.A further similarly structured comparative assessment of the data subsets by log position and location (not shown) also confirmed the equivalency between the data subsets.

Specifying Functional Forms and Parameterization Methods Utilized
Deploying the expanded acoustic-based inferential framework, the Silviscan-derived ( 1) m e , m a , w t , d r , d t , c o and s a estimates were expressed as a function of the Silviscan-derived wood density estimates and HM200-based acoustic velocity measures (m e , m a , w t , d r , d t , c o , s a = f w d v 2 l ), and (2) w d estimates were expressed solely as a function of HM200-based acoustic velocity measures (w d = f v 2 l ).Statistically, the model specification procedure consisted of assessing the results from extensive graphical and correlation analyses to determine the most appropriate functional form for each of these 8 relationships.Employing the full data set, bivariate scatterplots were used to determine the degree of linearity between the each dependent and independent variable (i.e., m e , m a , w t , d r , d t , c o , and s a versus w d v 2 l , and w d versus v 2 l ).A visual interpretation of the graphics for m e , w d , m a , w t , c o , and s a indicated mostly linear trends whereas there were no detectable trends for d r and d t .Consequently, log-linear, log-log and power-based transformed relationships were also examined for d r and d t using both scatterplots and statistical measures of association.However, the graphical trends and related measures of linear association (Pearson moment correlation coefficient) from this supplementary analysis did not reveal the presence of viable linear relationships.
Similar to the approach used for quantifying acoustic-attribute relationships for red pine logs (i.e., [15]), a two-level hierarchical mixed-effects linear model specification consisting of fixed and random effects was used (sensu [28]).Specifically, at the first level, log-specific mean attribute values were expressed as a linear function of the density-weighted or density-unweighted acoustic velocity (Equation ( 3)).
where F (k ) l is the cumulative area-weighted moving average value for the k'th attribute where k = m e , m a , w t , d r , d t , c o or s a for the lth log, F (w d )l is the cumulative area-weighted moving average wood density value for the lth log, v l l is the mean acoustic velocity value for the lth log, β 0 l (k) is a first-level random effects intercept parameter specific to the kth attribute (k = k attributes and w d inclusive) that is allowed to vary across the L logs, β 1(k) are first-level fixed effects slope parameter specific to the kth attribute, and ε (k) l is a random error term specific to the kth attribute for the lth log.The random errors (ε (k) l ) are assumed to be independent, uncorrelated and have constant variance.The second level expressed the first-level intercept parameter as a function of a grand mean plus a random error term whereas the slope parameter was considered a fixed effect (Equations (4a) and (4b), respectively).
where γ 0(k) and γ 1(k) are attribute-specific grand mean values and u 0(k) l is an attribute-specific random error term specific to the lth log.The final mixed-effects model specifications were derived by combining the level one-and two-level expressions into a single model (i.e., Equation ( 5)).
Given the complex error structures, the parameter estimates and model statistics were obtained via the deployment of the hierarchical linear and nonlinear modeling software program HLM7 [29].Statistically, the program provides empirical-Bayes first-level parameter estimates for each randomly varying coefficient (intercept term), generalized least squares-based estimates for second-level terms, and maximum likelihood estimates for the variance and covariance components.Procedurally, employing the observational pairs from the calibration data subset, the models were parameterized and subsequently assessed on their compliance with the underlying statistical assumptions.Based on the protocol developed by Raudenbush and Bryk [28], this evaluation included testing the (1) constant variance assumption among first stage residuals, and (2) presence of significant random effects as determined via testing the null hypothesis (u 0 = 0) versus the alternative hypothesis (u 0 = 0).Park's homogeneity test was used to evaluate the constant variance assumption for each of the significant regression relationships (sensu [30]).This involved regressing the first stage Bayes residual values (logarithmic square values; dependent variable) against the predictive variable values (logarithmic values; independent variable) using a simple linear regression model specification, and then testing the null hypothesis that the resultant slope parameter estimate was not significantly different from zero at the 0.01 probability level.Resultant slope values not significantly different from zero were supportive of the homoscedasticity assumption.Conversely, slope values significantly different from zero are suggestive of the presence of heteroscedasticity (non-constant variance).Additionally, given the potential effect of spatial correlation on statistical inference in terms of producing inefficient estimators without the minimum variance best linear unbiased estimator property when present, the spatial correlation among the empirical-Bayes residuals derived from the first-level models was assessed.The approach consisted of employing residual regressions in order to detect the presence of a first-order spatial correlation scheme among adjacent residuals for logs sampled from the same tree (i.e., either the 1st and 3rd log pair or the 2nd and 4th log pair): i.e., fitting the relationship e l+2 = ρe l + v a where e l and e l+2 are the residual values for the lth and lth + 2 ordinal positioned log, respectively (1st and 3rd ordinal log, or 2nd and 4th ordinal log, respectively), p is the first-order spatial correlation coefficient estimate (−1 ≤ p ≤ +1), and v a is a random error term [30].For significant (p ≤ 0.01) regression relationships, the resultant p estimate reflected the magnitude of the first-order spatial correlation.Otherwise, for regressions that were not significant, spatial correlation was assumed to be absent and the hence the independence assumption was not rejected.The presence of potential outliers or influential observations, systematic lack-of-fits, and non-constant variance was also assessed through the examination of predictor-residual error bivariate scatterplots.
These selected specifications acknowledge the potential log-to-log variation in relationships that may be present and hence the intercept term was allowed to vary randomly.However, the slope term was treated as fixed.This latter constraint was established during the preliminary model specification phase by initially treating both the intercept and slope terms as random.In cases were convergence could not be achieved, the model specification was changed by defining the intercept as fixed and the slope as random and vice versa.Overall, the results indicated that (1) convergence could not be achieved when both terms were treated as random irrespective of the relationship, (2) convergence could not be achieved when the intercept term was treated as fixed and the slope was treated as random for all relationships but that for wood density, and (3) convergence for all relationships was only achieved when the slope was treated as fixed and the intercept term was treated as random.These statistical results were used to inform the selection of the final model specifications (i.e., Equations (5a) and (5b)) and hence were considered the most applicable for the sample population used.

Goodness-of-fit, Lack-of-fit, and Predictive Ability of Fitted Models
The parameterized models were evaluated based on goodness-of-fit, lack-of-fit, and predictive ability employing the validation data subset.More specifically, the proportion of variability in the dependent variable (kth attribute) explained by the parameterized model as measured by the index-of-fit squared statistic (I 2 (k) ), was used as a goodness-of-fit metric (Equation ( 6)).Parameterized relationships which had mean absolute bias (B a(k) ; Equation ( 7)) or mean relative bias (B r(k) ; Equation ( 8)) that were significantly (p ≤ 0.05) different from zero as determined through the use of the 95% confidence intervals (Equation ( 9)), where considered as demonstrating a consequential lack-of-fit.The linear regression relationship between the observed and predicted values was employed as a secondary lack-of-fit measure in terms of detecting systematic bias (sensu [31]).Based on Reynolds [32] statistical framework in combination with Gribko and Waint [33] SAS macro program extension, prediction and tolerance intervals for absolute and relative error where used to quantify the predictive accuracy of the parameterized relationships (Equations ( 10) and (11), respectively).
where V (k) l and V(k) l are the observed and predicted values, respectively, for the kth attribute within the lth sample log belonging to the validation data subset, n (k) is the number of predicted-observed pairs specific to the kth attribute within the validation data subset, S a,r(k) is the standard deviation of the absolute (S a(k) ) or relative (S r(k) ) biases specific to the kth attribute, respectively, t (n (k) −1,0.975) is the 0.975 quantile of the t-distribution with n (k) − 1 degrees of freedom specific to the kth attribute, n p is the number of future predictions under evaluation (n p = 1; individual log level), and g is a normal distribution tolerance factor specifying the probability (λ; i.e., 0.05) that a specified minimum proportion (i.e., 95%) of the distribution of errors (P ) will be contained within the stated tolerance interval.

Predictive Performance when Deploying Acoustic-Derived Wood Density Estimates
Given the practical reality of not having access to Silviscan-equivalent wood density estimates when acoustic sampling, necessitated the evaluation of the density-weighted acoustic relationships when a surrogate acoustic-derived w d estimate is used.Procedurally, this involved (1) generating a density estimate for each log within the validation data subset utilizing the parameterized wood density prediction model in combination with the acoustic velocity measurement, (2) given (1), deploying the resultant density estimate along with its acoustic velocity measurement to generate attribute predictions for each of the successfully parameterized models, and (3) using the observed and predicted values to calculate (i) mean absolute and relative biases along with associated 95% confidence intervals, and (ii) prediction and tolerance error intervals.The computations used to generate the mean absolute and relative biases and associated confidence intervals, and the prediction and tolerance error intervals, are similar to those described above (i.e., Section 2.6; Equations ( 7)-( 11)).

Attribute-Acoustic Velocity Relationships: Parameter Estimates, Regression Statistics and Tenability of Associated Assumptions
The resultant parameter estimates, regression statistics and validation metrics for assessing compliance with underlying assumptions for the attribute-specific acoustic velocity relationships, as described by the mixed-effects model specifications (Equation ( 5)), are given in Table 7.The successful parameterization of 6 out of the 8 relationships assessed (m e , w t , m a , c o , s a = f γ0,1 , w d v 2 l and w d = f γ0,1 , v 2 l ) in terms of the (1) statistical significance of the parameter estimates and random effect term, (2) proportion of variability in the dependent variables explained (I 2 ), and (3) compliance with the independence and constant variance assumptions underlying the parameterization approach used, as assessed by the degree of spatial correlation among adjacent level Forests 2018, 9, 749 13 of 28 one residuals employing regression analysis and Park's homogeneity of variance test [30], respectively, confirmed the applicability of the chosen mixed-effects regression specification in quantifying the attribute-acoustic velocity relationships for jack pine.6)).d Non-significant (p > 0.05) and significant (p ≤ 0.05) random effects are denoted ns and *, respectively (i.e., testing the null hypothesis that u 0 = 0 versus the alternative hypothesis u 0 = 0 in Equation ( 5)).
e Non-rejection and rejection of the constant variance assumption at the 0.01 probability level are denoted H 0 and H 1 , respectively, as determined from Park's test for homoscedasticity (see text for details).f Non-rejection and rejection of the independence assumption (spatially uncorrelated residuals) at the 0.01 probability level are denoted H 0 and H 1 , respectively, as determined from the residual regression approach ( [30]; see text for details).
Specifically, for these fitted relationships, the intercept and slope parameter estimates and the random effect terms were significantly (p ≤ 0.05) different from zero, and the percentage of variation explained varied from a relatively low minimum value of 8% for m a to a relatively moderate maximum value of 60% for w t .The values for the remaining relationships ordered by magnitude of I 2 , were as follows: 32% for w d , 46% for c o , 47% for m e and 50% for s a .The lack of significant (p ≤ 0.01) regression relationships among adjacent residuals for the log pairs derived from the same sample trees were indicative of the absence of spatial correlation effects (sensu [29]).Furthermore, the constant variance assumption was not rejected given the results of Park's test for homoscedasticity along with the interpretation of the Bayes residual error-predictor bivariate scatterplots.Note, during the initial parameterizations, examination of the circular-shape data point clusters within the residual error-predictor scatterplots, did reveal the presence of 2 potential outliers or influential observations (i.e., these observational pairs were visually apart from the concentrated residual cloud within the scatterplots for attributes m e , w d , m a , w t , c o , and s a ).Thus a review of the field records and laboratory procedures in terms of identifying possible data recording errors, compiling errors or incorrect processing procedures for these suspect data pairs was implemented.Additionally, the attribute values for the other logs sampled from the suspect trees were also examined.Although this inquiry revealed the absence of recording or processing errors, the measurements for the suspect log within the calibration data subset and another one from the same sample tree but within the validation data subset, were substantially different from the measurements for the other logs within their respective data subsets.Based on the acknowledgment that these attributes may be occasionally influenced by uncontrollable and largely undetectable sources of variation such as complex internal knot distributions within logs, it was concluded that these suspect observations should be removed.Consequently, the calibration and validation data subsets were reduced by one observational pair each yielding a final total of 121 and 120 attribute-specific-acoustic velocity observational pairs within calibration and validation data subsets, respectively.The attribute-acoustic velocity observational pairs within the calibration data subset for each of the 8 relationships evaluated are presented in Figure 1A,B.The successfully parameterized models are also superimposed on the observational pairs where applicable.Examination of the subgraphs for the attributes that were not successfully parameterized by the mixed-effects model specification (i.e., d r and d t ; Table 7), reinforced the lack of a functional relationship between each of these attributes and density-weighted acoustic velocity: i.e., random cloud of the d r -w d v 2 l and d t -w d v 2 l observational pairs that were devoid of any obvious linear or nonlinear relationship.Conversely, for the attributes that were successfully parameterized by the mixed-effects model specification (i.e., m e , m a , w d , w t , c o , and s a ; Table 7) which are superimposed on their respective subgraphs confirmed the statistical findings (i.e., relationships were unbiased and in accord with the m e -w and s a -w d v 2 l linear trends).
Forests 2018, 9, x FOR PEER REVIEW 14 of 28 within the calibration data subset and another one from the same sample tree but within the validation data subset, were substantially different from the measurements for the other logs within their respective data subsets.Based on the acknowledgment that these attributes may be occasionally influenced by uncontrollable and largely undetectable sources of variation such as complex internal knot distributions within logs, it was concluded that these suspect observations should be removed.Consequently, the calibration and validation data subsets were reduced by one observational pair each yielding a final total of 121 and 120 attribute-specific-acoustic velocity observational pairs within calibration and validation data subsets, respectively.The attribute-acoustic velocity observational pairs within the calibration data subset for each of the 8 relationships evaluated are presented in Figure 1A,B 7).

Goodness-of-fit and Lack-of-fit Assessment
The goodness-of-fit evaluation consisted of assessing the magnitude of the variability explained by the parameterized models when applied to the validation data subset.Specifically, the I 2 statistic was calculated for the acoustic-based me, wd, ma, wt, co and sa mixed-effects models (Table 8).Results indicated that the models explained a relatively low to moderate level variation (ordered by magnitude): 8% for ma, 30% for wd, 43% for co, 44% for sa, 53% for me and 58% for wt.In agreement with expectation, the values for 5 of 6 attributes were slightly less than those observed for the relationships parameterized using the calibration data subset (cf.Table 7 versus Table 8).The exception being that of the result for me which was slightly greater than that derived from the calibration data subset.Based on the regression results for the linear relationship between the observed and predicted values, there was insufficient evidence to indicate systematic lack-of-fit for any of the 6 parameterized relationships: i.e., intercept and slope values were not significant (p > 0.01) different from zero and l and s a = f (w d v 2 l relationships.The parameterized model is denoted by solid line (Table 7).

Goodness-of-fit and Lack-of-fit Assessment
The goodness-of-fit evaluation consisted of assessing the magnitude of the variability explained by the parameterized models when applied to the validation data subset.Specifically, the I 2 statistic was calculated for the acoustic-based m e , w d , m a , w t , c o and s a mixed-effects models (Table 8).Results indicated that the models explained a relatively low to moderate level variation (ordered by magnitude): 8% for m a , 30% for w d , 43% for c o , 44% for s a , 53% for m e and 58% for w t .In agreement with expectation, the values for 5 of 6 attributes were slightly less than those observed for the relationships parameterized using the calibration data subset (cf.Table 7 versus Table 8).The exception being that of the result for m e which was slightly greater than that derived from the calibration data subset.Based on the regression results for the linear relationship between the observed and predicted values, there was insufficient evidence to indicate systematic lack-of-fit for any of the 6 parameterized relationships: i.e., intercept and slope values were not significant (p > 0.01) different from zero and unity, respectively (sensu [31]).Similarly, based on the mean absolute biases and their associated 95% confidence intervals, there was no evidence of lack-of-fit for any of the 6 relationships (Table 8): mean absolute biases were not significantly (p ≤ 0.05) different from zero as inferred from the 95% confidence intervals.Although approximately equivalent results were obtained when assessing lack-of-fit using mean relative biases, the mean relative bias generated for the m a relationship was significantly (p ≤ 0.05) different from zero (Table 8).
The density-weighed and density-unweighted acoustic velocity-attribute observational pairs within the validation data subset for each of the relationships including those that did not exhibit a significant relationship with acoustic velocity, are graphically presented in Figure 2A,B.The parameterized mixed-effects regression relationships derived using the calibration data subset were also superimposed on the m e , w d , m a , w t , c o and s a subgraphs.It was visually evident that the parameterized relationships were representative of the linear trends between the m e -w l and s a -w d v 2 l observational pairs.The general concordance between the linear trends exhibited by the observational pairs and that established by the regressions, were also suggestive of the absence of systematic bias.For the attributes that were not successfully described by the mixed-effects model specification (i.e., d r and d t ; Table 7), an examination of the subgraphs reconfirmed the statistical result, that is, there was a random scatter of d r -w d v unity, respectively (sensu [31]).Similarly, based on the mean absolute biases and their associated 95% confidence intervals, there was no evidence of lack-of-fit for any of the 6 relationships (Table 8): mean absolute biases were not significantly (p ≤ 0.05) different from zero as inferred from the 95% confidence intervals.Although approximately equivalent results were obtained when assessing lackof-fit using mean relative biases, the mean relative bias generated for the ma relationship was significantly (p ≤ 0.05) different from zero (Table 8).
The density-weighed and density-unweighted acoustic velocity-attribute observational pairs within the validation data subset for each of the relationships including those that did not exhibit a significant relationship with acoustic velocity, are graphically presented in Figure 2A,B.The parameterized mixed-effects regression relationships derived using the calibration data subset were also superimposed on the me, wd, ma, wt, co and sa subgraphs.It was visually evident that the parameterized relationships were representative of the linear trends between the me- w v observational pairs.The general concordance between the linear trends exhibited by the observational pairs and that established by the regressions, were also suggestive of the absence of systematic bias.For the attributes that were not successfully described by the mixed-effects model specification (i.e., dr and dt; Table 7), an examination of the subgraphs reconfirmed the statistical result, that is, there was a random scatter of dr-  7) and superimposed contextual 95% prediction error limits are denoted by the parallel dashed lines.

Figure 2.
Scatterplots illustrating the relationship between each attribute and density-weighted and density-unweighted acoustic velocity (v l ) for the jack pine logs within the validation data subset: l and s a = f (w d v 2 l ) relationships.Note, where applicable, the superimposed parameterized model is denoted by solid line (Table 7) and superimposed contextual 95% prediction error limits are denoted by the parallel dashed lines.a I 2 is the index-of-fit squared (Equation ( 6)).b Testing the (1) null hypothesis that α 0 = 0 versus the alternative hypothesis α 0 = 0 where H 0 and H 1 denote the non-rejection and rejection of the null hypothesis at the 0.05 probability level, and (2) null hypothesis that α 1 = 1 versus the alternative hypothesis α 1 = 1 where H 0 and H 1 denote the non-rejection and rejection of the null hypothesis at the 0.05 probability level, where α 0 and α 1 are derived from the simple linear ordinary least squares regression relationship between observed (y) and predicted (x) values (y = α 0 + α 1 x + ε where ε is an error term [31]).c Mean absolute (Equation ( 7)) and relative (Equation ( 8)) bias and the limits (CL) of the associated 95% confidence interval (Equation ( 9)) where mean values not significantly (p > 0.05) different from zero were indicative of an unbiased relationship (i.e., mean bias ± 95% confidence limits (CL)); note, underlined CL values denote approximations given the presence of significant (p ≤ 0.05) non-normality within the underlying error distribution; absolute error units are attribute-specific: GPa, kg/m 3 , µm, µg/m and m 2 /kg for m e , w d , m a , w t , c o and s a , respectively.d CL for the 95% prediction and tolerance error intervals for absolute and relative errors (Equations ( 10) and ( 11), respectively): mean bias ± 95% CL; specifically, there is a 95% probability that a future error will be within the stated prediction interval and that there is a 95% probability that 95% of all future errors will be within the stated tolerance interval [32,33]; note, underlined CL values denote approximations given the presence of significant (p ≤ 0.05) non-normality within the underlying error distribution.

Predictive Ability
The predictive ability of the successfully parameterized mixed-effects models was quantified via the use of prediction and tolerance error intervals (Equations ( 9) and (10), respectively): mean bias ± 95% confidence limit.The intervals were generated from the absolute and relative bias values generated from the validation data subset.Statistically, these error intervals attempt to quantify the performance of the models when they are actually deployed [32].Specifically, the prediction intervals quantify the boundaries of the expected range of absolute and relative errors when applying the models for a newly sampled jack pine log (e.g., there is 95% probability that the future error will be within the stated prediction interval).The tolerance intervals quantify the boundaries of the range of the expected population of absolute and relative errors which would be generated from the models when repeatedly deploying them to newly sampled jack pine log populations (e.g., there is 95% probability that 95% of all future errors will fall within the stated tolerance interval).

Predictive Performance of Parameterized Models When Deploying Acoustic Generated Wood Density Estimates
Employing the validation data subset, the empirical performance of the density-weighted models in cases where an acoustic-based wood density estimate was used as a surrogate for a Silviscan-based density estimate, was also evaluated.Specifically, based on the functional expressions, m e , m a , w t , c o , s a = f γ0,1 , ŵd v 2 l where ŵd is derived from the w d = f γ0,1 , v 2 l relationship, the following set of computations were carried out: (1) inputting the acoustic velocity measurement value for each log in the validation data subset into the parameterized w d model (Table 7) and generating associated estimates for each attribute; and (2) given (1), treating these resultant attribute estimates as predicted values and the corresponding Silviscan-determined estimates as observed values, absolute and relative biases and associated 95% confidence, and 95% prediction and tolerance intervals, were calculated via the deployment of the computational structure given by Equations ( 7)- (11).

Hierarchical Mixed-Effects Acoustic-Based Attribute Prediction Models for Jack Pine
The mixed-effects linear model for the primary relationship between the dynamic modulus of elasticity and density-weighted acoustic velocity (m e = f γ0,1 , w d v 2 l ) explained 47% of the variation in the m e for the logs within the calibration data subset (Table 7) and 53% of the variation in the m e for the logs within the validation data subset (Table 8).Based on the logs within the validation data subset, the model generated unbiased predictions when assessed using both absolute and relative error measures.The tolerance intervals indicated that 95% of the absolute and relative errors generated from repeatedly applying the parameterized model to a new sample population of jack pine logs would be expected (95% probability level) not to exceed ±2.8 GPa and ±25.0%respectively, when using Silviscan-equivalent w d estimate (Table 8).For the deployment scenario in which an acoustic-based wood density estimate is employed as a surrogate measure for its Silviscan-based counterpart, the model's predictions would be similarly unbiased (Table 9).However, the predictions would be less precise as evident from the wider tolerance intervals for absolute and relative error that were generated when using the acoustic-based density estimate (i.e., 95% tolerance limits for error of ±3.1 GPa and ±27.3%, respectively; Table 9).Newton [15] reported slightly better results in a similarly structured analysis of red pine logs in terms of explanatory power (e.g., explaining 71% of variation in m e for red pine logs versus 50% of variation in m e for jack pine logs) and precision of predictions when deploying either a Silviscan-equivalent estimate of wood density (e.g., relative tolerance error intervals of ±19% for red pine versus ± 25% for jack pine) or acoustic-based estimate (e.g., relative tolerance error intervals of ±19% for red pine versus ± 27% for jack pine).
The results for the mixed-effects linear regression model for the relationship between wood density and acoustic velocity (w d = f γ0,1 , v 2 l ), revealed that the parameterized model generated unbiased predictions (Tables 7 and 8).Overall, however, the model performed moderately in terms of (1) explanatory power given that only 32% and 31% of the w d variation for the logs within the calibration and validation data subsets, respectively, was explained (Tables 7 and 8), and (2) predictive performance given that 95% of absolute and relative errors generated from repeatedly applying the equation to a new sample population of logs would be expected (95% probability level) not to exceed ±51.6 kg/m 3 or ±12.5%, respectively (Table 8).Considerably better results were reported previously for red pine logs [15] in terms of explanatory power (e.g., explaining 80% of variation in w d for red pine logs versus 31% of variation in w d for jack pine logs).However, differences in the precision of the predictions were only marginally better (e.g., relative tolerance error intervals of ±10% for red pine versus ± 13% for jack pine).a Mean absolute (Equation ( 7)) and relative (Equation ( 8)) bias and the limits (CL) of the associated 95% confidence interval (Equation ( 9)) where mean values not significantly (p > 0.05) different from zero were indicative of an unbiased relationship; note, underlined CL values denote approximations given the presence of significant (p ≤ 0.05) non-normality within the underlying error distribution; absolute error units are attribute-specific: GPa, • , µm, µg/m and m 2 /kg for m e , m a , w t , c o and s a , respectively; and * denotes significant bias.b Confidence limits for the 95% prediction and tolerance error intervals for absolute and relative errors (Equations ( 10) and (11), respectively): mean bias ± 95% CL; specifically, there is a 95% probability that a future error will be within the stated prediction interval and that there is a 95% probability that 95% of all future errors will be within the stated tolerance interval (sensu [32]); notes, (1) absolute error units are attribute-specific as stated above, (2) underlined CL values denote approximations given the presence of significant (p ≤ 0.05) non-normality within the underlying error distribution, and (3) stand-level CL intervals denote the 95% prediction limits for the mean error generated when group sampling (e.g., assigning a mean attribute estimate to 30-log sample (pile) derived from a single stand).
Although the results for the relationship where microfibril angle is expressed as a function of density-weighted acoustic velocity (m a = f γ0,1 , w d v 2 l ) generated unbiased predictions (Tables 7 and 8), the parameterized model exhibited a low level of explanatory power and a relatively high degree of imprecision.Specifically, the model only explained 8% of the variation in m a for the logs within both the calibration and validation data subsets, respectively (Tables 7 and 8).The predictive intervals indicated that 95% of absolute and relative errors generated from repeatedly applying the equation to a new sample population of logs would be expected (95% probability level) not to exceed (1) ±6.2 • and ±55.5%, respectively, when using Silviscan-equivalent w d estimates (Table 8), and (2) ±6.2 • and ±55.1%, respectively, when deploying acoustic-based wood density estimates.Although these weak results are superior to those obtained previously for red pine logs in which no viable acoustic-based relationship was found for microfibril angle [15], the results are inferior to the level of variation described by the acoustic-based relationship reported in other studies (e.g., 86% for radiata pine logs [6]).
The relationship for cell wall thickness (w t = f γ0,1 , w d v 2 l ) exhibited the highest levels of explanatory power as evident from the level of variation explained: 60% and 58% for the calibration and validation data subsets, respectively (Tables 7 and 8).Furthermore, the parameterized model exhibited no lack-of-fit issues and produced unbiased predictions at a relatively high level of precision.Specifically, 95% of absolute and relative errors generated from repeatedly applying the equation to a new sample population of jack pine logs would be expected (95% probability level) to not exceed ±0.3 µm and ±11.8%, respectively, when using Silviscan-equivalent or acoustic generated w d estimates (Tables 8 and 9).Newton [15] reported a superior result for red pine logs in terms of explanatory power (e.g., explaining 90% of variation in w t for red pine logs versus 59% of variation in w t for jack pine logs), but much less so in terms of predictive precision.For example, relative tolerance error intervals of ±9% for red pine versus ±12% for jack pine when deploying a Silviscan-equivalent estimate of wood density, and ±12% for red pine versus ±12% for jack pine when using an acoustic-based estimate.The graphical and statistical results for the relationships in which radial and tangential diameters were expressed as functions of density-weighted acoustic velocity, did not support the existence of viable relationships for these attributes.This result is similar to that obtained for red pine logs (i.e., [15]).Combining these regression results with the observed lack of linear associations exhibited in the scatterplots (Figures 1 and 2 for jack pine logs; and Figures 1 and 2 for red pine logs in [15]), suggest that these attributes may be intrinsically unrelated to acoustic velocity for these species.
The regression relationships for fiber coarseness and specific surface area (c o , s a = f γ0,1 , w d v 2 l ) exhibited moderate levels of explanatory power (Tables 7 and 8): (1) 46% and 43% of the variation in fiber coarseness was explained for the logs within the calibration and validation data subsets, respectively; and (2) 50% and 44% of the variation in specific surface area was explained for the logs within the calibration and validation data subsets, respectively.These parameterized models produced unbiased and relative precise predictions when compared with the other attribute models.For the fiber coarseness model, 95% of absolute and relative errors generated from repeatedly applying the parameterized equation to a new sample population of logs would be expected (95% probability level) not to exceed (1) ±49.8 µg/m and ±13.0%, respectively, when using Silviscan-equivalent wood density estimates (Table 8), and (2) ±57.0 µg/m and ±14.9%, respectively, when deploying acoustic-based wood density estimates (Table 9).Similarly, for the model developed for specific surface area, the corresponding precision limits for absolute and relative error were (1) ±34.3 m 2 /kg and ±10.5%, respectively, when using Silviscan-equivalent wood density estimates (Table 8), and (2) ±40.7 m 2 /kg and ±12.5%, respectively, when deploying acoustic-derived wood density estimates (Table 9).
In relation to acoustic-based prediction of fiber coarseness, Newton [15] reported much improved results for red pine logs in terms of explanatory power (e.g., explaining 81% of variation in c o for red pine logs versus 44% of variation in c o for jack pine logs) but again much less so in terms of predictive precision (e.g., relative tolerance error intervals of ±12% for red pine versus ±13% for jack pine when deploying a Silviscan-equivalent estimate of wood density, and ±14% for red pine versus ±15% for jack pine when using an acoustic-based estimate).In comparison of the model developed for predicting specific surface area for red pine logs, the result was considerably superior in terms of explanatory power (e.g., 84% of the variation in s a explained for red pine logs versus 47% of the variation in s a explained for jack pine logs).However, similar to the interspecies comparison of the coarseness model, the difference in predictive ability was minimal between the species: (1) relative tolerance error intervals of ±8% for red pine versus ±11% for jack pine when deploying a Silviscan-equivalent estimate of wood density, and (2) ±10% for red pine versus ±13% for jack pine when using an acoustic-based estimate.

Potential Utility of the Expanded Acoustic-Based Inferential Framework for Jack Pine
Estimation of the modulus of elasticity is of primary importance in log-segregation operations where the objective is to forecast the potential type and grade of extracted solid wood products (e.g., dimension lumber).However, other wood attributes (e.g., wood density, microfibril angle, tracheid dimensions) and derived composite metrics (e.g., fiber coarseness and specific surface area), are also prime determinates of end-product potential (sensu Table 1).Consequently, the expanded acoustic inferential framework presented in this study along with its empirical validation for six of the eight jack pine attributes examined (dynamic modulus of elasticity, wood density, microfibril angle, cell wall thickness, fiber coarseness and specific surface area), yields a more comprehensive system for non-destructive forecasting of end-product potential than systems based solely on wood stiffness (modulus of elasticity).Operationally, however, the in-forest deployment of the parameterized relationships will be largely dependent on the objectives and precision requirements of the end-user.For example, the relative error intervals were quite large for m e and m a as exemplified by the resultant tolerance errors of ±25% and ±56%, respectively, when using a Silviscan-equivalent wood density estimate (Table 8), and ±27% and ±55%, respectively, when using an acoustic-based wood density estimate (Table 9).Thus, the ability to stratify individual jack pine logs into narrow-grade class-based acoustic estimates of m e and m a would be difficult if not impossible.However, in situations for which the objective is to segregate groups of logs into end-product categories according to their average attribute values, then the magnitude of estimation error associated with the acoustic-based population mean estimates may be more operationally acceptable.For example, the mean absolute and relative error to be expected when collectively sampling 30 new jack pine logs deploying an acoustic-based wood density estimate would be, respectively, 0.6 GPa and 5% for m e , and 1 • and 10% for m a (Table 9).
Comparably, the relative tolerance error intervals for w d , w t , c o and s a were considerably less (i.e., errors in the order of ±10 to ±15% (Table 8)).Thus, depending on the specified width of the w d , w t , c o and s a intervals the end-user employs to designate logs into specific end-product categories or a binary grade class, differentiating logs according to these attribute predictions may be feasible.Similarly, at the stand-level where the objective is to classify a population of logs in a specific end-product category based on their mean w d , w t , c o or s a values, the precision of the estimates would be improved considerably.For example, the mean absolute and relative error to be expected when collectively sampling 30 new jack pine logs deploying an acoustic-based wood density estimate would be respectively, 0.1 µm and 2% for w t , 10.4 µg/m and 3% for c o , and 7.5 m 2 /kg and 2% for s a (Table 9).
More generally however, based on the explanatory power and precision of acoustic-based attribute prediction models developed to date, increasing the precision of the point-estimates derived from acoustic-based models will be required if they are to be used to stratify individual logs into grade categories given the relatively narrow range of static modulus of elasticity values that currently delineate product grade classes.For example, the mean difference between the 14 machine stress-rated (MSR) lumber grade categories established for Canadian softwood species [34] is approximately 0.7 GPa and hence the explicit stratification of individual logs into such narrow-width grade classes would not be possible given the empirical results to date.Even if it is assumed that the dynamic estimate is equivalent to its static counterpart which in turn is reflective of the actual end-product-based value, the precision of most acoustic-based estimates would be inadequate for stratifying logs into MSR-based classes (e.g., ±3 GPa expected error for individual jack pine logs; Table 9).Thus reducing the number of grade categories from 14 to a smaller number of discrete quality classes based on the expected prediction error range(s) through clustering, and (or) using multiple individual log estimates to generate a mean site, landing or log pile value, may represent a viable alternative when using the acoustic approach in log-segregation operations for jack pine.
Further research into analytical advancements that increase the amount of variability explained and decrease the standard error of estimate, possibly by including covariates within the acoustic-attribute model specification as suggested by Bérubé-Deschênes et al. [35] and Butler et al. [36], may be warranted.Acoustic velocity is known to be affected by a multitude of internal and external factors which includes knot distributions, embedded voids, environmental conditions during acoustic sampling (temperature and moisture conditions), tree size, local competition, and overall site conditions (e.g., [35,37,38]).Provision of operational solutions for minimizing these effects could improve the predictive accuracy of acoustic-based attribute estimation.Recent results examining the effect of xylem temperature and moisture on acoustic velocity within standing semi-mature jack pine trees during the vegetative growing season indicated that acoustic velocity declined in linear fashion with increasing temperature; however, moisture had no appreciable influence [39].However, the temperature effect was not of consequential significance except when temperatures approached their seasonal extremities.Assuming similar inferences apply to jack pine logs, suggest that acoustic log samplers could treat such variation as a source of random error of minimal importance providing xylem temperatures did not approach their seasonal minimums (<5 • C) or maximums (>30 • C).In cases were temperature effects were of concern, the standardization equation for adjusting acoustic velocities to a reference xylem temperature of 20 • C could be deployed [39].

Similarities and Differences between Tree and Log Acoustic-Based Attribute Relationships
The time-of-flight acoustic approach to estimating the internal attributes within standing trees employs the dilatational stress wave velocity (v d ) which transverses the breast-height (1.3 m) region of the main stem [40].This acoustic velocity measurement when weighted by wood density is also related to the modulus of elasticity as described by a functional specification, similar to that given by Equation (1) (i.e., replacing v l by v d ).This primary relationship for standing trees can be also potentially expanded to include a similarly structured suite of secondary relationships as that presented in this study for logs (e.g., [16]).
Conceptually, the resonance-based acoustic velocity approach is considered to yield a more representative estimate of the fiber attributes within the xylem tissue than those estimates derived from the time-of-flight approach.This is principally because the mean velocity estimate is derived from many acoustic pulses resonating longitudinally throughout the log rather than the velocity of a single wave front moving between two points within the main stem of a standing tree [41].The empirical validity of this inference can be partially tested by assessing the differences between the tree and log acoustic-based attribute prediction models parameterized using the same sample tree population.Contrasting the results derived from a companion analysis which assessed the time-of-flight acoustic relationships deploying the same suite of attributes and sample tree population which were measured at the same time [42], revealed that the time-of-flight-based relationships were mostly superior in terms of explanatory power.Specifically, the percentage of attribute variation explained by the dilatational (tree) and longitudinal (log) acoustic relationships were respectively: 71% versus 50% for m e ; 30% versus 31% for w d ; 19% versus 8% for m a ; 66% versus 59% for w t ; 42% versus 44% for c o ; and 61% versus 47% for s a .However, conversely, in terms of predictive precision, the comparison between the log and tree models revealed minimal differences.The error arising from using a newly sampled acoustic velocity measurement used to estimate m e , w d , m a , w t , c o and s a for standing trees and derived logs would be expected at the 95% confidence level to fall within the respective relative error limits of ±25% versus ±21% for m e , ±13% versus ±14% for w d , ±56% versus ±47% for m a , ±12% versus ±11% for w t , ±13% versus ±12% for c o , and ±10% versus ±10% for s a .These similarities in terms of the predictive performance between the tree and log-based acoustic approaches for estimating attributes within the same species (jack pine) while attaining contrasting results between species (jack pine versus red pine for a given time-of-flight-or resonance-based relationship), suggest that acoustic-based attribute relationships may be intrinsic to a given species irrespective of either of these wave types.

Advancing Acoustic-Based Attribute Estimation
The resonance-based acoustic approach to estimating internal fiber attributes within logs and associated demonstrations of its utility in segregation, allocation and merchandizing operations have been presented to the forest management community through various case studies (e.g., [41,43,44] and comprehensive literature reviews (e.g., [11,17]).These contributions have mostly focused on the merits of estimating the modulus of elasticity and using it, or its surrogates, to forecast end-product potential.These include (1) treating wood density as a constant and inferring wood stiffness indirectly from acoustic velocity (e.g., implicitly based on the underlying theoretical relationship expressed by Equation (1)), (2) estimating wood stiffness using parameterized regression models in which the dynamic modulus of elasticity is expressed as a function of acoustic velocity and wood density (e.g., explicitly based on the conceptual formulation (Equation ( 1)), and (3) directly forecasting wood stiffness using relationships that explicitly relate the static modulus of elasticity of a given end-product (e.g., measured using static bending tests of sawn boards) to acoustic velocity.This last approach while the most informative is also the least developed.Consequentially, the full potential of in-forest acoustic-based forecasting of end-product potential has not been fully realized.Determination and quantification of the explicit relationship between log-based acoustic velocity measures and actual performance-based metrics within the derived manufactured end-product are also required.Preliminary results from several studies in which the static modulus of elasticity within dimensional lumber products have been explicitly related to log-based acoustic velocity measures have been positive.These include statistically viable relationships reported for balsam fir (Abies balsamea (L.) Mill.) and white spruce (Picea glauca (Moench) Voss.) [45], radiata pine [46], Douglas fir ((Pseudotsuga menziesii (Mirb.)Franco) [47], and loblolly pine (Pinus taeda L.) [48].
Similarly, previous research has shown that several important performance metrics associated with pulp and paper end-products are also empirically related to acoustic velocity measurements taken on logs.Clark et al. [48] found that fiber length, pulp strength and various handsheet properties derived from radiate pine logs varied systematically with acoustic velocity.Bradley et al. [49] demonstrated that the variation in the pulp quality (freeness) for radiata pine logs could be considerably reduced by stratifying the logs according to their acoustic velocity values before processing.Although beyond the scope of this study, establishment of a broader range of more explicit relationships between log-based acoustic velocity measures and attribute-derived performance metrics within actual manufactured end-products, would be an area worthy of further investigation.As these past studies have shown, in-forest acoustic velocity stratification of log populations not only improves segregation efficiency at the time of harvest but also assists in reducing the variation in end-product quality during the manufacturing and merchandizing stages.
In-forest acoustic log grading has been increasingly used to compliment visual-based approaches in the pursuit of deriving economically efficient sorting networks.This remains an iterative and ongoing process in wood allocation decision-making where the objective is to increase operational efficiency within the upstream portion of the forest products supply chain (sensu [50]).Mechanized acoustic sampling has also advanced to the point where onboard resonance tools have been installed directly on harvesting machines where logs are immediately sorted upon bucking according to their end-product potential (e.g., [51]).Although these innovations have advanced the in-forest non-destructive approach to log grading and sorting, the acoustic approach has been largely limited to providing a single measure of internal wood quality, that being wood stiffness.Thus, the confirmatory empirical results presented in this study for jack pine and previous for red pine logs [15] are collectively supportive of an expanded acoustic-based inferential framework in which estimates for a multitude of end-product-based attribute determinates can be attained: e.g., wood density, microfibril angle, tracheid wall thickness, fiber coarseness, and specific surface area, in addition to the dynamic modulus of elasticity (wood stiffness).Although further research is required to fully realize the benefits of this empirical-based framework, the ability to non-destructively forecast a wide array of commercially relevant attributes could have consequential potential utility in advancing the acoustic approach in forest operations.

Conclusions
Given the wide spread occurrence of jack pine across the boreal landscape combined with the vast array of potential end-products it can produce, inclusive of solid wood products (e.g., dimensional lumber) and associated mill-work derivatives (window frames, doors, shelving, moldings, and paneling, and composite lumber products such as glulam-based beams, headers and heavy trusses and finger-jointed joists and rafters), and pulp-derived products such as paperboards, newsprint, facial tissues and specialized coated papers [52], the species has become the dominant feedstock species for numerous industrial conversion facilities.However, this diversity of end-products complicates in-forest segregation, allocation, and merchandizing decision-making.Hence the provision of enhanced operational intelligence arising from in-forest forecasts of end-product potential of harvested logs through non-destructive acoustic-based methods, may yield increased efficiencies within the upper portion of the jack pine forest products supply chain.Consequentially, the development and evaluation of a suite of acoustic-based models for predicting the principal attributes governing end-product potential for jack pine as presented in this study, represents an incremental contribution towards more informed decision-making.Specifically, deploying a mixed-effects linear modeling approach combined with cross-validation techniques, viable forecasting models for predicting the dynamic modulus of elasticity, wood density, microfibril angle, cell wall thickness, fiber coarseness and specific surface area were developed.Although these positive results confer additional empirical support for the proposed acoustic-based inferential framework, further research in the areas of accounting for environmentally induced wave variation, specifying end-product-based design thresholds, and explicitly establishing linkages between log-based attribute estimates and those within recoverable end-products, would be beneficial.Collectively, the results presented here for jack pine not only provides the prerequisite parameterized relationships for improving in-forest segregation and allocation decision-making but also contributes to solidifying the empirical foundation of the expanded acoustic-based inferential framework.
Funding: Canadian Wood Fiber Centre, Canadian Forest Service, Natural Resources Canada.

Figure 1 .
Figure1.Scatterplots illustrating the relationship between each attribute and density-weighted and density-unweighted acoustic velocity (v l ) for the jack pine logs within the calibration data subset:(A) m e = f (w d v 2 l ), w d = f (v 2 l ), m a = f (w d v 2 l ), and w t = f (w d v 2 l ) relationships; and (B) d r = f (w d v 2 l ), d t = f (w d v 2l ), and c o = f (w d v 2 l and s a = f (w d v 2 l relationships.The parameterized model is denoted by solid line (Table7).

Table 1 .
Product-based performance measures and their relationship with fiber attributes for boreal conifers.

Table 2 .
Linear associations between the dynamic modulus of elasticity and other Silviscan-derived wood quality attributes for 3 boreal conifers and the resultant acoustic-based relationship inferred.

Table 3 .
Descriptive statistical summary of the mensuration characteristics of the 61 sample trees by site (n = 31 and 30 for Sewell and Tyrol, respectively).

Table 4 .
Descriptive statistical summary of the dimensions of the 243 sample logs and associated acoustic velocity measurements by log position (n = 243 of which 61, 61, 61 and 60 were first, second, third and fourth positioned logs, respectively).
b Coefficient of variation.

Table 5 .
Descriptive statistical summary of the fiber attributes of by log position (n = 243 of which 61, 61, 61 and 60 were first, second, third and fourth positioned logs, respectively).
a Proceeding upwards from the stump to the stem apex, 1st, 2nd, 3rd and 4th denotes the ordinal position of the first, second, third and fourth log extracted from the main stem of the 61 jack pine sample trees.b Coefficient of variation.

Table 6 .
Descriptive statistical summary of the calibration (n = 122 logs) and validation (n = 121 logs) data subsets.
a Coefficient of variation.

Table 7 .
Attribute-specific regression results for the mixed-effects model specification (Equation (5)): parameter estimates, regression statistics and compliance indicators.a random effect intercept parameter estimate specific to the kth attribute and γ1 is fixed effect slope parameter estimate specific to the kth attribute.Note, parameter estimates not significantly (p > 0.05) different from zero are superscripted by ns.b Degrees of freedom for regression (n reg ) and residual error (n res ).
c I 2 is the index-of-fit squared (Equation ( . The successfully parameterized models are also superimposed on the observational pairs where applicable.Examination of the subgraphs for the t ; µm) (A) 2 l and d t -w d v 2 l observational pairs devoid of discernable linear or nonlinear trends.

Table 8 .
Evaluation metrics for the parameterized models: goodness-of-fit, lack-of-fit statistics and predictive ability metrics.

Table 9 .
Predictive performance of parameterized models when deploying acoustic-based wood density estimates.