Equations of bark thickness and volume profiles at different heights with easy-measurement variables

The objective of this work was to develop equations of thickness profile and bark volume at different heights with easy-measurement variables, taking as a study case Nothofagus pumilio forests, growing in different site qualities and growth phases in Southern Patagonia. Data was collected from 717 harvested trees. Three models were fitted using multiple, non-lineal regression and generalized linear model, by stepwise methodology, iteratively reweighted least squares method for maximum likelihood estimation and Marquardt algorithm. The dependent variables were diameter at 1.30 m height (DBH), relative height (RH) and growth phase (GP). The statistic evaluation was made through the adjusted determinant coefficient (r2-adj), standard error of the estimation (SEE), mean absolute error and residual analysis. All models presented good fitness with a significant correlation with the growth phase. A decrease in the thickness was observed when the relative height increase. Moreover, a bark coefficient was made to calculate volume with and without bark of individual trees, where significant differences according to site quality of the stands and DBH class of the trees were observed. It can be concluded that the prediction of bark thickness and bark coefficient is possible using DBH, height, site quality and growth phase, common and easymeasurement variables used in forest inventories.


Data measurement
Data was collected in San Justo Ranch (54° 06' LS, 68° 37' LO) in Tierra del Fuego, Argentina.Five stands were selected along the site quality range: forests growing in SQ I (site index at a base age of 60 years -SI 60 = 19.8 to 23.2 m) present total volumes with bark (TVwb) of more than 1,100 m 3 .ha - and more than 27.5 m height in forests in maturing phase; in SQ II (SI 60 = 16.5 to 23.2 m) present up to 900 m 3 .ha - TVwb and height of 24.0 to 27.5 m in maturing phase; in SQ III (SI 60 = 13.1 to 16.5 m) TVwb reaches 700 m 3 .ha - and heights between 20.5 and 24.0 m in maturing phase; in SQ IV (SI 60 = 9.8 to 13.1 m) TVwb is up to 550 m 3 .ha - and trees with total heights between 17.0 and 20.5 m in maturing phase; and in SQ V (SI 60 = < 9.8 m) TVwb is less than 400 m 3 .ha - and trees present a total height of less than 17.0 m in maturing phase (Martínez Pastur et al., 2009).Since SI 60 is a variable of complex measurement, because it requires age determination, SQ was used instead.SQ can be estimated with dominant height on maturing phase making it an easy-measurement variable (Martínez Pastur 2005).
A total of 717 trees were cut, and the following variables were measured: normal diameter at 1.3 m -DBH (cm), paired data of diameters -Di (cm) and height -Hi (m) along the stem, total height of the tree -TH (m), bark thickness for each Di-Hi -BTi (cm).Also, data of Hi was expressed as a relative height -RH, considering as a ratio of each Hi and TH of each tree.Measures were taken with metric tapes over cut faces of each log, with a precision of ± 1 mm.Heights where BTi were measured and the number of measurements in each tree (1 to 3) were randomly determined to avoid autocorrelation.Moreover,

Introduction
Bark is defined as the layer between the cambium and the outer limit of the stems (Esau, 1969).Different from heartwood, bark changes its properties continuously with age, where thickness varies due to environmental, genetic and biological factors (Esau, 1969;Trockenbrodt, 1991).The bark proportion can reach up to 10% to 20% of the volume of harvested stems.Bark amount is a relevant variable in the planning of primary transformation of wood, and could be a problem for some industries.However, during the last years, it has become more important for biomass energy production (Adler, 2007).
Volume is the most known way to express the amount of harvested timber, where the amount of bark is usually estimated by differences between volumes obtained from diameters measured with and without its inclusion.It has been proposed that the amount of bark and its thickness across the stem might be related with easy-measurement variables in standing trees (Stayton and Hoffman, 1970), with the existence of some models that estimate bark thickness (Hamilton and Chikumbo, 1997;Laasasenaho et al., 2005).However, these models include age determination during forest inventories, making difficult its use.
Moreover, thickness is affected by site quality, diameter class and age.For example, thickness is proportionally lower in younger trees than in older ones (Gea et al., 2004).A model of bark profile including easy-measurement variables would allow more adjusted estimations of harvested volume without bark.The objectives of this work were to analyze the behaviour of bark thickness, to develop equations of thickness profile as a function of easy-measurement variables, and to generate a tool to predict volume without bark, taking as study case the species N. pumilio, growing in different site qualities and growth phases in southern Patagonian forests.The tested hypothesis were that bark thickness changes with tree height, since diameters of the stem at higher heights present lower ages and in consequence less bark thickness.Also, trees with greater diameter at 1.30 m height (DBH) would have thicker barks than smaller trees, while trees in worst quality sites with lower increment in DBH would have more bark thickness than trees growing in better quality sites.

Material and methods
Nothofagus pumilio (Poepp.et Endl.)Krasser forests are considered the most important timber resource in Equations of bark thickness and volume profiles at different heights ing to Schmidt and Urzúa (1982) since it is a useful estimation of the age influence, where 4 categories were defined: initial growth phase (IGP), final growth phase (FGP), maturing growth phase (MGP) and decaying or senescence growth phase (DGP).Finally, crown class of each tree in 4 categories were also defined according to Hildebrand-Vogel et al. (1990) and Martínez Pastur et al. (1994;1997): dominant (D), codominant (C), intermediate (I) and suppressed (S).Both categorizations allow including these variables into the studied models as dependent variables expressed as entire numbers or as factors.

Model construction
The influence of the measured variables (TH: total height of the tree, DBH: diameter at 1.30 m height, Di: diameter at a Hi, RH: relative height (Hi/TH), SQ: site quality, GP: growth phase, C: crown class, TVwb: total volume with bark) over bark thickness at different heights (BTi) was studied by graphic analyses of mean dispersion, and Pearson's correlation matrix with p < 0.05.Multiple and non-lineal regression techniques were used to model bark thickness through the stem (BTi) (Hamilton and Chikumbo, 1997;Laasasenaho et al., 2005), as well as generalized linear model (Mglm) (Hastie and Tibshirani, 1986).To select the variables for the multiple lineal regression model (MM), a forward stepwise method was performed (F = 50).Based on the variables selected by Pearson's correlation matrix and stepwise methodology, a number of non lineal models (Mnl) was obtained using these mathematical forms: where: a, b, c and d are parameters of the model, and X, Y and Z are variables selected using a forward stepwise method (F = 50).The model with the best fitting was chosen.The estimation of the parameters in the non lineal regression model was made with the algorithm proposed by Marquardt (1963).The selected variables were also used in the Mglm, using maximum likelihood estimation for inference of parameters (Hastie and Tibshirani, 1986).

Bark volume estimation
Using the paired data of diameters -Di (cm), and bark thickness -BTi (cm) along the stem, the rela-tionship between the area of the stem section without and with bark was calculated for each tree.The average per tree of all the measurements of each area of the section of the stem with and without bark is BR, which represents, in an individual tree, the relationship between volume without and with bark.Using the variables measured in the standing tree (TH, DBH, RH, SQ, GP and C) multiple and non-lineal regression as well as generalized linear models were generated to predict BR.The behaviour of this relationship was described by analysis of the variance using DBH classes, site and growth phases as factors.To select the variables for Bark volume estimation, the same methodology presented in Model construction was used.

Statistics and model validation
Statistic evaluation of the models was made by determination coefficient (r 2 -adj.), standard error of the estimation (SEE), absolute mean error (AME) and residual analysis.With the adjusted models an autovalidation was carried out, an analysis of mean errors (%ē) and absolute values (%|ē|) along diameter frequencies, growth phases and site qualities of the stands were conducted and expressed as: where: n is the data number, e i is the difference between observed and predicted values.All statistical analyses were conducted with R (R Development Core Team, 2011).
Assumptions of residual normality and homoscedasticity were analyzed by means of graphic analysis.Durbin-Watson statistic was calculated for serial correlation analysis.

Model of bark thickness through the stem
Mean, minimum and maximum values and standard deviation of the observed variables in the sample were covered satisfactorily throughout their distributions (Table 1).These values presented a slight tendency towards the inferior classes, although variations were small when standard deviation was analyzed.
In the Pearson's correlation matrix the association between the studied variables was determined (Table 2).Most variables presented a significant correlation with BTi, GP, DBH, TH and TVwb with a direct relationship, which means that bark thickness increases with tree size.In the case of C, SQ and RH the relationship was negative, which means that in lower site quality stands the bark thickness is lower, and bark thickness diminished when RH and C increased (e.g., dominant trees present higher bark thickness).
The variables resulting from the stepwise methodology were DBH (T = 23.24p < 0.0001), GP (T = 8.32 p < 0.0001) and RH (T = 9.62 p < 0.0001), and the final model was obtained with an adjusted r 2 of 0.54 and a standard error of 0.25, with a decrease in the fitness comparing with the model using TH, DBH, RH, SQ, GP, C and TVwb (r 2 = 0.58; SEE = 0.2425).
The multiple lineal regression model (MM) presents the following mathematic expression: where: BTi is bark thickness at each height (cm), DBH is diameter at 1.30 m height (cm), RH is relative height, and GP is growth phase (1 to 4).The model presents a mean absolute error of 0.20 cm DBH and growth phase increased with BTi and decreased with RH, when Pearson's correlation matrix was analyzed.Based on these results, a non-lineal model was proposed which included DBH, GP and RH as independent variables.
where: BTi is bark thickness at each height (cm), DBH is diameter at 1.30 m height (cm), RH is relative height, GP is growth phase (1 to 4).
The non-linear model presents similar fitness values as the lineal proposal (r 2 adj.= 0.52, SEE = 0.26 and AME = 0.20).BTi = 0.420301 -0.161387 • GP(1) -0.0700347 • GP(2) + + 0.00888027 • GP(3) -0.210024 • RH + 0.0148199 • DBH where: BTi is bark thickness at each height (cm), DBH is diameter at 1.30 m height (cm), RH is relative height, and GP is growth phase (1 to 4).GP Table 1.Statistic characterization of the sample.TH: total height of the tree (n = 717), DBH: diameter at 1.30 m height (n = 717), BTi: bark thickness at each height (n = 1325), Di: diameter at each sampled height (n = 1325), RH: relative height (Hi/TH) (n = 1325), SQ: site quality (n = 717), GP: growth phase (n = 717), C: social class (n = 717), TVwb: total volume with bark (n = 717), SD: standard deviation    Residual analyses were done using the same data base (autovalidation analysis), in order to analyse the behaviour of the models along site quality, DBH frequencies and growth phase (Fig. 1).Similar mean errors (%ē) and absolute values (%|ē|) distributions were observed when site quality, DBH and growth phase classes were considered in the models, except for decaying phase in the MM where %ē had larger errors in linear compared to non-linear model.The Mglm presented the best fittings in %ē along all growth phases.Graphical techniques of residual analysis (Fig. 2) showed that error was normally distributed, the error variance was constant and errors were independent.

Bark volume model
All variables presented a significant negative correlation with BR, except GP, TH, DBH and TVwb which present a positive correlation, meaning that BR increases with tree size (Table 3).
The multiple lineal regression model (MM), and the non-lineal model (Mnl) presents the following mathematic expression with SQ, TH and GP as independent variables:    where: BR is relationship between volume without and with bark, TH is total height of the tree (m), SQ is site quality, GP is growth phase (1 to 4).
where: BR is relationship between volume without and with bark, TH is total height of the tree (m), SQ is site quality, GP is developmental phase (1 to 4).GP Multiple regression (r 2 adj.= 0.19, SEE = 0.018 a n d A M E = 0 .0 1 4 ) , n o n l i n e a l m o d e l s ( r 2 adj.= 0.18, SEE = 0.020 and AME = 0.015) and M g l m ( r 2 a d j .= 0 . 1 7 , S E E = 0 .0 2 0 a n d AME = 0.015) presented r² values never higher than 0.2.Distributions of the observations according to the predicted and observed values showed a non explicative pattern of the dependant variable.The D-W value is near 2 for all the models analyzed (MM D-W = 2.03245, Mnl D-W = 1.93688,Mglm D-W = 1.99173).However, differences between means were found by Analysis of Variance for BR, according to diameter class (p < 0.0001; F = 8.17) and site quality (p < 0.0001; F = 25.22) (Fig. 3).The observed tendency showed an increase in BR while DBH of the tree is higher.In poor site qualities (SQ = V), BR presented the lowest values, while in SQ = I and II these values were significantly higher.

Discussion
Bark volume depends on its thickness and stem diameter, being affected by tree age and stand conditions (environmental factors), where bark thickness increases (Laasasenaho et al., 2005;Sonmeza et al., 2007) as DBH increases (with the age of the individuals).Dimitrov (1976) found that the best models for bark volume estimation for Picea included DBH, total height, volume and site quality of the stands.On the other hand, Johnson (1966) observed in Pseudotsuga that variables determining bark volume are different in the higher part of the stem than in the lower part, and suggested a polynomial for determining the upper stem bark.In the present study we found that volume is influenced by site class of the stands and diameter class of the trees, being the inclusion of these variables fundamental for a good TVwb estimation.Stayton and Hoffman (1970) found for Acer saccharum that diameter relationship with and without bark (k) decreases with tree height.Moreover, this relationship is influenced by tree age.These authors reported high variations to establish a model that explains bark profile, probably because A. saccharum presents different types of bark according to tree age (Sajdak, 1968).The model generated in this work proposes the use of growth phase as independent variable, where each category differentially influences bark.This categorization increases the accuracy of BTi estimation, being a variable of easy estimation during forest inventories.In Picea orientalis it was observed that bark thickness at DBH is not only correlated with tree age and diameter class but also with stand location in sunny or shadow hillsides (Sonmeza et al., 2007).
Autovalidation of the models might lead to wrong conclusions since the assumptions of the statistic models are not confirmed with an independent sample.However, it gives a first impression of how effective can a model be through the different analyzed gradients.Nonetheless, we suggest, in future applications of these models by potential users, the realization of independent validations to confirm the accuracy of the estimations.
In the present work models presented a similar fitness of the determinant coefficient and in the SEE and AME through SQ, GP and DBH, having the Mglm a better distribution of errors through GP.Even though Mnl and Mglm presented errors in the order of 0.26 cm of BT, it is feasible the use of this kind of models in standing trees.The search of a model to predict BR presented encouraging results using a value of BR for different DBH or site classes.

Conclusions
It can be concluded that it was possible to predict BTi in N. pumilio using DBH, RH and growth phase as independent variables, which are easy to measure during forest inventories.All models developed were representative of the mean real values obtained in the sampling, which can adapt to the gradients of the analyzed variables.Also, models present a biological significance maintaining a separation for different growth phases in the used range.

Figure 2 .
Figure 2. Distribution of the residuals according to observed and predicted values.MM: multiple lineal regression model; Mnl: non lineal model; Mglm: generalized linear model.

Table 2 .
Correlation matrix to analyse