A Preliminary System of Equations for Predicting Merchantable Whole-Tree Volume for the Decurrent Non-native Quercus rubra L. Grown in Navarra (Northern Spain)

1. Introduction

Estimation of tree wood volume is one of the most important tasks for forest managers. Tree volume is usually considered to equal stem volume, and taper equations are used when product demand is expressed in terms of tree dimensions for different industrial uses. Taper equations predict the variation in diameter along the stem, and they therefore characterize stem form (Clutter et al. 1983; Avery and Burkhart 2002). Integration of the taper equation from the ground to any height provides an estimate of the wood volume up to that height. However, in contrast to the almost always excurrent form in softwoods, hardwoods often have decurrent form (Assmann 1970; Burkhart and Tomé 2012); and not all trees have a well-defined stem and they may be forked or lack apical dominance. In the development of taper equations, trees that differ from the idealized excurrent form are generally excluded (MacFarlane and Weiskittel 2016). However, some authors have incorporated different stem forms in stem taper models; thus e.g. Adu-Bredu et al. (2008) considered three general forms of Tectona grandis L.f. in West Africa, in order to take into account the effect of forks in the taper equation: zero-forked trees, one-forked trees and two-forked trees. For forked trees, the height(s) of fork(s) must be measured by the forest manager. MacFarlane and Weiskittel (2016) developed a generalized stem taper which can be calibrated by adding random-effects parameters depending on the species (22 species were considered) and merchantable form types (10 types were considered, including the combination of the presence/absence of merchantable stem/branches with the destination of saw/pulp).

In the other hand, some branches may be of sufficiently large diameter to be merchantable (Figure 1) with a main destination of pulp products (MacFarlane and Weiskittel 2016), and their volume should be known. However, most research is focused on modelling stem volume, and estimation of branch volume has received less attention (MacFarlane 2011; Ver Planck and MacFarlane 2014; Corral-Rivas et al. 2017). Quantification of tree branches has typically considered branch biomass, as quantifying branch volume is challenging owing to the large numbers of pieces of different sizes and shapes (MacFarlane 2011).

The northern red oak (Quercus rubra L.) was introduced to northern Spain for productive purposes in the 20th century, and Navarra region currently has the largest stocks of this species (Pemán García et al. 2013). Because of the decurrent form, the Regional Forest Service, in addition to diameter at breast height (D, measured 1.3 m above the ground) and total tree height (H), usually also measures the merchantable height (h_m, Figure 1). Merchantable height is the subjective height up to which the stem has a form similar to excurrent form. When only taper equations are developed using D and H, the following are assumed: (i) excurrent form occurs and the stem is considered by the minimum top diameter fixed by the forest manager (d_min, Figure 1 right); and (ii) merchantable whole-tree volume (V_t) would be equal to stem volume, and V_t could therefore be underestimated. It has been found by Regional Forest Service that V_t is generally underestimated in trees classified as “poorly formed” because they have a large branch volume. Stem volume equations have already been developed for different stem forms and the main forest species considered in the Spanish National Forest Inventory (SNFI) (MITECO 2023). However, the SNFI approach does not classify stem products or consider the branch volume.

The overall objective of the present study was to develop a preliminary system of equations that allows estimation of the merchantable whole-tree volume and classification of wood products for Q. rubra growing in Navarra. The specific objectives were as follows: (i) to develop a taper equation that enables classification of stem products; (ii) to develop an equation that estimates the merchantable whole-tree volume (V_t), for both stem and branches; and (iii) to determine the compatibility of both equations and their applicability in routine forest inventories in Navarra. The nomenclature used in this study is summarized in Table 1.

2. Materials and Methods

2.1. Data

The study was conducted in the region of Navarra (North of Spain). A total of 30 oak trees were destructively sampled in 2001 in Etxalar in order to develop a taper equation for the Navarra Government. In 2021, sampling of 34 trees was carried out in Altsasu. Diameter at breast height (D -cm-, 1.3 m above ground level) was measured to the nearest 0.1 cm for each tree. Subsequently, the trees were felled, resulting in stumps of average height of 0.10 m. The maximum length of the felled trees was measured to the nearest 0.01 m to determine the total tree height (H, m). The bole was subsequently sectioned into logs at 2 m intervals in the trees from Etxalar and 1 m in those from Altsasu, up to a minimum diameter of 7 cm in the trees which had a well-defined dominant stem or, otherwise, until the merchantable height (h_m) was reached. Two perpendicular diameters over bark were measured in each cross section (at height h -m- from ground level), to the nearest 0.5 cm, and the values were then averaged (d, cm). Moreover, in Altsasu the merchantable branches also were measured at 1-m intervals. For the branches, the same minimum diameter of 7 cm was defined unless branch architecture suggested a larger top diameter. The volume of a section (stem or branches) was calculated by Smalian´s formula.

2.2. Taper Equation

The scatter plot of relative diameter (d/D) against relative height (h/H) was visually examined for the 64 available trees to detect any possible anomalies in the data. Because inconsistent data were detected, a systematic method was taken to identify abnormal data points (Bi 2000), enhancing process efficiency. This involved employing local quadratic nonparametric fitting (LOESS) with a smoothing factor of 0.3, determined through iterative fitting and visual assessment of the smoothed taper curves overlaid on the data. This approach detected 1.28% anomalies from the initial data points, and these outliers were removed before parameterizing the taper models. Some of the outliers could be due to errors in measuring the bole sections or mistakes in transcribing field notes. However, others might be the result of deformations caused by fire damage, large knots, and other types of physical damage, such as partial death, growth deformities, and so on. For the data ultimately utilized, the summary statistics are displayed in Table 2, while Figure 2 illustrates the scatter plot of d/D against h/H.

The taper equations tested include the models developed by Max and Burkhart (1976), Biging (1984), Stud (Daquitaine et al., 1999), Bi (2000), Fang et al. (2000), Sharma and Zhang (2004), and Kozak (2004). Firstly, the models were fitted by Ordinary Nonlinear Least Squares (ONLS) using the Gauss-Newton iterative method (Hartley 1961). The initial values used in the iteration process were obtained from previous studies.

2.3. Merchantable Whole-Tree Volume

Merchantable whole-tree volume equation (V_t) includes stem volume (v_s or V_s) and merchantable branch volume (V_b). We chose to develop a merchantable whole-tree volume equation instead of a merchantable branch volume equation because merchantable branch volume will be estimated by subtraction once stem volume is estimated up to the merchantable height (h_m). All 64 trees were used for development of the taper equation, but only the 34 trees from Altsasu (summary statistics in Table 3) were used to develop the merchantable whole-tree volume equation (V_t) because merchantable branch volumes were calculated.

Linear, power and exponential models were fitted during the development of the V_t equation. The linear models were fitted by linear regression using the least squares method, and the non-linear models were fitted by ONLS using the Gauss-Newton iterative method (Hartley 1961). The initial values used in the iteration process were obtained by linear regression on the logarithmic transformation of the models. This approach is useful for solving many fitting problems in nonlinear models (Draper and Smith 1998). Diameter at breast height (D), total tree height (H), the combined variable D²H, the merchantable height (h_m), and the h_m/H ratio were tested as independent variables. In the present study, the h_m/H ratio ranged from 0.36 to 0.91, with a mean value of 0.72.

Merchantable whole-tree volume (V_t) should be always greater than or equal to stem volume (V_s). However, this restriction was not included in the model formulation. A simulation was carried out to verify that V_t ≥ V_s. A database of D in the range 10-100cm was simulated at intervals of 5 cm and all possible combinations were crossed with a database of H in the range 5-60m, at intervals of 5 m. V_s was selected in the comparison because it is the highest value estimated by the taper equation.

2.4. Model Fitting

Statistical and graphical analyses were conducted to compare the performance of the equations. Plots of residuals against predicted values and plots of observed values against predicted values were examined to evaluate heteroscedasticity, anomalous trends and bias. Moreover, in development of the taper equation, the residuals were evaluated across relative height and D classes (McTague and Weiskittel 2021), both in absolute values of residuals and in relative values ($\left(Y_i-{\widehat{Y}}_i\right)/Y_i$, see notation below). The numerical analysis consisted of comparing three statistical criteria obtained from the residuals (Myers 1990): the coefficient of determination (R²), which measures the amount of the observed variability explained by the model, the mean bias (e), which is the average difference between the values estimated by a model and the actual observed values, and the root mean squared error (RMSE), which provides a measure of the precision of the estimates in the same units as the dependent variable. Despite the recognized limitations attributed to the utilization of R² in nonlinear regression analysis, it is pertinent to underscore that the broader applicability of a comprehensive measure of model adequacy may supersede such constraints (Ryan 1997). The mean bias (e) and RMSE are also expressed as percentages (e (%) and RMSE (%)). The expressions of these statistics are summarized as follows:

where $Y_i,{\widehat{Y}}_i$ and $\bar{Y}$ are the observed, estimated and the mean 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.

In regression analysis, the error terms should be independent, identically distributed, normal, random variables. However, in developing taper equations, multiple observations for each tree are used (i.e. longitudinal data). Thus, it is reasonable to expect that the observations within each tree will be spatially correlated and the assumption of independent error terms will therefore be violated (Burkhart and Tomé 2012). Autocorrelation can be taken into account by modelling the error term. To correct for the autocorrelation between residuals from the same tree, we used a continuous autoregressive error structure (CAR(x)), which accounts for the distance between measurements. This error structure expands the error term to (Zimmerman and Núñez-Antón 2001):

$$e_{ij}={\sum }_{k=1}^{k=x}{I}_k{\rho }_k^{h_{ij}-h_{ij-k}}{e}_{ij-k}+{\epsilon }_{ij}$$

(4)

where e_ij is the jth ordinary residual of the ith individual, e_ij−k is the j-kth ordinary residual of the ith individual, I_k = 1 when j > k and 0 when j ≤ k, ρ_k is the k-order continuous autoregressive parameter to be estimated, and h_ij-h_ij−k is the distance separating the jth from the j-kth observation within each tree i, with h_ij > h_ij−k. In this case ε_ij is an independent normal distributed error term with mean value of zero.

To assess the existence of autocorrelation and determine the appropriate order of the continuous autoregressive error structure (CAR(x)), plots of residuals versus lag-residuals (residuals from previous observations) within each tree were visually inspected. Appropriate fits for the selected model with correlated errors were obtained by including the CAR(x) error structure in the MODEL procedure of SAS/ETS^® software (SAS Institute Inc. 2008), which allows for dynamic updating of the residuals. In developing the total merchantable tree volume equation (V_t), the linear models were fitted using the REG procedure of SAS/STAT^® (SAS Institute Inc. 2009) and the non-linear models were fitted using the NLIN procedure of the same software.

3. Results

3.1. Taper Equation

Although all models tested satisfactorily estimated diameters at different heights, the model of Fang et al. (2000) performed slightly better than the others and/or all of its parameters were significant at the α=5% significance level. This model assumes that the stem has three sections, each with a different form factor:

where d is the diameter (cm) at height h; h is the height above ground (m) to diameter d; D is the diameter at breast height (1.3 m above ground, cm); H is the total tree height (m); h_st is the stump height (m); k is π/40,000, a metric constant for converting from diameter squared in cm² to cross-section area in m²; b_i is the form factor of tree section i and a₀–a₂, b₁–b₃, p₁ and p₂ are parameters to be estimated.

Integration of Equation (5) to any height limit (h) produces the merchantable volume Equation (6), m³. Assuming the tree has an excurrent form, integration of Equation (1) over total tree height (H) produces the stem volume Equation (7), m³.

$$v_s=c_1^2{H}^{k/b_1}\left(b_1{r}_0+\left(I_1+I_2\right)\left(b_2-b_1\right)r_1+I_2\left(b_3-b_2\right){\alpha }_1{r}_2-\beta {\left(1-q\right)}^{k/\beta }{\alpha }_1^{I_1+I_2}{\alpha }_2^{I_2}\right)$$

(6)

$$V_s=a_0{D}^{a_1}{H}^{a_2}$$

(7)

The model of Fang et al. (2000) was initially fitted using ONLS, without expanding the error terms to account for autocorrelation. A trend in residuals as a function of lag-residuals within the same tree was apparent, as expected because of the longitudinal nature of the data used for model fitting (Figure 3, first row). To correct autocorrelation, the error term was modelled with a second order CAR structure, CAR(2), since the first order structure, CAR(1), did not totally correct autocorrelation with more than one lag (Figure 3, third and second rows, respectively). Although accounting for autocorrelation does not improve the predictive ability of the taper equation, it prevents underestimation of the covariance matrix of the parameters, thereby enabling use of the usual statistical tests (West et al. 1984), i.e. it allows valid interpretation. All parameter estimates, including autocorrelation parameters, were significant at the 5% level (Table 4).

The sole purpose of correction for autocorrelation was to obtain unbiased and efficient estimates of the parameters (Parresol and Vissage 1998), and it has no use in practical application. Therefore, to obtain a better idea of the true error expected, for practical use of the taper equation, we calculated the residuals obtained without autocorrelation parameters. The values of the statistical criteria were R² = 0.901, e = 0.369 cm, e(%) = 1.42, RMSE = 3.54 cm and RMSE(%) = 13.7. The plot of residuals against estimated values showed a random pattern of residuals around zero, with no detectable significant trends (Figure 4).

The behaviour of residuals across the diameter at breast height classes was homogeneous both in absolute values and relative values of residuals (Figure 5, second row). The behaviour of residuals in absolute values across the relative height class was also homogeneous (Figure 5, first row, left), but the relative residuals increased with relative height (Figure 5, first row, right).

3.2. Merchantable Whole-Tree Volume

The merchantable whole-tree volume equation (V_t) includes stem volume (v_s) and merchantable branch volume (V_b, with d_min= 7cm). The values of the statistical criteria for the model selected (Equation 8) were R² = 0.861, e = -0.00625 m³, e(%) = -0.389, RMSE = 0.270 m³ and RMSE(%) = 16.8.

$$V_t=t_0\mathrm{e}\mathrm{x}\mathrm{p}\left(t_{1}D\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(t_{2}H\right)$$

(8)

where V_t is the total merchantable tree volume (m³); D is the diameter at breast height (cm); H the total tree height (m); and t₀, t₁ and t₂ are the fitted parameters, all of which were significant at the 5% level (Table 5). The plot of residuals against estimated values for Equation (8) showed a random pattern of residuals around zero, with no detectable significant trends (Figure 6).

3.3. System of Equations

Both the taper equation and merchantable whole-tree volume (V_t) equation have been integrated into the Cubica Navarra software (https://administracionelectronica.navarra.es/CubicacionMadera/), to be used as a system of equations. The input variables are diameter at breast height (D), total tree height (H) and merchantable height (h_m). If h_m is not specified, the system considers a well-defined excurrent form, with the stem top height defined by the minimum top diameter (d_min) selected by the forest manager. For merchantable whole-tree volume the d_min was fixed at 7 cm. The system also classifies products in the stem, with the remaining merchantable volume assumed to be in the branches, estimated by subtracting the stem volume from the merchantable whole-tree volume. Simulations performed with all possible combinations of D = 10-100 cm (at 5 cm intervals) and H = 5-60 m (at 5 m intervals) produced logical results (V_t > V_s) across the whole range of values.

4. Discussion

We developed a merchantable whole-tree volume model system, including stem and merchantable branches, based on the current measurements available from the Navarra Regional Forest Service. In accordance with Ver Planck and MacFarlane (2014), this study did not exclude forked trees or trees of poor form, typically seen in hardwood species, and which would exclude trees with large branches. MacFarlane and Weiskittel (2016) describe some defects in the ideal excurrent stem form which lead to trees or even stands being excluded in studies of stem taper equation development. The objective of the present study was to improve economic assessment of wood products, focused on stem volume because this prioritizes economic profitability while also considering the volume of merchantable branches.

The data used to develop the present system of equations is not widely representative regarding either the spatial distribution or number of trees; moreover, the ranges of diameters (25-58.5 cm) and heights (17.5-31.4 m) did not include small or overmature trees. However, in developing the taper equation, we used 930 observations of stem diameters along the bole, which is more than the minimum of 825 observations reported by Kitikidou and Chatzilazarou (2008) to be needed to parameterize taper equations correctly. We also used 64 trees, which is greater than the minimum 40-85 sample trees needed, and measurements of diameters along the stem were taken at least every 2 m (McTague and Weiskittel 2021). For development of merchantable whole-tree volume equations, it is advisable to use several stands because tree form varies between stands as a result of several factors such as site productivity, density, genetics, stand management and historical biotic or abiotic damage (Gevorkiantz and Olsen 1955; Neilsen and Gerrand 1999; Kerr and Boswell 2001; Gomat et al. 2011; MacFarlane and Weiskittel 2016; McTague and Weiskittel 2021). Taking biomass into account, Montagu et al. (2005) associated the proportion of mass allocated to branches with location and D. Stem form and partitioning of biomass among components have been shown to be affected by density, thinning, pruning and fertilization (Bi et al. 2010). In this study, only two locations were included; therefore, the developed model system is considered preliminary and independent data from other locations will be useful for refitting or validating it. In the taper equation development, we did not consider the discontinuity of diameter produced by the crown basis (Le Moguédec and Dhôte 2012) or the forks (Adu-Bredu et al. 2008). However, a new variable -already measured by the Navarra Regional Forest Service- is considered, the merchantable height (h_m), simplifying consideration of the rest of the merchantable volume (including the wood above h_m) as branches.

In the prediction of merchantable whole-tree volume, inclusion of the available merchantable height (h_m) or the h_m/H ratio as independent variables did not improve the predictive capacity of the model. MacFarlane (2011) found that variables such as crown ratio (CR) and the diameter at the base of the largest branch were good predictors of large branch volume. Crown ratio can be calculated by measuring height to the crown base, in addition to total tree height. Digital hypsometers can be used to measure several heights. However, variables such as diameter at the base of the largest branch or crown diameter are not easily measured. MacFarlane (2011) selected D and CR as predictor variables in developing species-specific large branch volume equations. Mäkelä and Valentine (2006) also indicated that CR is useful for tree allometry prediction. Crown length can be also used (Özçelik and Bal 2013). McTague and Weiskittel (2021) mentioned that the exploration of crown variables is encouraged in the development of taper equations. Moreover, the widely used crown classes (dominant, codominant, intermediate and suppressed trees) could be tested (Sanquetta et al. 2020). The use of emerging technologies such as terrestrial laser scanning (TLS) could improve the development and/or utilization of taper and volume equations (McTague and Weiskittel 2021). Moreover, non-destructive technologies such as ultrasonic and thermal images could be used to detect internal defects (Taskhiri et al. 2020).

Large branches can be defined as those with a significant impact on stem form (MacFarlane 2010) or through a minimum top diameter. MacFarlane (2011) defined large branches in hardwoods in Michigan state (USA) as those of diameter ≥ 9 cm and observed that the proportion of sampled trees with large branches changed as a function of the D class, ranging from 0% in the 10-cm class up to almost 100% in ≥60 cm classes. In our study, minimum top diameter for merchantable branches is 7 cm and only one tree (D= 38 cm) did not have merchantable branch volume (V_b). The V_b/V_t ratio ranged from 0 to 43% (mean 20%), and there was a non-significant positive correlation between V_b/V_t and D (Pearsons r = 0.25, p-value = 0.159). A weak positive correlation between V_b/V_t and D was also observed by MacFarlane (2011) and Ver Planck and MacFarlane (2014). It should also be taken into account that the data used in this study did not include small trees; and Gómez-García et al. (2016) reported that the contribution of each tree biomass component to the total biomass tends to stabilize as the tree diameter increases in Quercus robur L. growing in NW Spain.

Minimum top diameter (d_min) for merchantable volume is not a global standard and includes values of e.g. 5 cm (Rytter 2002), 7 cm (Segura and Kanninnen 2005; Gómez-García et al. 2016), 9 cm (Lambert et al. 2005; MacFarlane 2011) and 10 cm (Jenkins et al. 2003; MacFarlane and Weiskittel 2016). Moreover, regarding the stem form and the merchantability of branches, values other than the minimum top diameter are considered, and our data included values of top diameter per whole-tree (stem and branches) in the range 7.0-10.5 cm. On the other hand, new biomass harvesting methods and industrial uses have increased the demand for previously non-merchantable products (Nicholls et al. 2009), leading to a change in the minimum top diameter.

More informative models of whole-tree volume could be developed by measuring and testing ancillary variables. For instance, it may be possible to develop equations modelling the vertical allocation of branch volume (Ver Planck and MacFarlane 2014), equations with a variable minimum top diameter (Gómez-García et al. 2015) and considering different merchantable form types (MacFarlane and Weiskittel 2016). Such approaches should also include trees with decurrent form. Further research is also necessary to consider variables related to stem form, branch architecture and volume, with the purpose of developing more biologically based models.