Compatible System for Predicting Total and Merchantable Stem Volume over and under Bark , Branch Volume and Whole-Tree Volume of Pine Species

Accurate quantification of branch volume in trees is important for sustainable forest management, especially as these fractions are increasingly used for bioenergy, and for precise forest CO2 quantification. Whereas a large focus has been placed on the compatible estimation of tree taper and bole volume with and without bark, little effort has been made to develop models that allow a simultaneous prediction of these variables together with tree branch volume. In this study, 595 Pinus cooperi trees and 700 Pinus durangensis trees were sampled in pine-oak forests in the Sierra Madre Occidental, Mexico. A compatible system for predicting two segmented taper functions, over and under bark; the corresponding merchantable volumes; coarse branch volume and whole-tree volume was fitted using a modified continuous autoregressive structure to account for autocorrelation. The proposed compatible equations explained more than 97% of the observed variability in diameter over and under bark, volume over and under bark, and total tree volume and more than 64% of the observed variability in branch volume in both species. The method described can theoretically be replicated for any tree species, thus providing a better understanding of the patterns of volume distribution by components, potentially improving carbon accounting system and forest bioenergy planning.


Introduction
In contrast to a large body of literature focusing on modelling the merchantable stem volume (e.g., [1,2]), the modelling of the volume of branches has received much less attention [2,3].Whole-tree volume models that include the volume of fractions such as tree branches are required in order to provide information about the sustainable use of whole trees [2][3][4][5][6].This could lead to an improvement of our understanding of synergies and tradeoffs between different types of forest ecosystem services, such as carbon storage, versus timber production [7].Such models would help in decision-making associated with the sustainable extraction of more biologically renewable fuels from forests (e.g., [2,8]).This is particularly relevant considering the potential impacts of the removal of such fractions, increasingly used as a bioenergy feedstock, on nutrient and carbon cycling (e.g., [9][10][11][12][13]). Accurate information about bark and branch volume is also very important for planning silvicultural practices, transportation and storage logistics in relation to use of these fractions for bioenergy production (e.g., [14][15][16][17]).
Most of the research on tree volume has focused on the modelling of stem total volume and volume profiles, commonly through the use of taper equations (e.g., [18][19][20][21][22][23]). Such equations predict the upper stem diameter at any point on the tree bole, with the aim of determining possible end uses for a single log.Therefore, they are important tools for establishing the most appropriate type of silvicultural management for the desired tree dimensions and associated end uses (e.g., [18][19][20][21][22][23]). Different stem profile modelling approaches have been developed during the last decades for this goal.The majority of the models generate a solid of revolution from a single continuous function with an exponent that may change from ground to top (e.g., [24][25][26][27][28][29]). A number of inflection points are defined in this approach to describe the species-specific point at which changes occur and to establish different functions for each segment (e.g., [30][31][32][33]).Other approaches are based on considering tree taper within a cause-effect context and then, assume a specific stem development theory to link the growth processes and the resulting realization of tree spatial structures (e.g., [34][35][36]).
Compatible taper and stem volume equations are readily obtained by including a constraint that ensures that the total stem volumes obtained by integrating the taper function are equal to those obtained with a volume equation (e.g., [33,37,38]).This is desirable because of the simplicity and usefulness of total volume equations when classification of the products is not required [22,39,40] and because the parameters of such algebraically consistent systems are often meaningful and better reflect the underlying biological and physical relationships exhibited by the system [41].Consequently, several studies have focused on the compatibility between the diameter reduction described by the taper function and the total stem volume (e.g., [33,39,42,43]) and some works have also extended compatibility to bark (e.g., [22,32]).Compatible adjustment is also considered important for biomass fractions prediction (e.g., [13,44,45].In contrast to the relatively large body of literature on compatible stem and taper prediction (e.g., [42,43]), little effort has been made to develop compatible models for predicting the volume of important tree fractions such as branches (e.g., [4][5][6][7]).Clearly, the next phase in tree volume modelling should be to describe the distribution of the volume of the different fractions of the whole tree, not only of tree boles [3].This approach should aim for compatibility between the estimates of tree biomass and volume of the different fractions, potentially improving carbon accounting systems [12,[46][47][48][49].
Although it has long been recognized that crown size is one of the main determinants of tree bole shape [50], there are few quantitative theories relating crown characteristics and stem form development [51].Various approaches have been used to model the effect of crown architecture on stem taper variation such as the inclusion of crown characteristics as independent variables in stem taper curves [25,52,53]; using different models above and below a relative crown height [3,54] or modelling the effect of crown size on the centroid of volume for different trees [2].
In Mexico, previous studies have compared different taper equations for several species (e.g., [55][56][57][58]).Corral-Rivas et al. [56] and Quiñonez-Barraza et al. [58] compared several taper equations for the simultaneous prediction of diameter and total volume in several pine species in Mexico.Both studies found that the model proposed by Fang et al. [33] yielded the best predictions and that species-specific parametrization of the model was more effective than genus-level parametrization for improving the estimates obtained.However, none of the above-mentioned studies in Mexico attempted neither to predict taper and over and under bark stem volume simultaneously with branch volume and total tree volume nor to consider the effect of crown architecture on stem taper variation.
The objectives of the present study were to simultaneously fit a system of equations for predicting taper and merchantable over and under bark volume, coarse branch volume and total tree volume, for two pine species growing in the state of Durango, Mexico: Pinus cooperi and Pinus durangensis and to analyse the effect of crown size on the values of the inflexion points and the stem factors of the taper equations.

Data
The study was conducted in the part of the Sierra Madre Occidental that lies within the State of Durango, Mexico, and covers an area of about 6.3 million ha.The elevation above sea level varies between 363 and 3190 m.The prevailing climate in the area is mostly temperate, with rainfall in summer.The annual average precipitation ranges between 443 and 1452 mm and the annual average temperature between 8 and 26 • C [59].The predominant type of forest is uneven aged pine-oak forest, often mixed with Arbutus spp., Juniperus spp., Pseudotsuga menziesii and other tree species [56].According to Wehenkel et al. [60], P. cooperi typically grows on very strongly acidic soils (pH 4.5 ± 0.2 (standard deviation (SD))) and is mostly found on areas with a Julian date of the last freezing date of spring (Sday) of 142 ± 11 (SD) days, while the P. durangensis is located on strongly to moderately acidic soils (pH 5.3 ± 0.3 (SD)) in areas with a Sday of 146 ± 26 (SD) days.
A total of 595 P. cooperi trees and 700 P. durangensis trees were subjectively selected to represent the existing range of site qualities and densities in the study area.Diameter at breast height (1.3 m above ground level) (D, in cm) was measured to the nearest 0.1 cm in each tree.In each tree, the height to the base of the live crown (HBLC, m), defined as the point on the stem of the lowest live branch above which there were at least two consecutive live branches, and the largest crown diameter (D crown , m, two measures taken at right angles) were also measured to the nearest 0.1 m in each tree.The trees were later felled, leaving stumps of average height 0.1 m, and total tree height (H, in m) was measured to the nearest 0.01 m.The trees were cut into logs as follows.The first two logs were of constant length 0.3 m and the third log was of variable length, because the upper diameter coincided with D. Subsequent logs were 2.5 m long.All coarse tree branches, defined as those with a diameter at base of trunk greater or equal to 5 cm, were also sectioned at variable length.Branches of diameter <5 cm were not included in the sampling because their contribution to the total volume is often considered negligible [61] and are generally not considered worth being utilized for economic and biomass quality considerations [9], in addition to potential impacts on soil nutrients [16].Two perpendicular diameters over and under bark were measured in each cross section (of trunk or branch) to the nearest 0.1 cm and then arithmetically averaged.Log volumes (m 3 ) were calculated using Smalian's formula.The top section was treated as a cone.Total stem and branch volumes under and above bark were obtained by summing the volume of the sections and the volume of the tops.Summary statistics including number of observations, mean, minimum and maximum values and standard deviation for the main tree variables of both species under study are shown in Table 1.A plot of relative height against relative diameter for both species is shown in Figure 1.
branch) to the nearest 0.1 cm and then arithmetically averaged.Log volumes (m 3 ) were calculated using Smalian's formula.The top section was treated as a cone.Total stem and branch volumes under and above bark were obtained by summing the volume of the sections and the volume of the tops.Summary statistics including number of observations, mean, minimum and maximum values and standard deviation for the main tree variables of both species under study are shown in Table 1.A plot of relative height against relative diameter for both species is shown in Figure 1.

Compatible System Fitting
The system is based on the compatible model proposed by Fang et al. [33].This model assumes that the stem has three sections, each with a different form factor.
The system comprises five equations: The first (Equation ( 1)) is a stem taper function that estimates diameter over bark (d ob , cm) to a height limit (h, m).Integration of this equation to any height limit (h) produces the merchantable volume over bark Equation (2), and integration of Equation (1) over total tree height produces the stem volume over bark Equation (3).
(1) Over bark taper function: where q = h/H and I 1 = 1 if p 1 ≤ q ≤ p 2 ; 0 otherwise I 2 = 1 if p 2 < q ≤ 1; 0 otherwise p 1 and p 2 are relative heights from ground level where the two inflection points are assumed by the model to occur, dividing the stem in three sections Forests 2017, 8, 417 where k is π/40,000, a metric constant for converting from diameter squared in cm 2 to cross-section area in m 2 ; b i is the form factor of tree section i; D is the diameter at breast height over bark (1.3 m above ground, cm); d ob is the top diameter over bark at height h (cm); H is the total tree height (m); h is the height above ground to top diameter d (m); h st is the stump height (m) and a 0 -a 2 , b 1 -b 3, p 1 and p 2 are parameters to be estimated.
(2) Equation for estimating merchantable volume over bark (3) Equation for estimating stem volume over bark where v ob is the merchantable volume over bark (m 3 ), i.e., the volume from stump to the point where diameter = d, and V ob is the total stem volume over bark (m 3 ).Diameter under bark (d ub ) was also estimated using the model described by Fang et al. [33] (Equation ( 4)).
The merchantable volume under bark (v ub ) was therefore effectively obtained by integrating Equation (4) to a height limit (h) (Equation ( 5)), and bark thickness (bt) was estimated by subtracting the diameter under bark (d ub ) from the diameter over bark Integration of Equation ( 4) over total tree height produces the stem volume under bark (V ub , Equation ( 6)).
(4) Under bark taper function where Equations ( 1) and ( 4) differ in relation to the parameters affecting the compatible volume equation (a 0 -a 3 and e 0 -e 3 ), while the form factors (b 1 -b 3 ) and the inflection points (p 1 and p 2 ) are common to both equations.
(5) Equation for estimating merchantable volume under bark (5) (6) Equation for estimating stem volume under bark In an initial analysis, allometric models were fitted to estimate the branch volume (V branches ) for both species.The following different tree variables were combined: diameter at breast height (D), total height (H), crown diameter (D crown , m), crown length (cl, m) and crown ratio (cr), defined as the ratio between crown length and total height.The best results were obtained for both species with the diameter at breast height and the crown ratio.Therefore, two new equations were included in the system to estimate branch volume (V branches ) and total tree volume (V total ) as the sum of stem volume over bark (V ob ) and branch volume (V branches ).

Model Fitting
The system comprises eight equations including endogenous and exogenous variables.The endogenous variables are included on the left-hand side of the equations (d ob , V ob , v ob , bt, v ub , V bark , V branches and V total ) and they are assumed to be determined by the model structure.The remaining variables (D, H, h, h st and cr) are exogenous, i.e., they are independent variables.An endogenous variable on one equation of the system can also appear on the right-hand side of the other equation.
As the compatibility between the different equations does not depend on the parameter estimation process, three different methods can be used to fit the system: (i) estimation of the parameters of the taper function (Equations ( 1) and ( 4) for diameter over and under bark, respectively), estimation of the branch volume (Equation ( 7)) and recovery of the implied volume Equations ( 2), ( 3), ( 5), ( 6) and ( 8); (ii) estimation of the parameters of some of the volume equations (Equations ( 2) and (3) and Equations ( 5) and ( 6) for diameter over and under bark, respectively, and Equation ( 7)), substitution of the estimated parameters in the system, and fitting the remaining parameters; or (iii) estimation of all the parameters of the system simultaneously.
With options (i) and (ii) it is easier to achieve convergence of the parameter estimates and thus obtain the best estimate of diameter at a certain height or of volume equation, depending on which equation is prioritized; however, this may increase the bias and the standard error of the other equation [39].Option (iii) reduces squared error for the total system, i.e., it simultaneously minimizes diameter at different heights and volume prediction errors.The forest manager must select the fitting option by considering whether the system will mainly be used to estimate total volume and then to estimate volumes in size assortments (or vice versa) or will be used for a mixture of these.In this study, we selected the simultaneous fitting option for diameter over and under bark and option (ii) for stem volumes.
In order to ensure that predictions of d ub (Equation ( 4)) are lower than d ov (Equation ( 1)) for the ordinary range of tree conditions, a restriction was included in the fitting process.Provided that Equations ( 1) and ( 4) differ only in relation to the parameters affecting the compatible volume equation (a 0 -a 3 and e 0 -e 3 ), the logarithm of the ratio between Equations ( 3) and ( 6) was restricted to be greater than 0, considering maximum values of diameter at breast height (D) and tree height (H) of 200 cm and 50 m, respectively.log For simultaneous fitting of the system, the number of observations of the endogenous variables must be equal.However, there is more than one diameter observation for each tree but only one observation for total stem volume.To solve this problem, we created a special structure for the data: the total volume of the tree was assigned to each diameter observation on the same tree.The inverse of the number of observations in each tree (n i ) was then used to weight the total volume in the fitting process.

Autocorrelation and Heteroscedasticity
Because the data includes multiple observations for each tree (i.e., it is hierarchical data), autocorrelation within the residuals for each individual might be expected, which would violate the assumption of independent error terms.To overcome possible autocorrelation, the error was modelled using a continuous autoregressive error structure (CAR(x)), which accounted for the distance between measurements.The CAR(x) error structure was programmed to enable dynamic updating of the residuals.The Durbin-Watson (DW) test was used to check for the presence of autocorrelation.Based on a preliminary inspection, some problems of heteroscedasticity related to the use of the allometric Equations (( 7) and ( 8)) were identified, which would violate the general assumptions of least squares regression.In order to account for this, we used weighted regression, by weighting each observation by the inverse of its variance (σ i 2 ).Although this variance is unknown, it is commonly assumed that the variance of the error of the ith individual can be modelled as a power function of a subset of the independent variables, i.e., σ i 2 = (X i ) r [62].The exponential term r was optimized to yield the most homogeneous studentized residual plot.
The MODEL procedure of SAS/ETS ® [63] was used to include the modified k-order autoregressive error structure and the weighting factor for heteroscedasticity 1/(X i ) k and difference in number of observations.Because the data used to fit the models were artificially modified (by weighting), the approximate standard errors of the coefficients were also affected.The following expression was used to calculate the standard errors [64]: where N weight is the number of observations used to fit the models, N true = N weight ∑ i=1 w i , w i is the weighting factor used to modify the data base and p is the number of model parameters.
The residuals were numerically and graphically analyzed to evaluate the goodness of fit by using the coefficient of determination for nonlinear regression (R 2 ) as a statistical criterion (see Ryan 1997, pp.419 and 424).R 2 is defined as the square correlation coefficient between the measured and estimated values (r Y i Ŷi ) and the root mean square error (RMSE), which can be defined as follows: where Y i and Ŷi are the observed and estimated values of the dependent variable, respectively, n is the total number of observations used to fit the model, and p is the number of model parameters.

Effect of Crown Variables on Stem Form Variation
To analyze the effect of crown metrics on stem form, Equations ( 1), ( 2), ( 4) and ( 5) were simultaneously fitted for each tree individually but considering only some parameters as specific for each tree.The parameters associated with the total stem volume equations (a 0 -a 2 for stem volume over bark and e 0 -e 2 for stem volume under bark) were common for all the trees of each species, and their values were the estimates obtained in the previous fitting of the system to the complete database for each species.The inflexion points (p 1 and p 2 ) and the form factors for each stem section (b 1 -b 3 ) were considered as specific for each tree and, therefore, estimated in the new simultaneous fitting for each tree.Different variables related to crown architecture were calculated for each tree: the crown ratio (cr); the ratio between crown diameter and crown length and the slenderness coefficient (s).Pearson's correlation coefficients between tree-specific parameters and crown metrics were calculated and when strong correlation was observed, linear models were fitted to relate the tree-specific parameters with the crown metrics.Finally, these relationships where included in Equations ( 1), ( 2), ( 4) and ( 5) to generalize the stem taper models including crown characteristics.

Results and Discussion
The system was initially fitted without expanding the error term to account for autocorrelation.As expected, because of the hierarchical nature of the data, the values of the Durbin-Watson test for the taper functions (Equations ( 1) and ( 4)) and for the merchantable volumes (Equations ( 2) and ( 5)) of both species and the graphical inspection indicated that the residual trend could be described as a function of lag1-and lag2-residuals within the same tree.The error term was therefore expanded as a continuous autoregressive process, and a second order autoregressive structure was sufficient to eliminate the autocorrelation of d ov and d uv predictions in the system for both species (Durbin-Watson values close to 2 for the four equations and for both species: Table 2).Scatter plots of residuals against predicted branch and whole tree volumes (Equations ( 7) and ( 8)) showed an increasing variance in residuals as the diameter values increased.Therefore, weighting factors of D 4.95 and D 4.65 were included in Equations ( 7) and ( 8) respectively, to correct the heteroscedasticity.Inclusion of these factors allowed minimum variance estimates and reliable prediction intervals to be obtained.The goodness of fit statistics of the simultaneous fit of the variables of study are shown in Table 2 for both species under study.
The taper equations predicted diameter reduction and associated volume over and under bark with high precision, with R 2 values of 0.98-0.99 for both species.High R 2 values were also obtained for the simultaneous prediction of whole tree volume, close to 0.98 for both species under study.The RMSE values were low, approximately 2.1-2.2cm for d ob and d ub, 0.02-0.03m 3 for v ob, v ub , 0.05-0.06m 3 for V branches and 0.15-0.16m 3 for V total for both species.These values are similar to the R 2 values of 0.94-0.95obtained by Corral-Rivas et al. [56] in a study performed in the forest region of El Salto, Durango for predicting taper and total volume for the same two species, and with similar RMSE values.The precision is also similar to that obtained by Quiñonez-Barraza et al. [58] also with the Fang et al. [33] taper equation-for several pine species in Mexico.As commonly found for branch Forests 2017, 8, 417 9 of 18 components, the precision was lower, with R 2 values of respectively 0.65 and 0.7 for P. cooperi and P. durangensis.Similar R 2 values for branch biomass prediction are commonly reported in the literature (e.g., [11][12][13]49,[64][65][66]).Tree branches are widely acknowledged to be difficult to model [67].Moreover, modelling branch volume is particularly challenging because of the variety of sizes and shapes present in tree crowns and because tree crown attributes are expected to vary widely between trees and between locations [3,68].
Scatter plots of the unweighted residuals against diameter classes are shown in Figures 2 and 3.For most of the volume variables in both species, a narrower range in the unweighted residuals could be observed for the higher diameter classes (above 60-70 cm).We speculate this could be due to a limited number of trees sampled in this large diameter range, which are indeed very scarce in harvested forest stands in the area of study.Although crown ratio was considered in predicting branch volume, branch volume was slightly underestimated in the largest diameter classes for both species.Very mature trees tend to have lower branch and leaf allocation (e.g., [13,69], together with some possible crown senescence.Caution must therefore be taken in using the derived equations in very large trees of diameters above 60-70 cm.
The estimated values of the parameters for the simultaneous system fitted are shown in Table 3.
Forests 2017, 8, 417 9 of 17 branch volume, branch volume was slightly underestimated in the largest diameter classes for both species.Very mature trees tend to have lower branch and leaf allocation (e.g., [13,69], together with some possible crown senescence.Caution must therefore be taken in using the derived equations in very large trees of diameters above 60-70 cm.The estimated values of the parameters for the simultaneous system fitted are shown in Table 3. The p1 parameter for P. cooperi suggests a first inflection point at 4.6% of the total tree height for this species, similar to the values of 4.1 to 4.7% reported by Quiñonez-Barraza et al. [58] for several pine species in Mexico, and to the value of 4.1% reported by Corral-Rivas et al. [56] for both of the pine species considered here.In other studies, use of the same taper function yielded higher first inflection points for other pine species: e.g., Castedo-Dorado et al. [70] reported a p1 value at approximately 6% for Pinus radiata plantations in NW Spain; Diéguez-Aranda et al. [39] reported a species; e.g., Quiñonez-Barraza et al. [58] reported values of approximately 60% of the total height for P. leiophylla in Mexico; Diéguez-Aranda [39] reported a value of 61% for Scots pine plantations in northwestern Spain; and Fang et al. [33] reported values of respectively 54% and 57% for P. elliottii and P. taeda, in plantations in the USA.These differences in the inflection point seem to confirm the importance of species-specific fitting to account for differences in the diameter profile of different pine species.The form factors for the three segments (lowest to highest) as found by the ratio bi/k were 0.086, 0.550 and 0.395 for P. cooperi and 0.068, 0.538 and 0.377 for P. durangensis.As expected, the form factors are smallest at bottom, largest in the middle and moderate at the top, and the values for the middle and the top stem are very similar to those of respectively a paraboloid and a cone (0.500 and 0.333).The form factors observed are similar to those obtained for other pine species in Mexico [56,58] and for other pine species in different study areas (e.g., [33,39,70]).
To analyze the effect of crown metrics on stem form, Equations ( 1), ( 2), ( 4) and ( 5) were fitted simultaneously for each tree individually but fixing the values of the parameters associated with the total stem volume equations (a0, a1, a2 and e0, e1 and e2) obtained when the complete database is used to fit the equations system.The mean, maximum, minimum and standard deviation of the estimated values of the inflexion points (p1 and p2) and the form factors (b1, b2 and b3) are shown in Table 4.The p 1 parameter for P. cooperi suggests a first inflection point at 4.6% of the total tree height for this species, similar to the values of 4.1 to 4.7% reported by Quiñonez-Barraza et al. [58] for several pine species in Mexico, and to the value of 4.1% reported by Corral-Rivas et al. [56] for both of the pine species considered here.In other studies, use of the same taper function yielded higher first inflection points for other pine species: e.g., Castedo-Dorado et al. [70] reported a p 1 value at approximately 6% for Pinus radiata plantations in NW Spain; Diéguez-Aranda et al. [39] reported a value of 10% for Pinus sylvestris in Spain and Corral-Rivas et al. [56] reported values of 7-8% for other Mexican pine species such as P. englemanii, P. leiophilla and P. teocote.In this study, the first inflection point was lower for P. durangensis than the aforementioned species and P. cooperi, at approximately 2% of the total height.
The second inflection point, defined by parameter p 2 , was obtained at respectively 71% and 74% of the total tree height for P. cooperi and P. durangensis.These values are similar to those obtained by Corral-Rivas et al. [56], for the same species in Mexico, and by Quiñonez-Barraza et al. [58], for P. arizonica, P. ayacahuite (P.strobiformis) and P. teocote, with the same taper function.
However, the same taper function also yielded earlier second inflection points for other pine species; e.g., Quiñonez-Barraza et al. [58] reported values of approximately 60% of the total height for P. leiophylla in Mexico; Diéguez-Aranda [39] reported a value of 61% for Scots pine plantations in northwestern Spain; and Fang et al. [33] reported values of respectively 54% and 57% for P. elliottii and P. taeda, in plantations in the USA.These differences in the inflection point seem to confirm the importance of species-specific fitting to account for differences in the diameter profile of different pine species.
The form factors for the three segments (lowest to highest) as found by the ratio b i /k were 0.086, 0.550 and 0.395 for P. cooperi and 0.068, 0.538 and 0.377 for P. durangensis.As expected, the form factors are smallest at bottom, largest in the middle and moderate at the top, and the values for the middle and the top stem are very similar to those of respectively a paraboloid and a cone (0.500 and 0.333).The form factors observed are similar to those obtained for other pine species in Mexico [56,58] and for other pine species in different study areas (e.g., [33,39,70]).
To analyze the effect of crown metrics on stem form, Equations ( 1), ( 2), ( 4) and ( 5) were fitted simultaneously for each tree individually but fixing the values of the parameters associated with the total stem volume equations (a 0 , a 1 , a 2 and e 0 , e 1 and e 2 ) obtained when the complete database is used to fit the equations system.The mean, maximum, minimum and standard deviation of the estimated values of the inflexion points (p 1 and p 2 ) and the form factors (b 1 , b 2 and b 3 ) are shown in Table 4. Table 4. Summary statistics of the estimated parameters for simultaneous fitting of Equations ( 1), ( 2), (4) and ( 5) for each tree individually.The Pearson's correlation coefficients between the tree-specific parameters and the tree and crown variables only indicated a strong and significant linear relationship (p < 0.001) between the form factor of the upper stem segment (b 3 ) and two variables: the stem slenderness and the crown ratio.

Mean
Similar results were found by Garber and Maguire [71]; these authors showed that the inclusion of functions of diameter-height ratio were important variables for distinguishing among tree forms in plots of varying spacing and stand structure for Abies grandis Dougl.ex D. Don, Pinus contorta Dougl.ex Loud., and Pinus ponderosa Dougl.ex Laws.In that study, crown ratio added also significant predictive power to the taper function for P. contorta; however, this variable was dropped from the final model.Other previous research also established the desirability of incorporating some stem form surrogate into taper models, typically a diameter-height ratio, height-diameter ratio, or crown ratio (e.g., [24,26,72,73]).
The equations system was refitted for the complete database of each species including linear relationships of the form factor b 3 with these two variables; however, the parameters associated to these variables were not significant and the generalization of the taper equations by including tree or crown variables was dismissed.The reason could be the simplified stand conditions of these conifer species where crown attributes may not add additional information into the model; moreover, it should be taken into account that the taper functions fitted without including crown metrics explained more than 97% of the observed variability.
The models obtained enable prediction of bark volume and bark thickness.Bark thickness increased with increasing tree size (Figure 4, upper figures) in both species, in accordance with studies reporting a positive correlation between bark thickness and diameter at breast height in other conifer species (e.g., [74,75]).For the larger diameter classes, P. cooperi bark was thicker than P. durangensis bark.Average percentages of bark volume relative to stem were respectively 17.2% and 19.8% for P. cooperi and P. durangensis, ranging from approximately 10 to 30% for both species.Figure 4 (lower figures) shows the percentage of volume for tree diameter classes.As expected (e.g., [76]), lower values of relative bark volume were found as diameter at breast height increased.Lower percentage values were found for P. durangensis.Such variation is an important criterion in considering the target diameter in uneven-aged management of pinewoods, with important implications for log merchandising (e.g., [77]).
4 (lower figures) shows the percentage of volume for tree diameter classes.As expected (e.g., [76]), lower values of relative bark volume were found as diameter at breast height increased.Lower percentage values were found for P. durangensis.Such variation is an important criterion in considering the target diameter in uneven-aged management of pinewoods, with important implications for log merchandising (e.g., [77]).
The proposed equations enable bark thickness to be modelled across tree height.Figure 5 represents the variation in predicted bark thickness along the stem heights hi for two diameter classes (i.e., above and below 35 cm) for both species.As expected, a decreasing pattern of bark thickness with height is predicted for both species, with the thinnest values at the tree top for both diameter classes of the two species of study.The P. cooperi bark was thicker than the P. durangensis bark for both diameter classes, and this difference was more marked for the smaller diameter (<35 cm) trees.The thickness of bark relative to diameter variations along the stem is known to be related to the ability of pines to resist fires of different intensity (e.g., [78]).Therefore, P. cooperi may potentially be more resistant to intense fires.Unlike previous studies that focused on the modelling of specific volume components such as stem [56] or bark [79], the system of equations presented here enables compatible prediction of stem, bark and branches volume, thus providing a more biologically consistent understanding of the mechanisms of total tree volume allocation to different fractions.This may help to improve monitoring of the carbon dynamics of these ecosystems.Compatible systems that help clarify the relationships between branch and bark to bole volume could improve forest carbon inventories by reducing error in whole-tree volume estimation [5].
In addition, precise estimation of the volume of bark and branches, combined with speciesspecific piling coefficients for these fractions (e.g., [16,80], may be of great value for planning harvesting and estimating the cost of transporting whole trees and residual biomass (e.g., [14,81,82]  The proposed equations enable bark thickness to be modelled across tree height.Figure 5 represents the variation in predicted bark thickness along the stem heights hi for two diameter classes (i.e., above and below 35 cm) for both species.As expected, a decreasing pattern of bark thickness with height is predicted for both species, with the thinnest values at the tree top for both diameter classes of the two species of study.The P. cooperi bark was thicker than the P. durangensis bark for both diameter classes, and this difference was more marked for the smaller diameter (<35 cm) trees.The thickness of bark relative to diameter variations along the stem is known to be related to the ability of pines to resist fires of different intensity (e.g., [78]).Therefore, P. cooperi may potentially be more resistant to intense fires.
Forests 2017, 8, 417 13 of 17 future studies may expand this approach by simultaneously quantifying the density of bole and branches [5,84,85], to yield a compatible, integral understanding of carbon and biomass allocation and volume distribution in different tree fractions.

Conclusions
A compatible system for simultaneously predicting tree taper, stem volume over and under bark, branch volume and total tree volume was developed for two pine species grown commercially for timber in the Sierra Madre Occidental, Durango, Mexico.The observed variability explained by the compatible system was more than 97% for over and under bark diameters, close to 99% for merchantable stem volumes, more than 64% for volume of branches and more than 97% for total tree volume.Diameter-height ratio and crown ratio were significant and apparently adequate variables for accounting for tree forms variability, although these variables were not finally included into the taper functions.Use of the compatible system provides a better understanding of the volume distribution in different fractions, thus potentially improving carbon accounting in these ecosystems.These models may also be useful for optimizing harvesting, transportation and storage logistics and for sustainable planning of forest operations involving the use of residual forest biomass or whole trees for bioenergy production.Future work might explore the compatibility of volume and biomass predictions for a better understanding of the mechanisms of carbon and biomass allocation and volume distribution in these ecosystems.
Acknowledgments: The data used were collected as part of a program fomenting regional social planning and organization and forest development, funded by the National Forestry Commission of Mexico.
Author Contributions: Study conception, experimental design and field data collection: Jose Javier Corral-Rivas; Analysis and interpretation of data: Juan Gabriel Álvarez-González, Ana Daría Ruiz-González, Jose Javier Corral-Rivas, Daniel Jose Vega-Nieva; Drafting of manuscript: Daniel Jose Vega-Nieva, Roque Rodríguez-Soalleiro, Carlos Antonio López-Sánchez, Christian Wehenkel, Benedicto Vargas-Larreta.Unlike previous studies that focused on the modelling of specific volume components such as stem [56] or bark [79], the system of equations presented here enables compatible prediction of stem, bark and branches volume, thus providing a more biologically consistent understanding of the mechanisms of total tree volume allocation to different fractions.This may help to improve monitoring of the carbon dynamics of these ecosystems.Compatible systems that help clarify the relationships between branch and bark to bole volume could improve forest carbon inventories by reducing error in whole-tree volume estimation [5].
In addition, precise estimation of the volume of bark and branches, combined with species-specific piling coefficients for these fractions (e.g., [16,80], may be of great value for planning harvesting and estimating the cost of transporting whole trees and residual biomass (e.g., [14,81,82] as well as for dimensioning volume requisites for storage of the residual biomass in the biomass processing plants (e.g., [16]).
Future research should focus on the compatibility between volume and biomass estimates for different tree fractions.As in studies that have simultaneously predicted tree bole biomass and volume (e.g., [49,83]), sometimes by considering the vertical changes in wood density (e.g., [11]), future studies may expand this approach by simultaneously quantifying the density of bole and branches [5,84,85], to yield a compatible, integral understanding of carbon and biomass allocation and volume distribution in different tree fractions.

Conclusions
A compatible system for simultaneously predicting tree taper, stem volume over and under bark, branch volume and total tree volume was developed for two pine species grown commercially for timber in the Sierra Madre Occidental, Durango, Mexico.The observed variability explained by the compatible system was more than 97% for over and under bark diameters, close to 99% for merchantable stem volumes, more than 64% for volume of branches and more than 97% for total tree volume.Diameter-height ratio and crown ratio were significant and apparently adequate variables for accounting for tree forms variability, although these variables were not finally included into the taper functions.Use of the compatible system provides a better understanding of the volume distribution in different fractions, thus potentially improving carbon accounting in these ecosystems.These models may also be useful for optimizing harvesting, transportation and storage logistics and for sustainable planning of forest operations involving the use of residual forest biomass or whole trees for bioenergy production.Future work might explore the compatibility of volume and biomass predictions for a better understanding of the mechanisms of carbon and biomass allocation and volume distribution in these ecosystems.

Figure 1 .Figure 1 .
Figure 1.Scatter plots of relative diameter (d/D) against relative height (h/H) for the Pinus cooperi trees (a) and for the Pinus durangensis trees (b) under study.

Figure 2 .
Figure 2. Box-and-whisker plots of the unweighted residuals for the variables of study against diameter class (cm) for Pinus cooperi.The boxes represent the interquartile ranges.The upper and lower ends of the whiskers represent the maximum and minimum height prediction errors respectively.where: dob: diameter over bark (cm); dub: diameter under bark (cm); vob: merchantable volume over bark (m 3 ); vub: merchantable volume under bark (m 3 ); Vbranches: coarse branches (>5 cm) volume (m 3 ); Vtotal: total tree volume (m 3 ).

Figure 2 .
Figure 2. Box-and-whisker plots of the unweighted residuals for the variables of study against diameter class (cm) for Pinus cooperi.The boxes represent the interquartile ranges.The upper and lower ends of the whiskers represent the maximum and minimum height prediction errors respectively.where: d ob : diameter over bark (cm); d ub : diameter under bark (cm); v ob : merchantable volume over bark (m 3 ); v ub : merchantable volume under bark (m 3 ); V branches : coarse branches (>5 cm) volume (m 3 ); V total : total tree volume (m 3 ).

Figure 3 .
Figure 3. Box-and-whisker plots of the unweighted residuals for the variables of study against diameter class (cm) for Pinus durangensis.The boxes represent the interquartile ranges.The upper and lower ends of the whiskers represent the maximum and minimum height prediction errors respectively.where: dob: diameter over bark (cm); dub: diameter under bark (cm); vob: merchantable volume over bark (m 3 ); vub: merchantable volume under bark (m 3 ); Vbranches: coarse branches (>5 cm) volume (m 3 ); Vtotal: total tree volume (m 3 ).

Figure 3 .
Figure 3. Box-and-whisker plots of the unweighted residuals for the variables of study against diameter class (cm) for Pinus durangensis.The boxes represent the interquartile ranges.The upper and lower ends of the whiskers represent the maximum and minimum height prediction errors respectively.where: d ob : diameter over bark (cm); d ub : diameter under bark (cm); v ob : merchantable volume over bark (m 3 ); v ub : merchantable volume under bark (m 3 ); V branches : coarse branches (>5 cm) volume (m 3 ); V total : total tree volume (m 3 ).

Figure 4 .
Figure 4. Observed bark thickness and percentage of bark volume for the diameter classes of the studied trees of Pinus cooperi (left figures) and Pinus durangensis (right figures).BarkThickness: thickness of bark in cm; dclass: normal diameter class (cm); VolBarkPerc: percentage of volume of bark to volume of stem under bark (%).

Figure 4 .
Figure 4. Observed bark thickness and percentage of bark volume for the diameter classes of the studied trees of Pinus cooperi (left figures) and Pinus durangensis (right figures).BarkThickness: thickness of bark in cm; dclass: normal diameter class (cm); VolBarkPerc: percentage of volume of bark to volume of stem under bark (%).

Figure 5 .
Figure 5. Predicted bark thickness (cm) across tree height (m) for two diameter (D) classes above and below 35 cm for trees of both species under study.hi: measurement height (m); D: diameter at breast height (cm): sp1: Pinus cooperi; sp2: Pinus durangensis.

Figure 5 .
Figure 5. Predicted bark thickness (cm) across tree height (m) for two diameter (D) classes above and below 35 cm for trees of both species under study.hi: measurement height (m); D: diameter at breast height (cm): sp1: Pinus cooperi; sp2: Pinus durangensis.

Table 1 .
Summary statistics of the data on the tree species under study.

Table 2 .
Goodness of fit statistics for simultaneous fitting of the equation system.

Table 3 .
Parameters for simultaneous fitting of the equation system.