Modelling the variation of bark thickness within and between European silver fir (Abies alba Mill.) trees in southwest Germany

Modelling the variation of bark thickness within and between European silver fir (Abies alba... Abstract This study examined bark thickness variability of Silver fir (Abies alba Mill.) in southwest Germany over time and space by comparing a dataset of bark measurements from the 1970s with more recent assessments. Within-tree variability of bark thickness was analysed to estimate the required number of sampling locations per tree for defined accuracy levels. A range of models from the literature that predict bark thickness were compared for their predictive performance and Monte-Carlo simulations were used to estimate the effect of the number of sample trees and plots on the precision of the predictions. In addition, net log volume after bark subtraction was compared for logs of varying lengths to assess the influence of log assortments on calculated sales volume. Results show that several sampling locations are needed per tree and that at least five trees from at least 35 plots should be selected for measurements in the study region. For practical applications, diameter outside bark and breast height diameter are suggested as explanatory variables for models that predict double bark thickness. Additionally, relative tree height and age – and therefore growth rate – significantly improved predictions; however, environmental factors could not explain the variation between stands. Log lengths from 5 to 21 m only slightly influenced bark thickness equations that were fit on measurements at log midpoints. The findings highlight the need to consider bark thickness variability at different levels when developing bark thickness equations. In general, bark thickness was found to be smaller in more recent assessments and this indicates the need to regularly review existing bark equations for their validity. Introduction The diameter inside bark is a crucial measure in forestry that is used to calculate wood volume of logs and trees and to optimize bucking for defined inside-bark diameters. However, tree and log diameters are usually measured outside bark and bark thickness has to be estimated. Incorrect bark thickness or bark volume estimates could lead to inaccurate wood volume estimates in forest inventories, in increment studies or in the log trade (Marshall et al., 2006). Furthermore, the interest in accurate bark thickness estimates is driven by the need to obtain not only accurate inside-bark diameters and volumes, but also accurate estimates of bark volume. The importance of accurate bark thickness estimates has increased with the shift in the commercial relevance of bark from an unwanted residue to a valuable fuel and feedstock for high-value biomaterial (Doruska et al., 2009), such as bark tannin based foams (e.g. Pizzi, 2016). The estimation of available bark biomass is important to assess the potential of such technologies for generating possible additional income for the forestry sector. Additionally, bark volume estimates are needed for biomass estimates, which are becoming more important for quantifying carbon stocks (Temesgen et al., 2015). In this study, we use a common definition of bark that includes all tissues outside the vascular cambium, and comprises the secondary phloem up to the last formed periderm and the rhytidome, which comprises all layers of dead tissue outside the currently active periderm (Martin and Crist, 1970). For many coniferous species, it has been shown that bark thickness can be described well by diameter outside bark at the sampling location and additional variables, as it often correlates with diameter at breast height, total tree height and height of the sampling location within the tree (e.g. Li and Weiskittel, 2011). Several external factors that correlate with the relative bark thickness of trees, i.e. the proportion of the outside-bark diameter that constitutes bark, have been reported. Most of these can be interpreted as a consequence of slower tree growth, which leads to a larger relative bark thickness. A correlation between a low site index or yield class of stands with a larger relative bark thickness was shown for Norway spruce (Picea abies (L.) Karst) (Hoffmann, 1958; Dimitrov, 1976; Schmidt-Vogt, 1986) and Silver fir (Abies alba Mill.) (Božić et al., 2007). For Norway spruce (Laasasenaho et al., 2005, Stängle et al., 2017) and Silver fir (Božić et al., 2007), tree age was shown to have an influence on the diameter-to-bark ratio. Sonmez et al. (2007) also reported a positive effect of tree age on bark thickness for Picea orientalis (L. Link). Contrasting results were reported for the broadleaved Nothofagus pumilio (Poepp. & Endl.) Krasser by Cellini et al. (2012), where relative bark thickness was lower on low productivity sites. For Pinus sylvesteris (L.), the latitude of forest stands – which mainly determines differences in temperature and length of the growing season – could help to explain variation in bark thickness (Wilhelmsson et al., 2002). For Pinus taeda (L.) in the United States, there were no observed effects of silvicultural treatments, such as initial planting density and thinning, on whole-tree bark percentage. Only very intense treatments including fertilization and early-age weed control affected bark volume (Antony et al., 2015). Besides growth conditions, genetics can also determine bark thickness. Provenance has been shown to influence the phenotypic development of bark thickness of Pinus contorta var. latifolia Engelm. (Persson and Downie, 1992) and Pseudotsuga menziesii (Mirb.) Franco (Kohnle et al., 2012). To develop accurate equations for predicting regional bark thickness, the intraspecific within-tree and between-tree variability of bark thickness have to be considered. Bark thickness variation within trees occurs at two levels. First, variation at the disc level around the circumference of the tree can be caused by a non-uniform bark thickness that is correlated with elliptical stem growth (e.g. Niklas, 1999) or a bark surface pattern with fissures, cracks or scales. The second level of within-tree variation occurs between sampling locations along the stem axis. To capture both sources of within-tree variability, several bark thickness measurements are usually taken at each of several sampling locations along the stem. For Norway spruce, variability within and between trees and between sites indicates the importance of multiple sampling points within each tree of a range of selected trees covering a range of sampling sites throughout the study region (Stängle et al., 2016b). Roundwood volume calculation in Germany is performed according to the Huber formula, which is based on a cylinder volume calculation using log length and the rounded diameter inside bark at log midpoint, which is located at exactly halfway along the log’s length (Fonseca, 2005). In the study region, long logs (6–21 m) are mostly measured manually, whereas shorter logs from cut-to-length operations are usually measured at the mill site using electronic scanning devices. Measurement instructions for manual log measurement in Germany require the operator to take two diameter measurements at log midpoint (approximately perpendicularly) and round them to the lower full centimetre. The mean of these rounded values is also rounded down and species-specific fixed integer bark deduction values are subtracted to estimate the inside-bark diameter (Anonymous, 2015). Geospatial patterns of bark thickness within a study region could lead to skewed wood volume estimations if only one bark thickness equation and one set of fixed deduction values is used for the whole region. Additionally, relative bark thickness can vary along the stem and if the applied bark thickness function does not consider the position in the stem, volume calculations could be biased for different assortments. If long logs of up to 21 m length are produced, for example, most log midpoints would be positioned in the lower half of the stem, whereas midpoints of 5-m segments can basically be positioned all along the stem. The integer bark deduction values that are currently in use have been extracted from bark thickness equations that report double bark thickness as a function of mid-diameters of logs of up to 26 m in length (Altherr et al., 1978). However, different bark thickness equations could be useful for different assortments. Furthermore, bark thickness equations should be up-to-date as allometric relationships can change with changes in growth conditions. For Silver fir, growth conditions in the study area have changed in the last century with the average growth rate having increased since the 1990s (Kohnle et al., 2014). Currently, bark thickness equations applied in Germany are based on measurements from the 1970s and probably overestimate bark volume as has already been shown for Norway spruce in southwest Germany (Stängle et al., 2016a). Using Silver fir also from southwest Germany, the objectives of this study were to: estimate sample sizes required to capture bark thickness variation within trees and to examine the effect of number of sample plots and trees per plot on the precision of bark thickness equations; compare and evaluate established model forms for estimating bark thickness and test how geospatial variation could be explained by environmental factors; and estimate how log length and different datasets (measurements from the 1970s and 2010s) affect bark thickness equation coefficients for roundwood volume calculation. Materials and methods Data Bark thickness data were obtained from two sources. Dataset 1, an extensive dataset with 3008 Silver fir trees, was collected by Altherr et al. (1978) in 82 clusters. Several clusters were sampled in the same stands, resulting in 50 sampled stands. Each stand was treated as one plot in the analysis. Dataset 2 consisted of newly acquired measurements on 217 trees in 27 plots. Measurements of dataset 2 were taken in the years 2015 and 2016. Dimensional details of all trees are listed in Table 1. Trees of both datasets were sampled throughout the 35 752 km2 state of Baden-Württemberg in temperate forests with a broad range of site and stand conditions. Elevation of the plots from dataset 2 varied between 306 and 974 m above sea level (Figure 1). Mean annual air temperature varied between 5.8°C and 9.4°C and the annual precipitation was between 868 and 1944 mm. All stands were mixed-species stands with a range of ~20–80 per cent fir trees. Trees from a large range of merchantable diameter classes were selected from regularly harvested trees in selection cuts. Table 1 Characteristics of the measured trees from dataset 1 and 2. Dataset 1 originates from Altherr et al. (1978). Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Note: dbh = diameter at breast height; DBT1.3 = double bark thickness at breast height; SD = standard deviation. Table 1 Characteristics of the measured trees from dataset 1 and 2. Dataset 1 originates from Altherr et al. (1978). Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Note: dbh = diameter at breast height; DBT1.3 = double bark thickness at breast height; SD = standard deviation. Figure 1 View largeDownload slide Distribution of sample plots in the study region and number of sampled trees per plot for dataset 2. The distribution map of Silver fir originates from (EUFORGEN, 2009). Figure 1 View largeDownload slide Distribution of sample plots in the study region and number of sampled trees per plot for dataset 2. The distribution map of Silver fir originates from (EUFORGEN, 2009). The sampling for both datasets was performed as described by Altherr et al. (1974; 1978): trees were felled with chainsaws, delimbed, and measured in the forest before any further log manipulation was performed. Measurement locations were at breast height (1.3 m above ground) and along the tree bole in 2 m increments up to a top diameter of ~10–15 cm. Both diameter and bark thickness were measured twice (approximately perpendicularly) at each location using a calliper and a Swedish bark gauge, respectively. Double bark thickness was calculated as the sum of the two gauge measurements. Additional measurements for the new data were tree age, and total tree height, which was measured on the felled tree. Tree height for dataset 1 and for trees with broken tips in dataset 2 was modelled with mixed B-spline regression describing tree taper using the package TapeR (Kublin and Breidenbach, 2015) in R (R Core Team, 2015). Two predefined bucking patterns were applied to produce virtual logs of different length from both datasets. The first bucking rule aimed at long logs. From each tree one virtual log between 6 and 21 m was generated. According to the second rule, several segments of 5 m length were generated from each tree. Due to convergence problems of the taper functions, only trees with a diameter at breast height (dbh) larger than 18 cm could be used for this analysis. For each log, diameter and bark thickness at midpoint and at the small end was calculated using taper equations based on inside-bark and outside-bark diameter. To evaluate the effect of log length on equation coefficients and roundwood volume, the harvesting statistics of the State Forest Enterprise of Baden-Württemberg (ForstBW) were used. These data, hereafter referred to as dataset 3, included all manually measured logs of length longer than 6 m that were produced in the state forest of Baden-Württemberg in the period 2011–2015. Bark thickness variability at different hierarchical levels Bark thickness measurements within trees and within plots can be expected to be similar (e.g. Li and Weiskittel, 2011). The variance within trees determines the sampling design with the number of required sampling locations per tree and the number of measurements per location. If the variance between groups (such as trees and plots in our case) is large, additional information on group-level predictors could help explain this variance when modelling bark thickness (Schielzeth and Nakagawa, 2013). To quantify the degree of variation that can be accounted for by different grouping levels, we fitted the established bark thickness prediction equation (1), which does not consider any predictor variables at the tree or plot level, but only the diameter outside bark at the measurement location. DBT=β0+β1dob+ε (1) where DBT is the double bark thickness at any location along the stem (mm), dob is the diameter outside bark at that location (mm), ε is the residual of the model, β0, and β1, are the regression coefficients. Thereafter, we calculated the intraclass correlation coefficient (ICC) as the ratio of between-group variances to total residual variance applying equation (2). ICC(%)=σˆ2groupσˆ2total⁎100 (2) where ICC is the intraclass correlation coefficient, σˆ2total is the total residual variance, and σˆ2group is the variance within groups (trees or plots in our case). Equation (1) was introduced by Loetsch et al. (1973) and used by Zacco (1974) in Sweden. Li and Weiskittel (2011, equation (7)) suggested that it is useful when applied to spruce species. For model validation, 10-fold block cross-validation was performed. In order to do this, the plots were randomly assigned to one of 10 subsamples with each of them containing the same number of plots. In each of 10 runs, all measurements of one subsample were used as validation data and the measurements in the remaining nine subsamples were used as training data. Fit statistics include mean absolute error (MAE), root mean square error (RMSE) and mean percentage bias (PB). Results from all runs were averaged. Calculating the required number of sampling locations per tree For Norway spruce, Stängle et al. (2016b) suggested taking five and three measurements per sampling location for an allowable error of 15 and 20 per cent, respectively. The bark surface pattern of Silver fir is similar to that of Norway spruce, being smooth in earlier years and becoming scaly with flaking scales with increasing tree age (Freund and Grehn, 1951). Hence, the two species may be expected to have similar variability of bark thickness per cross-section, and the same number of measurements per sampling location seems to be adequate to reach similar accuracy. The bark thickness variability along each tree bole was expressed by the variation of relative bark thickness within each tree. For each tree of dataset 1 (1970s) the number of required sampling locations per tree was calculated iteratively at predefined allowable errors of 15 and 20 per cent at a 95% confidence level using equation (3): n=(CV⋅tallowableerror(%))2 (3) where n is the required sample size, CV=SD(x)/x¯ (SD, standard deviation) and t refers to the 95 per cent-quantile of the two-tailed t distribution with n−1 degrees of freedom. Monte-Carlo simulations to estimate the required number of sample trees and plots To assess the influence of different sample sizes at the plot and regional levels on the predictive accuracy of equation (1), a Monte-Carlo simulation was performed on dataset 1 according to the method suggested in Stängle et al. (2016b): in each of 1000 runs, data were split in 35 plots for model fitting and 15 plots for model validation and equation (1) was fitted to the full training dataset and to thinned datasets representing different sample sizes at the tree and plot level. The tested sample sizes were 35, 20, 10 and 5 plots each combined with 20, 10, 5 and 1 trees. The plots and trees were drawn randomly with replacement from the training data as some plots contained fewer than 20 trees. The predictions from equations that were fit to the largest possible dataset, i.e. the full training data, were assumed to best capture the variability and were therefore considered the most optimal predictions. We used equivalence testing to determine the sample size required to have a high probability of fitting an equivalent equation as achieved with the full-data. To do so, in each of 1000 runs the most optimal predictions were compared with predictions from model fits achieved with thinned datasets using a robust two one-sided t-test (TOST) for equivalence (Robinson, 2016) with a region of equivalence of one millimetre (P < 0.05). Evaluation of different bark thickness models Six different models from the literature, which had different numbers of predictor variables at the disc and tree level, were compared for their fit quality and their predictive capacity on both datasets (1970s and 2010s). To account for the hierarchical data structure, mixed-effects modelling with random deviation of the intercept on the tree and plot level was applied, leading to equations (4–9): DBT=β0+β1dob+bi+bij+εijk (4) DBT=β0+β1dbhob+β2dob+bi+bij+εijk (5) log(DBT)=β0+β1dob+β2logdob+β3inddbh+bi+bij+εijk (6) DBT=dob(β1+β2hH+β3(hH)2+β4H)+bi+bij+εijk (7) DBT=β0+β1dob+β2dob2+β3dob3+bi+bij+εijk (8) log(DBTdob)=β0+β1(1−hH)β2+β3(hH)β4H+β5dbhob+β6Hdbhob+bi+bij+εijk (9) where inddbh is an indicator for breast height (if the measurement was taken at breast height = 1, else = 0), h is the height of the sampling location above ground, H is the total tree height (m), bi is the random effect for the ith plot, bij is the random effect for the jth tree in plot i, εijk is the residual error for the kth measurement in the jth tree in plot i, β0−6 are the regression coefficients, and other variables are defined as above. Equation (4) basically corresponds to equation (1) with the additional random effect terms. Equation (5) was introduced by Hannrup (2004) and is implemented in the harvester protocol Standard for Forest Data and Communication StanForD (Skogforsk, 2012). Equation (6) was found to be suitable for spruce by Wilhelmsson et al. (2002 equation (S5b)) and equation (7) was introduced by Cao and Pepper (1986 equation (4)) and recommended by Li and Weiskittel (2011 equation (4)). Model forms of equations (8) and (9) were described by Gordon (1983, equations (9 and 3)). The dependent variable of equations (6), (7) and (8) in the original publications was the diameter inside bark but was changed for this study so that the equations describe double bark thickness. To account for positional autocorrelation within each tree, which was not eliminated by introducing the tree-level random effect, a first-order continuous autoregressive correlation structure (CAR1) was applied (Pinheiro and Bates, 2000), which has been used widely in forestry (Weiskittel et al., 2011). This accounts for the similarity of bark thickness values within the tree along the stem, with closer positions being more similar to each other. A power variance function with dob as covariate was introduced to account for an observed larger residual spread for increasing diameters. The Bayesian information criterion (BIC) was computed to quantify an improvement of the model fit by introducing the above-described within-group correlation and heteroscedasticity structures. Fit statistics include MAE, RMSE and PB. Predictions of the models (6) and (9) were corrected after logarithmic back-transformation using the correction factor CF   =   exp(SEE2/2), where SEE is the standard error of the estimate (Baskerville, 1972; Sprugel, 1983). For the evaluation of the predictive capacity of the models (4)–(9), 10-fold block cross-validation was performed as described above for the validation of model (1). To obtain population-level predictions for the validation data, predictions were calculated without considering the nesting levels introduced in model fitting. Fit quality from the validation process was reported using the average of the MAE, the RMSE, and PB from the 10 runs. Mixed model fitting was performed using the lme function of the R package nlme (Pinheiro et al., 2016). To check if the two datasets from the 1970s and the 2010s led to significantly different parameter estimates, the model determined to best describe both datasets was fitted to each dataset separately. Prediction bias was calculated in each case and residuals from both fits were compared with equivalence testing to see whether the model form could describe both datasets with similar accuracy. Subsequently, both fitted equations were applied to calculate double bark thickness predictions for both datasets and results were compared with equivalence testing. Geospatial variation of bark thickness To interpret the variation in bark thickness between stands or regions, we included spatial and environmental predictor variables in the modelling process. As exact stand location and tree age were only known for the newly collected data, this step could only be performed for dataset 2. We selected the predictors dob, dbhob, and h/H from equations (4–9), as they showed to significantly explain bark thickness variation. The predictor inddbh was ignored as it strongly correlates to h/H. To account for climatic effects, we included annual precipitation (mm) and mean annual air temperature (°C). Elevation of the sample plots (metres above sea level) and geographic coordinates (Gauss–Krueger projection, in m) were included to evaluate spatial effects that could be explained by climatic as well as by biotic or abiotic factors such as site index and soil materials, respectively. Additionally, we included tree age to account for the growth rate, which can be influenced by site quality and tree-individual factors, such as suppression by larger trees in a tree’s early growth phase. Generalized additive mixed models (GAMMs) were chosen for this analysis and predictor variables were introduced as thin-plate-regression-splines using the gamm function of the R package mgcv (Wood, 2011). Geographic coordinates were introduced as tensor product splines to model a two-dimensional surface as suggested in Wood (2006). Thus, the applied equation (10) is similar to equation (4) plus relative height, age, the geographic surface and climatic variables: DBT=β+fs(dob)+fs(dbhob)+fs(hH)+fs(Age)+fte(XCoord,YCoord)+fs(Elev.)+fs(Prec.)+fs(Temp.)+bi+bij+εijk (10) where fs stands for spline smooth functions, Elev. is the elevation, Prec. is the annual precipitation, Temp. is the mean annual temperature, fte(XCoord, YCoord) is a tensor product spline of geographic coordinates, bi is the random effect for the ith plot, bij is the random effect for the jth tree in plot i, εijk is the residual error for the kth measurement in the jth tree in plot i, and other variables are defined as above. Fit statistics of sub-models of equation (10) were compared after dropping single explanatory variables to test for the significance of each model parameter using the dredge function of the R package MuMIn (Barton, 2016). BIC was used to rank the fits achieved by this procedure. Bark thickness equations for roundwood volume determination To test the effect of different log lengths on the accuracy of bark thickness equations, taper curves of trees from dataset 2 were used to create a set of virtual logs of different length according to the above-described bucking rules (long logs and 5-m segments). From those logs, rounded mid-diameters and small end diameter were extracted as sets of reference diameters, on which bark thickness equations were fit. Additionally, the diameter and bark measurements from all sampling locations (every 2 m in each tree) were used as a control set of reference diameters. Equations (1) and (11) were fit on the different sets of reference diameters and the better fit for each set was selected using BIC. Both models only have dob as the predictor variable; however equation (11) also includes it as a quadratic term: DBT=β0+β1dob+β2dob2+ε (11) Depending on the species, one of equations (1) and (11) is currently used to estimate bark thickness at log midpoints of manually measured logs in Germany. For Silver fir, the currently used equation corresponds to equation (11) with a negative coefficient for the quadratic term leading to a flattening of the equation curve at larger diameters (Altherr et al., 1978). After evaluating the effect of log length, the effect of different fitting datasets on the developed bark thickness equations was assessed. Therefore, equations (1) and (11) were fitted on mid diameters of virtual long logs from each dataset, and again, the better fit was chosen using BIC. From the developed bark thickness equations, integer bark deduction values of 1–4 cm were assigned to corresponding diameter classes. At the class midpoints 15, 25 and 35 mm double bark thickness, the class boundaries were set for deduction values of 2, 3 and 4 cm. The economic implication of updated bark deduction values was evaluated by comparing the already established bark deduction values, which correspond to the calculated values from dataset 1 (Anonymous, 2015), and the newly calculated values from dataset 2. In order to do so, net log volume after bark subtraction was calculated using outside-bark diameter measurements of the logs from dataset 3. Results Sampling design and sample size With increasing total tree height, more sampling locations within each tree were required to describe the variance of relative bark thickness (Figure 2). The calculated average number of sampling locations that was required for an allowable error of 15 per cent was 6.5 ± 3.3 (mean ± SD). A 20 per cent allowable error resulted in an average of 4.2 ± 2.3 measurements. Figure 2 View largeDownload slide Density distribution of required number of measurement locations per tree grouped by tree height classes. The number of trees per class is shown at the top. The mean per group is indicated by bold horizontal lines. Figure 2 View largeDownload slide Density distribution of required number of measurement locations per tree grouped by tree height classes. The number of trees per class is shown at the top. The mean per group is indicated by bold horizontal lines. Fit statistics of the linear model (1) revealed larger errors and higher variances for the smaller dataset 2 (Table 2). A strong autocorrelation of bark thickness within plots and within trees could be shown with the between-tree ( σˆtree2) and within-tree ( σˆε2) variation each making up about one-third of total variation ( σˆtotal2). The fewer trees were used for model fitting of equation (1), the less runs resulted in equation fits leading to equivalent predictions as the equation fit on the full data (Figure 3). In each sample size simulation with the same number of plots, the number of trees hardly affected the precision of the equation coefficients as long as five or more trees were sampled per plot. Predictions were equivalent to predictions of the full model in more than 94 per cent of the iterations, as long as five trees were sampled in each of 35 plots. Table 2 Parameter estimates with their standard error in parentheses, summary of fit statistics, and correlation structure of both datasets using equation (1) Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Note: MAE is the mean absolute error, RMSE is the root mean square error, PB is the percentage bias, σˆ2total is the total residual variance, σˆ2plot is the variance between plots, σˆ2tree is the variance between trees, and σˆ2ε is the variance within trees. Table 2 Parameter estimates with their standard error in parentheses, summary of fit statistics, and correlation structure of both datasets using equation (1) Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Note: MAE is the mean absolute error, RMSE is the root mean square error, PB is the percentage bias, σˆ2total is the total residual variance, σˆ2plot is the variance between plots, σˆ2tree is the variance between trees, and σˆ2ε is the variance within trees. Figure 3 View largeDownload slide Results from the equivalence test of equation (1) for the 16 different sample size combinations. Those shown are proportion of runs (out of 1000) that resulted in equivalent (P < 0.05, region of equivalence = 1 mm) predictions from the full training dataset (larger values are better). Labels on the x-axis refer to the number of plots and the number of trees per plot, respectively. Figure 3 View largeDownload slide Results from the equivalence test of equation (1) for the 16 different sample size combinations. Those shown are proportion of runs (out of 1000) that resulted in equivalent (P < 0.05, region of equivalence = 1 mm) predictions from the full training dataset (larger values are better). Labels on the x-axis refer to the number of plots and the number of trees per plot, respectively. Model evaluation Testing the different model forms from the literature showed that equation (6) could best describe both datasets in terms of MAE, RMSE and PB (Table 3). Blocked cross-validation revealed that equations (5) and (9) were best for predictions for datasets 1 and 2, respectively. Large errors in the validation process showed that equation (7) was not useful for predicting double bark thickness of Silver fir. Parameter estimates of all tested models can be found in Table 4. The prediction bias (expressed by the residuals) of equation (5) was equivalent for the two datasets (P < 0.05, region of equivalence: 0.5 mm), showing that the model could describe both datasets equally well. However, bark thickness was larger in dataset 1. Predictions for the newly assessed data (dataset 2) were significantly higher (P < 0.05, region of equivalence: 0.5 mm) when the model had been fit on dataset 1 compared with predictions made with a model fit from dataset 2. Table 3 Summary of fit statistics for mixed-effects models (4)–(9) and the respective prediction bias using only the fixed terms Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Best values per dataset are printed in bold. Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; MAE = mean absolute error; RMSE = root mean square error; PB, percentage bias. Table 3 Summary of fit statistics for mixed-effects models (4)–(9) and the respective prediction bias using only the fixed terms Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Best values per dataset are printed in bold. Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; MAE = mean absolute error; RMSE = root mean square error; PB, percentage bias. Table 4 Parameter estimates with their standard error in parentheses of equations (4–9) for both datasets Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; δ = random effect variance; CAR1 = first-order continuous autoregressive term; σ = residual standard error. Table 4 Parameter estimates with their standard error in parentheses of equations (4–9) for both datasets Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; δ = random effect variance; CAR1 = first-order continuous autoregressive term; σ = residual standard error. Geospatial effects The GAMM fitting of equation (10) and its sub-models obtained the best fit (BIC = 12,109.4) when only the tree-inherent variables dob, dbhob, h/H and Age were included as predictors, which is why only these predictors were kept for further analysis. Removing h/H as well only slightly reduced fit quality (BIC = 12,119.1). Total tree height, geographic coordinates, elevation and the climatic factors of temperature and precipitation could not significantly contribute to explaining the variation of bark thickness. The estimated partial effects of the remaining predictors on double bark thickness are shown with their smooth curves in Figure 4. Double bark thickness showed an increasing trend with increasing dob (Figure 4a) and with dbhob (Figure 4b). The effect of h/H was small, but obvious at the lowest measurement point at breast height (Figure 4c). A significant increase of DBT with increasing age was observed (Figure 4d). Figure 4 View largeDownload slide Estimated smooth functions describing the partial effect of the four covariates on double bark thickness (DBT): (a) the smooth of diameter outside bark, (b) the smooth of dbh outside bark, (c) the smooth of relative tree height, and (d) the smooth of age. Shaded areas are 95%-confidence intervals. Figure 4 View largeDownload slide Estimated smooth functions describing the partial effect of the four covariates on double bark thickness (DBT): (a) the smooth of diameter outside bark, (b) the smooth of dbh outside bark, (c) the smooth of relative tree height, and (d) the smooth of age. Shaded areas are 95%-confidence intervals. Bark thickness equations for roundwood volume determination Relative bark thickness was quite constant in the lower three quarters of tree height and increased slightly further up in the tree. Depending on the bucking pattern with varying log length, midpoints of logs are positioned at different relative tree heights. Nevertheless, results show that bark thickness equations referring to mid-diameters of long logs and of 5-m segments differ only slightly from each other (Table 5, Figure 5). When the equation coefficients for long logs were used as opposed to coefficients for 5-m segments, the calculated log volume was 0.8 and 0.7 per cent higher for long logs (n = 192) and 5-m segments (n = 807), respectively. The volumes of the individual logs were equivalent (paired TOST, P < 0.05, region of equivalence 0.01 m3) using the two sets of coefficients for the 5-m segments, but not equivalent for the long logs. Predicted bark thickness for small ends of 5-m segments with a diameter of up to 40 cm was higher than for the other tested reference diameters (Figure 5). Table 5 Best fit model coefficients with their standard error in parentheses and coefficient of determination of the functions (1) or (11) predicting double bark thickness with different reference diameters for datasets 1 and 2. All sampling locations means every 2 m in each tree. Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Table 5 Best fit model coefficients with their standard error in parentheses and coefficient of determination of the functions (1) or (11) predicting double bark thickness with different reference diameters for datasets 1 and 2. All sampling locations means every 2 m in each tree. Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Figure 5 View largeDownload slide Regression lines describing double bark thickness against the rounded diameter using different reference diameters, which were gained from different bucking variants of virtual logs produced from dataset 2. For comparison, also all sampled locations (every 2 m along the tree bole) are used as reference diameters. Figure 5 View largeDownload slide Regression lines describing double bark thickness against the rounded diameter using different reference diameters, which were gained from different bucking variants of virtual logs produced from dataset 2. For comparison, also all sampled locations (every 2 m along the tree bole) are used as reference diameters. Bark thickness at midpoints of long logs was smaller for virtual logs of dataset 2 than those of dataset 1 (Table 5, Figure 6). In Germany, the currently suggested boundaries for the diameter classes referring to bark deduction values of 2, 3 and 4 cm have been developed using dataset 1 and are 23, 39 and 56 cm, respectively (Anonymous, 2015). The developed bark thickness equation from dataset 2 would lead to new class boundaries at 25, 43 and 61 cm. Assuming that harvesting intensity and assortments are similar to those of the last 5 years, the application of the new class boundaries would mean a 0.8 per cent annual increase in sales volume of manually measured Silver fir logs. Figure 6 View largeDownload slide Double bark thickness plotted against the rounded mid-diameter with 95% confidence intervals of the regression lines. Vertical lines indicate the class boundaries of integer bark deduction values (horizontal lines: class midpoints 10–20, 30–40 and 40–50 mm dbt). Shifts are resulting from more shallow-barked Silver firs in the more recent data (1970s vs 2010s). Figure 6 View largeDownload slide Double bark thickness plotted against the rounded mid-diameter with 95% confidence intervals of the regression lines. Vertical lines indicate the class boundaries of integer bark deduction values (horizontal lines: class midpoints 10–20, 30–40 and 40–50 mm dbt). Shifts are resulting from more shallow-barked Silver firs in the more recent data (1970s vs 2010s). Discussion The development of accurate bark thickness equations depends on how well within-tree and between-tree variability is captured by an appropriate sampling design, adequate sample sizes, and the predictor variables in the chosen model. This study showed that the sampling design with 2-m increments between the sampling locations starting from 1.3 m above ground was adequate to assess the variability within a stem in most cases. Thus, this design is recommended for further studies of species with similar bark morphology. The sampling design with regularly spaced sampling locations (every 2 m along the tree bole) will lead to more measurements with increasing tree height. In contrast, a fixed number of relative positions along the tree bole, as is suggested in other studies (e.g. Korell, 1972; Feduccia and Mann, 1976; Kozak and Yang, 1981; Gordon, 1983; Laasasenaho et al., 2005), could reduce the support for each data point and would require additional effort in the field to recompute relative positions for each single tree. As shown by Stängle et al. (2016b), the type of bark gauge used in this study overestimated bark thickness by 0.52 mm (SD = 1.59 mm) when compared with CT-derived measurements on Norway spruce logs. This could potentially lead to biased bark thickness predictions for Silver fir as well. If a more accurate measurement device was used, a smaller number of measurements per location would be sufficient to reach the same level of accuracy. The observed high intraclass correlation within the residual variance of equation (1) shows that the relationship between diameter and bark thickness varies strongly between individual trees and plots. This was also supported by the Monte-Carlo simulations representing different sample sizes. The number of sampled plots strongly determined the number of iterations with predicted bark thickness values that were equivalent to predictions from the full-data model (Figure 3). Model choice is an important factor influencing the response in predictive capacity of bark thickness equations on sample size. For Norway spruce it was shown that a more complex model with a higher fit quality required fewer sampled plots and trees per plot (Stängle et al., 2016b). The chosen method of drawing plots and trees with replacement from the training data was necessary to account for differing numbers of trees per sampled plot and differing diameter distributions between plots. The variation observed cannot, therefore, be attributed only to sample size but also to the data-derived uncertainty. For that reason, the suggested sample sizes might slightly overestimate those actually required to achieve the intended precision of predictions. Results showed that a rather complex nonlinear model (equation (9)) and quite a simple model (equation (5)) were most suited for predicting bark thickness of Silver fir. Model (6) described both datasets best, but it failed to predict bark thickness as precisely as models (5) and (9). The same results were reported for Norway spruce in the same region (Stängle et al., 2017). Depending on the purpose of application, different explanatory variables are available and different degrees of computational complexity are preferable. If bark thickness is to be modelled for inventory data, for example, a complex model with many explanatory variables, such as model (9) is suggested because taper curves outside bark are often available. If, in contrast, total and relative tree height are not known, as in the situation when a harvester measures a tree, the best possible model without these covariates should be selected (equation (5)). Effects of total tree height were generally absent and relative tree height only marginally explained variation in bark thickness. Equation (7), which is mainly based on relative and total tree height as predictors, cannot be recommended for Silver fir, although it was proposed for different Pinus species (Cao and Pepper, 1986). The high influence of absolute and relative tree height on bark thickness variation for pines was also confirmed by Murphy and Cown (2015), who reported bark volume variation within trees, between trees, and between sites for Pinus radiata (D. Don). For Norway spruce, effects of tree height and relative tree height on bark thickness variation were shown as well (Stängle et al., 2017). As Norway spruce and most Pinus species are less shade tolerant than Silver fir, a higher correlation between total tree height and age can be expected for those species. Therefore, the effect of total tree height for those species could, at least partially, be caused by collinearity. The variation observed between the sample plots could not be linked to large-scale spatial patterns, such as climatic zones and effects of geographic coordinates were generally absent. Moreover, the climatic factors of temperature and precipitation as well as elevation could not explain differences between plots. Tree age, however, was found to have a significant positive effect on bark thickness. These findings are in accordance to those of Laasasenaho et al. (2005), Sonmez et al. (2007) and Stängle et al. (2017) for Norway spruce. This leads to the conclusion that higher growth rates result in smaller relative bark thickness for those species. Other studies also show that increased site quality, which can usually be linked to faster growth of individual trees, reduced bark thickness of Silver fir (Božić et al., 2007). It seems that no single environmental predictor can explain differences in bark thickness. Rather, the combination of many factors influencing the productivity of a stand and the growth of individual trees is required. For future studies of bark thickness, the assessment of tree age or site index is, therefore, strongly recommended. The observation that faster growing trees have thinner bark might explain the observed thinner bark in the more recent measurements of dataset 2. Since growth of Silver fir in southwestern Germany has slightly accelerated in the twentieth century (Kohnle et al., 2014), it can be expected that the trees of the new dataset have grown faster than those trees measured 40 years earlier. The small difference in relative bark thickness along the stem was apparent when bark thickness equations were fitted on mid-diameters of long logs and of 5-m segments. As there were hardly any differences in volume calculation, the same equation coefficients could be used for different log assortments. However, the small ends of 5-m segments showed larger relative bark thickness as they are positioned further up in the tree. This difference, especially for small diameters, suggests that different equation coefficients should be used for bucking optimization that is based on minimum small end diameters. Alternatively, a model, which includes information on the relative height of the sampling location, could be chosen. Thinner bark in the more recent measurements resulted in different diameter classes for integer bark deduction for the two datasets. The values that are currently in use in Germany (Anonymous, 2015) are based on the older measurements. Therefore, they tend to overestimate bark thickness and lead to biased inside-bark volume estimations. Conclusion To parameterize bark thickness equations, within-tree and between-tree variation in Silver fir bark thickness requires several sampling locations per tree and a minimum of five sampled trees in at least 35 plots in the study region. Bark thickness of Silver fir can be described well using tree diameter at the sampling location and the breast height diameter as predictor variables. Faster growing trees were shown to have relatively thinner bark. From the demonstrated influence of tree age – and therefore growth rate – on bark thickness, it can be inferred that changing growth conditions can alter bark allometry, which is why bark thickness equations should be regularly reviewed for their validity. For measuring roundwood, the same bark deduction values can be suggested for long logs (6–21 m) and 5-m segments. For bucking optimization based on minimum small end diameters, different equation coefficients that reflect higher relative bark thickness at larger relative tree height are recommended. The findings show the need to update the current bark deduction values for manual log measurement in the state of Baden-Württemberg in order to estimate log volume more accurately. As the same values are currently used in large parts of Central Europe, further validation measurements in other federal states in Germany and in neighbouring countries are suggested. Acknowledgements The methods described in this manuscript were partially developed with support from Professor Aaron Weiskittel (University of Maine) and additional statistical advice from Dr Gerald Kändler (FVA Baden-Württemberg) was appreciated. We would like to thank two anonymous reviewers, the Editor-in-Chief Dr Gary Kerr, and especially the handling Editor Professor Alexis Achim for their comments and suggestions that helped to improve previous versions of this manuscript. Conflict of interest statement None declared. Funding This study was funded by the Ministry for Rural Affairs and Consumer Protection Baden-Württemberg, Germany (MLR) and was conducted at the Forest Research Institute Baden Württemberg (FVA). References Altherr , E. , Unfried , P. , Hradetzky , J. and Hradetzky , V. 1974 Statistische Rindenbeziehungen als Hilfsmittel zur Ausformung und Aufmessung unentrindeten Stammholzes: Teil I: Kiefer, Buche, Hainbuche, Esche und Roterle. In Mitteilungen der Forstlichen Versuchs- und Forschungsanstalt, Forstliche Versuchs- und Forschungsanstalt Baden-Württemberg, Freiburg im Breisgau, Germany, pp. 137. Altherr , E. , Unfried , P. , Hradetzky , J. and Hradetzky , V. 1978 Statistische Rindenbeziehungen als Hilfsmittel zur Ausformung und Aufmessung unentrindeten Stammholzes: Teil IV: Fichte, Tanne, Douglasie und Sitka-Fichte. In Mitteilungen der Forstlichen Versuchs- und Forschungsanstalt, Forstliche Versuchs- und Forschungsanstalt Baden-Württemberg, Freiburg im Breisgau, Germany, pp. 294. Anonymous . 2015 Rahmenvereinbarung für den Rohholzhandel in Deutschland (RVR). Deutscher Forstwirtschaftsrat e.V.; Deutscher Holzwirtschaftsrat e.V. http://www.saegeindustrie.de/rvr/docs/dynamisch/6205/rvr_gesamtdokument_2.auflage_stand_oktober_2015.pdf (accessed on 14 July, 2016) Antony , F. , Schimleck , L.R. , Daniels , R.F. , Clark , A. , Borders , B.E. , Kane , M.B. , et al. 2015 Whole-tree bark and wood properties of loblolly pine from intensively managed plantations . For. Sci. 61 , 55 – 66 . Barton , K. 2016 MuMIn: Multi-model inference: R package version 1.15.6. Baskerville , G.L. 1972 Use of logarithmic regression in the estimation of plant biomass . Can. J. For. Res. 2 , 49 – 53 . Google Scholar CrossRef Search ADS Božić , M. , Čavlović , J. , Vedriš , M. and Jazbec , M. 2007 Modeling bark thickness of silver fir trees (Abies alba Mill.) . Šumarski list 131 , 3 – 12 . Cao , Q.V. and Pepper , W.D. 1986 Predicting inside bark diameter for Shortleaf, Loblolly, and Longleaf pines . South. J. Appl. For. 10 , 220 – 224 . Cellini , J.M. , Galarza , M. , Burns , S.L. , Martinez-Pastur , G.J. and Lencinas , M.V. 2012 Equations of bark thickness and volume profiles at different heights with easy-measurement variables . Forest Sys. 21 , 23 – 30 . Google Scholar CrossRef Search ADS Dimitrov , E.T. 1976 Mathematical models for determining the bark volume of spruce in relation to certain mensurational characteristics (translated from Hungarian) . Gorstoskopanska Nauka (Forest Science) 13 , 52 – 63 . Doruska , P.F. , Patterson , D. , Hartley , J. , Hurd , M. and Hart , T. 2009 Newer technologies and bioenergy bring focus back to bark factor equations . J. For. 107 , 38 – 43 . EUFORGEN . 2009 Distribution map of Silver fir (Abies alba). http://www.euforgen.org (accessed on 14 July, 2016). Feduccia , D.P. and Mann , W.F. 1976 Bark thickness of 17-year-old loblolly pine planted at different spacings. USDA For. Serv., p. 2. Fonseca , M.A. 2005 The Measurement of Roundwood: Methodologies and Conversion Ratios . CABI , p. 269 . Google Scholar CrossRef Search ADS Freund , H. and Grehn , J. 1951 Handbuch der Mikroskopie in der Technik. 1 - Bd. 5, Mikroskopie des Rohholzes und der Rinden . Umschau-Verlag , p. 456 . Gordon , A. 1983 Estimating bark thickness of Pinus radiata . New Zeal. J. For. Sci 13 , 340 – 348 . Hannrup , B. 2004 Funktioner för skattning av barkens tjocklek hos tall och gran vid avverkning med skördare. (Functions for prediction of bark thickness of Norway spruce and Scots pine at CTL-harvesting). Skogforsk, Uppsala, p. 33. Hoffmann , J. 1958 Untersuchungen über die Rindenstärke der Fichte auf verschiedenen Standorten im südöstlichen Thüringer Wald . Wiss Z. TU Dresden 7 , 361 – 368 . Kohnle , U. , Albrecht , A. , Lenk , E. , Ohnemus , K. and Yue , C. 2014 Growth trends driven by environmental factors extracted from long term experimental data in southwest Germany . AFJZ 185 , 97 – 117 . Kohnle , U. , Hein , S. , Sorensen , F.C. and Weiskittel , A.R. 2012 Effects of seed source origin on bark thickness of Douglas-fir (Pseudotsuga menziesii) growing in southwestern Germany . Can. J. For. Res. 42 , 382 – 399 . Google Scholar CrossRef Search ADS Korell , U. 1972 Über Rindendicken der Fichte . Beitr. Forstwirtsch 2 , 54 – 55 . Kozak , A. and Yang , R.C. 1981 Equations for estimating bark volume and thickness of commercial trees in British Columbia . For. Chron. 57 , 112 – 115 . Google Scholar CrossRef Search ADS Kublin , E. and Breidenbach , J. 2015 TapeR: Flexible tree taper curves based on Semiparametric Mixed Models: R package version 0.3.3. Laasasenaho , J. , Melkas , T. and Aldén , S. 2005 Modelling bark thickness of Picea abies with taper curves . For. Ecol. Manage. 206 , 35 – 47 . Google Scholar CrossRef Search ADS Li , R. and Weiskittel , A.R. 2011 Estimating and predicting bark thickness for seven conifer species in the Acadian Region of North America using a mixed-effects modeling approach: comparison of model forms and subsampling strategies . Eur. J. Forest. Res. 130 , 219 – 233 . Google Scholar CrossRef Search ADS Loetsch , F. , Zöhrer , F. and Haller , K.E. 1973 Forest Inventory Vol. II . BLV Verl. Ges , p. 436 . Marshall , H.D. , Murphy , G.E. and Lachenbruch , B. 2006 Effects of bark thickness estimates on opimal log merchandising . For. Prod. J 56 , 87 – 92 . Martin , R.E. and Crist , J.B. 1970 Elements of bark structure and terminology . Wood Fiber 2 , 269 – 279 . Murphy , G. and Cown , D. 2015 Within-tree, between-tree, and geospatial variation in estimated Pinus radiata bark volume and weight in New Zealand . New Zeal. J. For. Sci. 45 , 18 . Google Scholar CrossRef Search ADS Niklas , K.J. 1999 The mechanical role of bark . Am. J. Bot. 86 , 465 – 469 . Google Scholar CrossRef Search ADS PubMed Persson , B. and Downie , B. 1992 Variation in bark thickness of young Pinus contorta var. latifolia Engelm. in Sweden . Scand. J. For. Res. 7 , 99 – 106 . Google Scholar CrossRef Search ADS Pinheiro , J. , Bates , D. , DebRoy , S. , Sarkar , D. and Team , R.C. 2016 nlme: Linear and Nonlinear Mixed Effects Models: R package version 3.1–124. Pinheiro , J.C. and Bates , D.M. 2000 Mixed-effects Models in S and S-PLUS . Springer , p. 528 . Google Scholar CrossRef Search ADS Pizzi , A. 2016 Wood products and green chemistry . Ann. For. Sci. 73 , 185 – 203 . Google Scholar CrossRef Search ADS R Core Team . 2015 R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing . Robinson , A.P. 2016 equivalence: Provides tests and graphics for assessing tests of equivalence: R package version 0.7.1. Schielzeth , H. and Nakagawa , S. 2013 Nested by design: model fitting and interpretation in a mixed model era . Methods Ecol. Evol. 4 , 14 – 24 . Google Scholar CrossRef Search ADS Schmidt-Vogt , H. 1986 Die Fichte–Ein Handbuch in zwei Bänden, Band II/1: Wachstum, Züchtung, Boden, Umwelt, Holz . Parey , p. 563 . Skogforsk . 2012 Standard for Forest Data and Communication Appendix: Definitions of variables - General and country specific, http://www.skogforsk.se/contentassets/b063db555a664ff8b515ce121f4a42d1/appendix1_eng_120418.pdf (accessed on 27 June, 2016). Sonmez , T. , Keles , S. and Tilki , F. 2007 Effect of aspect, tree age and tree diameter on bark thickness of Picea orientalis . Scand. J. For. Res. 22 , 193 – 197 . Google Scholar CrossRef Search ADS Sprugel , D. 1983 Correcting for bias in log-transformed allometric equations . Ecology 64 , 209 – 210 . Google Scholar CrossRef Search ADS Stängle , S.M. , Sauter , U.H. , Brüchert , F. and Kändler , G. 2016 a Überprüfung der Rindenabzugswerte für Fichten-Stammholz in Baden-Württemberg (A review of bark deduction values for Norway spruce logs in Baden-Württemberg) . Forstarchiv 87 , 162 – 169 . Stängle , S.M. , Sauter , U.H. and Dormann , C.F. 2017 Comparison of models for estimating bark thickness of Picea abies in southwest Germany: the role of tree, stand, and environmental factors . Ann. For. Sci. 74 , 16 . Google Scholar CrossRef Search ADS Stängle , S.M. , Weiskittel , A.R. , Dormann , C.F. and Brüchert , F. 2016 b Measurement and prediction of bark thickness in Picea abies: assessment of accuracy, precision, and sample size requirements . Can. J. For. Res. 46 , 39 – 47 . Google Scholar CrossRef Search ADS Temesgen , H. , Affleck , D. , Poudel , K. , Gray , A. and Sessions , J. 2015 A review of the challenges and opportunities in estimating above ground forest biomass using tree-level models . Scand. J. For. Res. 30 , 326 – 335 . Weiskittel , A.R. , Hann , D.W. , Kershaw , J.A. and Vanclay , J.K. 2011 Forest Growth and Yield Modeling , 415 . Wiley . Google Scholar CrossRef Search ADS Wilhelmsson , L. , Arlinger , J. , Spångberg , K. , Lundqvist , S.-O. , Grahn , T. , Hedenberg , Ö. and Olsson , L. 2002 Models for Predicting Wood Properties in Stems of Picea abies and Pinus sylvestris in Sweden . Scand. J. For. Res. 17 , 330 – 350 . Google Scholar CrossRef Search ADS Wood , S.N. 2006 Generalized Additive Models: An Introduction With R . Chapman & Hall/CRC , p. 392 . Wood , S.N. 2011 Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models . J. Roy. Stat. Soc. Ser. B (Stat. Method.) 73 , 3 – 36 . Google Scholar CrossRef Search ADS Zacco , P. 1974 Barktjockleken hos sågtimmer (The bark thickness of saw logs). In Rapport . The Swedish University of Agricultural Sciences, Department of Forest Products. Stockholm , p. 53 . © Institute of Chartered Foresters, 2017. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Forestry: An International Journal Of Forest Research Oxford University Press

Modelling the variation of bark thickness within and between European silver fir (Abies alba Mill.) trees in southwest Germany

Loading next page...
 
/lp/ou_press/modelling-the-variation-of-bark-thickness-within-and-between-european-GPluvAoLRC
Publisher
Oxford University Press
Copyright
© Institute of Chartered Foresters, 2017. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com.
ISSN
0015-752X
eISSN
1464-3626
D.O.I.
10.1093/forestry/cpx047
Publisher site
See Article on Publisher Site

Abstract

Abstract This study examined bark thickness variability of Silver fir (Abies alba Mill.) in southwest Germany over time and space by comparing a dataset of bark measurements from the 1970s with more recent assessments. Within-tree variability of bark thickness was analysed to estimate the required number of sampling locations per tree for defined accuracy levels. A range of models from the literature that predict bark thickness were compared for their predictive performance and Monte-Carlo simulations were used to estimate the effect of the number of sample trees and plots on the precision of the predictions. In addition, net log volume after bark subtraction was compared for logs of varying lengths to assess the influence of log assortments on calculated sales volume. Results show that several sampling locations are needed per tree and that at least five trees from at least 35 plots should be selected for measurements in the study region. For practical applications, diameter outside bark and breast height diameter are suggested as explanatory variables for models that predict double bark thickness. Additionally, relative tree height and age – and therefore growth rate – significantly improved predictions; however, environmental factors could not explain the variation between stands. Log lengths from 5 to 21 m only slightly influenced bark thickness equations that were fit on measurements at log midpoints. The findings highlight the need to consider bark thickness variability at different levels when developing bark thickness equations. In general, bark thickness was found to be smaller in more recent assessments and this indicates the need to regularly review existing bark equations for their validity. Introduction The diameter inside bark is a crucial measure in forestry that is used to calculate wood volume of logs and trees and to optimize bucking for defined inside-bark diameters. However, tree and log diameters are usually measured outside bark and bark thickness has to be estimated. Incorrect bark thickness or bark volume estimates could lead to inaccurate wood volume estimates in forest inventories, in increment studies or in the log trade (Marshall et al., 2006). Furthermore, the interest in accurate bark thickness estimates is driven by the need to obtain not only accurate inside-bark diameters and volumes, but also accurate estimates of bark volume. The importance of accurate bark thickness estimates has increased with the shift in the commercial relevance of bark from an unwanted residue to a valuable fuel and feedstock for high-value biomaterial (Doruska et al., 2009), such as bark tannin based foams (e.g. Pizzi, 2016). The estimation of available bark biomass is important to assess the potential of such technologies for generating possible additional income for the forestry sector. Additionally, bark volume estimates are needed for biomass estimates, which are becoming more important for quantifying carbon stocks (Temesgen et al., 2015). In this study, we use a common definition of bark that includes all tissues outside the vascular cambium, and comprises the secondary phloem up to the last formed periderm and the rhytidome, which comprises all layers of dead tissue outside the currently active periderm (Martin and Crist, 1970). For many coniferous species, it has been shown that bark thickness can be described well by diameter outside bark at the sampling location and additional variables, as it often correlates with diameter at breast height, total tree height and height of the sampling location within the tree (e.g. Li and Weiskittel, 2011). Several external factors that correlate with the relative bark thickness of trees, i.e. the proportion of the outside-bark diameter that constitutes bark, have been reported. Most of these can be interpreted as a consequence of slower tree growth, which leads to a larger relative bark thickness. A correlation between a low site index or yield class of stands with a larger relative bark thickness was shown for Norway spruce (Picea abies (L.) Karst) (Hoffmann, 1958; Dimitrov, 1976; Schmidt-Vogt, 1986) and Silver fir (Abies alba Mill.) (Božić et al., 2007). For Norway spruce (Laasasenaho et al., 2005, Stängle et al., 2017) and Silver fir (Božić et al., 2007), tree age was shown to have an influence on the diameter-to-bark ratio. Sonmez et al. (2007) also reported a positive effect of tree age on bark thickness for Picea orientalis (L. Link). Contrasting results were reported for the broadleaved Nothofagus pumilio (Poepp. & Endl.) Krasser by Cellini et al. (2012), where relative bark thickness was lower on low productivity sites. For Pinus sylvesteris (L.), the latitude of forest stands – which mainly determines differences in temperature and length of the growing season – could help to explain variation in bark thickness (Wilhelmsson et al., 2002). For Pinus taeda (L.) in the United States, there were no observed effects of silvicultural treatments, such as initial planting density and thinning, on whole-tree bark percentage. Only very intense treatments including fertilization and early-age weed control affected bark volume (Antony et al., 2015). Besides growth conditions, genetics can also determine bark thickness. Provenance has been shown to influence the phenotypic development of bark thickness of Pinus contorta var. latifolia Engelm. (Persson and Downie, 1992) and Pseudotsuga menziesii (Mirb.) Franco (Kohnle et al., 2012). To develop accurate equations for predicting regional bark thickness, the intraspecific within-tree and between-tree variability of bark thickness have to be considered. Bark thickness variation within trees occurs at two levels. First, variation at the disc level around the circumference of the tree can be caused by a non-uniform bark thickness that is correlated with elliptical stem growth (e.g. Niklas, 1999) or a bark surface pattern with fissures, cracks or scales. The second level of within-tree variation occurs between sampling locations along the stem axis. To capture both sources of within-tree variability, several bark thickness measurements are usually taken at each of several sampling locations along the stem. For Norway spruce, variability within and between trees and between sites indicates the importance of multiple sampling points within each tree of a range of selected trees covering a range of sampling sites throughout the study region (Stängle et al., 2016b). Roundwood volume calculation in Germany is performed according to the Huber formula, which is based on a cylinder volume calculation using log length and the rounded diameter inside bark at log midpoint, which is located at exactly halfway along the log’s length (Fonseca, 2005). In the study region, long logs (6–21 m) are mostly measured manually, whereas shorter logs from cut-to-length operations are usually measured at the mill site using electronic scanning devices. Measurement instructions for manual log measurement in Germany require the operator to take two diameter measurements at log midpoint (approximately perpendicularly) and round them to the lower full centimetre. The mean of these rounded values is also rounded down and species-specific fixed integer bark deduction values are subtracted to estimate the inside-bark diameter (Anonymous, 2015). Geospatial patterns of bark thickness within a study region could lead to skewed wood volume estimations if only one bark thickness equation and one set of fixed deduction values is used for the whole region. Additionally, relative bark thickness can vary along the stem and if the applied bark thickness function does not consider the position in the stem, volume calculations could be biased for different assortments. If long logs of up to 21 m length are produced, for example, most log midpoints would be positioned in the lower half of the stem, whereas midpoints of 5-m segments can basically be positioned all along the stem. The integer bark deduction values that are currently in use have been extracted from bark thickness equations that report double bark thickness as a function of mid-diameters of logs of up to 26 m in length (Altherr et al., 1978). However, different bark thickness equations could be useful for different assortments. Furthermore, bark thickness equations should be up-to-date as allometric relationships can change with changes in growth conditions. For Silver fir, growth conditions in the study area have changed in the last century with the average growth rate having increased since the 1990s (Kohnle et al., 2014). Currently, bark thickness equations applied in Germany are based on measurements from the 1970s and probably overestimate bark volume as has already been shown for Norway spruce in southwest Germany (Stängle et al., 2016a). Using Silver fir also from southwest Germany, the objectives of this study were to: estimate sample sizes required to capture bark thickness variation within trees and to examine the effect of number of sample plots and trees per plot on the precision of bark thickness equations; compare and evaluate established model forms for estimating bark thickness and test how geospatial variation could be explained by environmental factors; and estimate how log length and different datasets (measurements from the 1970s and 2010s) affect bark thickness equation coefficients for roundwood volume calculation. Materials and methods Data Bark thickness data were obtained from two sources. Dataset 1, an extensive dataset with 3008 Silver fir trees, was collected by Altherr et al. (1978) in 82 clusters. Several clusters were sampled in the same stands, resulting in 50 sampled stands. Each stand was treated as one plot in the analysis. Dataset 2 consisted of newly acquired measurements on 217 trees in 27 plots. Measurements of dataset 2 were taken in the years 2015 and 2016. Dimensional details of all trees are listed in Table 1. Trees of both datasets were sampled throughout the 35 752 km2 state of Baden-Württemberg in temperate forests with a broad range of site and stand conditions. Elevation of the plots from dataset 2 varied between 306 and 974 m above sea level (Figure 1). Mean annual air temperature varied between 5.8°C and 9.4°C and the annual precipitation was between 868 and 1944 mm. All stands were mixed-species stands with a range of ~20–80 per cent fir trees. Trees from a large range of merchantable diameter classes were selected from regularly harvested trees in selection cuts. Table 1 Characteristics of the measured trees from dataset 1 and 2. Dataset 1 originates from Altherr et al. (1978). Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Note: dbh = diameter at breast height; DBT1.3 = double bark thickness at breast height; SD = standard deviation. Table 1 Characteristics of the measured trees from dataset 1 and 2. Dataset 1 originates from Altherr et al. (1978). Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Dataset 1 (1970s) Dataset 2 (2010s) Mean SD Minimum Maximum Mean SD Minimum Maximum Age (years) 126.8 37.5 54 225 Tree height (m) 29.5 5.3 13.0 40.0 dbh (cm) 36.3 14.9 13.0 99.9 48.3 17.3 18.5 97.2 DBT1.3 (mm) 21.3 7.7 7.0 52.0 25.3 7.8 10.0 49.0 Note: dbh = diameter at breast height; DBT1.3 = double bark thickness at breast height; SD = standard deviation. Figure 1 View largeDownload slide Distribution of sample plots in the study region and number of sampled trees per plot for dataset 2. The distribution map of Silver fir originates from (EUFORGEN, 2009). Figure 1 View largeDownload slide Distribution of sample plots in the study region and number of sampled trees per plot for dataset 2. The distribution map of Silver fir originates from (EUFORGEN, 2009). The sampling for both datasets was performed as described by Altherr et al. (1974; 1978): trees were felled with chainsaws, delimbed, and measured in the forest before any further log manipulation was performed. Measurement locations were at breast height (1.3 m above ground) and along the tree bole in 2 m increments up to a top diameter of ~10–15 cm. Both diameter and bark thickness were measured twice (approximately perpendicularly) at each location using a calliper and a Swedish bark gauge, respectively. Double bark thickness was calculated as the sum of the two gauge measurements. Additional measurements for the new data were tree age, and total tree height, which was measured on the felled tree. Tree height for dataset 1 and for trees with broken tips in dataset 2 was modelled with mixed B-spline regression describing tree taper using the package TapeR (Kublin and Breidenbach, 2015) in R (R Core Team, 2015). Two predefined bucking patterns were applied to produce virtual logs of different length from both datasets. The first bucking rule aimed at long logs. From each tree one virtual log between 6 and 21 m was generated. According to the second rule, several segments of 5 m length were generated from each tree. Due to convergence problems of the taper functions, only trees with a diameter at breast height (dbh) larger than 18 cm could be used for this analysis. For each log, diameter and bark thickness at midpoint and at the small end was calculated using taper equations based on inside-bark and outside-bark diameter. To evaluate the effect of log length on equation coefficients and roundwood volume, the harvesting statistics of the State Forest Enterprise of Baden-Württemberg (ForstBW) were used. These data, hereafter referred to as dataset 3, included all manually measured logs of length longer than 6 m that were produced in the state forest of Baden-Württemberg in the period 2011–2015. Bark thickness variability at different hierarchical levels Bark thickness measurements within trees and within plots can be expected to be similar (e.g. Li and Weiskittel, 2011). The variance within trees determines the sampling design with the number of required sampling locations per tree and the number of measurements per location. If the variance between groups (such as trees and plots in our case) is large, additional information on group-level predictors could help explain this variance when modelling bark thickness (Schielzeth and Nakagawa, 2013). To quantify the degree of variation that can be accounted for by different grouping levels, we fitted the established bark thickness prediction equation (1), which does not consider any predictor variables at the tree or plot level, but only the diameter outside bark at the measurement location. DBT=β0+β1dob+ε (1) where DBT is the double bark thickness at any location along the stem (mm), dob is the diameter outside bark at that location (mm), ε is the residual of the model, β0, and β1, are the regression coefficients. Thereafter, we calculated the intraclass correlation coefficient (ICC) as the ratio of between-group variances to total residual variance applying equation (2). ICC(%)=σˆ2groupσˆ2total⁎100 (2) where ICC is the intraclass correlation coefficient, σˆ2total is the total residual variance, and σˆ2group is the variance within groups (trees or plots in our case). Equation (1) was introduced by Loetsch et al. (1973) and used by Zacco (1974) in Sweden. Li and Weiskittel (2011, equation (7)) suggested that it is useful when applied to spruce species. For model validation, 10-fold block cross-validation was performed. In order to do this, the plots were randomly assigned to one of 10 subsamples with each of them containing the same number of plots. In each of 10 runs, all measurements of one subsample were used as validation data and the measurements in the remaining nine subsamples were used as training data. Fit statistics include mean absolute error (MAE), root mean square error (RMSE) and mean percentage bias (PB). Results from all runs were averaged. Calculating the required number of sampling locations per tree For Norway spruce, Stängle et al. (2016b) suggested taking five and three measurements per sampling location for an allowable error of 15 and 20 per cent, respectively. The bark surface pattern of Silver fir is similar to that of Norway spruce, being smooth in earlier years and becoming scaly with flaking scales with increasing tree age (Freund and Grehn, 1951). Hence, the two species may be expected to have similar variability of bark thickness per cross-section, and the same number of measurements per sampling location seems to be adequate to reach similar accuracy. The bark thickness variability along each tree bole was expressed by the variation of relative bark thickness within each tree. For each tree of dataset 1 (1970s) the number of required sampling locations per tree was calculated iteratively at predefined allowable errors of 15 and 20 per cent at a 95% confidence level using equation (3): n=(CV⋅tallowableerror(%))2 (3) where n is the required sample size, CV=SD(x)/x¯ (SD, standard deviation) and t refers to the 95 per cent-quantile of the two-tailed t distribution with n−1 degrees of freedom. Monte-Carlo simulations to estimate the required number of sample trees and plots To assess the influence of different sample sizes at the plot and regional levels on the predictive accuracy of equation (1), a Monte-Carlo simulation was performed on dataset 1 according to the method suggested in Stängle et al. (2016b): in each of 1000 runs, data were split in 35 plots for model fitting and 15 plots for model validation and equation (1) was fitted to the full training dataset and to thinned datasets representing different sample sizes at the tree and plot level. The tested sample sizes were 35, 20, 10 and 5 plots each combined with 20, 10, 5 and 1 trees. The plots and trees were drawn randomly with replacement from the training data as some plots contained fewer than 20 trees. The predictions from equations that were fit to the largest possible dataset, i.e. the full training data, were assumed to best capture the variability and were therefore considered the most optimal predictions. We used equivalence testing to determine the sample size required to have a high probability of fitting an equivalent equation as achieved with the full-data. To do so, in each of 1000 runs the most optimal predictions were compared with predictions from model fits achieved with thinned datasets using a robust two one-sided t-test (TOST) for equivalence (Robinson, 2016) with a region of equivalence of one millimetre (P < 0.05). Evaluation of different bark thickness models Six different models from the literature, which had different numbers of predictor variables at the disc and tree level, were compared for their fit quality and their predictive capacity on both datasets (1970s and 2010s). To account for the hierarchical data structure, mixed-effects modelling with random deviation of the intercept on the tree and plot level was applied, leading to equations (4–9): DBT=β0+β1dob+bi+bij+εijk (4) DBT=β0+β1dbhob+β2dob+bi+bij+εijk (5) log(DBT)=β0+β1dob+β2logdob+β3inddbh+bi+bij+εijk (6) DBT=dob(β1+β2hH+β3(hH)2+β4H)+bi+bij+εijk (7) DBT=β0+β1dob+β2dob2+β3dob3+bi+bij+εijk (8) log(DBTdob)=β0+β1(1−hH)β2+β3(hH)β4H+β5dbhob+β6Hdbhob+bi+bij+εijk (9) where inddbh is an indicator for breast height (if the measurement was taken at breast height = 1, else = 0), h is the height of the sampling location above ground, H is the total tree height (m), bi is the random effect for the ith plot, bij is the random effect for the jth tree in plot i, εijk is the residual error for the kth measurement in the jth tree in plot i, β0−6 are the regression coefficients, and other variables are defined as above. Equation (4) basically corresponds to equation (1) with the additional random effect terms. Equation (5) was introduced by Hannrup (2004) and is implemented in the harvester protocol Standard for Forest Data and Communication StanForD (Skogforsk, 2012). Equation (6) was found to be suitable for spruce by Wilhelmsson et al. (2002 equation (S5b)) and equation (7) was introduced by Cao and Pepper (1986 equation (4)) and recommended by Li and Weiskittel (2011 equation (4)). Model forms of equations (8) and (9) were described by Gordon (1983, equations (9 and 3)). The dependent variable of equations (6), (7) and (8) in the original publications was the diameter inside bark but was changed for this study so that the equations describe double bark thickness. To account for positional autocorrelation within each tree, which was not eliminated by introducing the tree-level random effect, a first-order continuous autoregressive correlation structure (CAR1) was applied (Pinheiro and Bates, 2000), which has been used widely in forestry (Weiskittel et al., 2011). This accounts for the similarity of bark thickness values within the tree along the stem, with closer positions being more similar to each other. A power variance function with dob as covariate was introduced to account for an observed larger residual spread for increasing diameters. The Bayesian information criterion (BIC) was computed to quantify an improvement of the model fit by introducing the above-described within-group correlation and heteroscedasticity structures. Fit statistics include MAE, RMSE and PB. Predictions of the models (6) and (9) were corrected after logarithmic back-transformation using the correction factor CF   =   exp(SEE2/2), where SEE is the standard error of the estimate (Baskerville, 1972; Sprugel, 1983). For the evaluation of the predictive capacity of the models (4)–(9), 10-fold block cross-validation was performed as described above for the validation of model (1). To obtain population-level predictions for the validation data, predictions were calculated without considering the nesting levels introduced in model fitting. Fit quality from the validation process was reported using the average of the MAE, the RMSE, and PB from the 10 runs. Mixed model fitting was performed using the lme function of the R package nlme (Pinheiro et al., 2016). To check if the two datasets from the 1970s and the 2010s led to significantly different parameter estimates, the model determined to best describe both datasets was fitted to each dataset separately. Prediction bias was calculated in each case and residuals from both fits were compared with equivalence testing to see whether the model form could describe both datasets with similar accuracy. Subsequently, both fitted equations were applied to calculate double bark thickness predictions for both datasets and results were compared with equivalence testing. Geospatial variation of bark thickness To interpret the variation in bark thickness between stands or regions, we included spatial and environmental predictor variables in the modelling process. As exact stand location and tree age were only known for the newly collected data, this step could only be performed for dataset 2. We selected the predictors dob, dbhob, and h/H from equations (4–9), as they showed to significantly explain bark thickness variation. The predictor inddbh was ignored as it strongly correlates to h/H. To account for climatic effects, we included annual precipitation (mm) and mean annual air temperature (°C). Elevation of the sample plots (metres above sea level) and geographic coordinates (Gauss–Krueger projection, in m) were included to evaluate spatial effects that could be explained by climatic as well as by biotic or abiotic factors such as site index and soil materials, respectively. Additionally, we included tree age to account for the growth rate, which can be influenced by site quality and tree-individual factors, such as suppression by larger trees in a tree’s early growth phase. Generalized additive mixed models (GAMMs) were chosen for this analysis and predictor variables were introduced as thin-plate-regression-splines using the gamm function of the R package mgcv (Wood, 2011). Geographic coordinates were introduced as tensor product splines to model a two-dimensional surface as suggested in Wood (2006). Thus, the applied equation (10) is similar to equation (4) plus relative height, age, the geographic surface and climatic variables: DBT=β+fs(dob)+fs(dbhob)+fs(hH)+fs(Age)+fte(XCoord,YCoord)+fs(Elev.)+fs(Prec.)+fs(Temp.)+bi+bij+εijk (10) where fs stands for spline smooth functions, Elev. is the elevation, Prec. is the annual precipitation, Temp. is the mean annual temperature, fte(XCoord, YCoord) is a tensor product spline of geographic coordinates, bi is the random effect for the ith plot, bij is the random effect for the jth tree in plot i, εijk is the residual error for the kth measurement in the jth tree in plot i, and other variables are defined as above. Fit statistics of sub-models of equation (10) were compared after dropping single explanatory variables to test for the significance of each model parameter using the dredge function of the R package MuMIn (Barton, 2016). BIC was used to rank the fits achieved by this procedure. Bark thickness equations for roundwood volume determination To test the effect of different log lengths on the accuracy of bark thickness equations, taper curves of trees from dataset 2 were used to create a set of virtual logs of different length according to the above-described bucking rules (long logs and 5-m segments). From those logs, rounded mid-diameters and small end diameter were extracted as sets of reference diameters, on which bark thickness equations were fit. Additionally, the diameter and bark measurements from all sampling locations (every 2 m in each tree) were used as a control set of reference diameters. Equations (1) and (11) were fit on the different sets of reference diameters and the better fit for each set was selected using BIC. Both models only have dob as the predictor variable; however equation (11) also includes it as a quadratic term: DBT=β0+β1dob+β2dob2+ε (11) Depending on the species, one of equations (1) and (11) is currently used to estimate bark thickness at log midpoints of manually measured logs in Germany. For Silver fir, the currently used equation corresponds to equation (11) with a negative coefficient for the quadratic term leading to a flattening of the equation curve at larger diameters (Altherr et al., 1978). After evaluating the effect of log length, the effect of different fitting datasets on the developed bark thickness equations was assessed. Therefore, equations (1) and (11) were fitted on mid diameters of virtual long logs from each dataset, and again, the better fit was chosen using BIC. From the developed bark thickness equations, integer bark deduction values of 1–4 cm were assigned to corresponding diameter classes. At the class midpoints 15, 25 and 35 mm double bark thickness, the class boundaries were set for deduction values of 2, 3 and 4 cm. The economic implication of updated bark deduction values was evaluated by comparing the already established bark deduction values, which correspond to the calculated values from dataset 1 (Anonymous, 2015), and the newly calculated values from dataset 2. In order to do so, net log volume after bark subtraction was calculated using outside-bark diameter measurements of the logs from dataset 3. Results Sampling design and sample size With increasing total tree height, more sampling locations within each tree were required to describe the variance of relative bark thickness (Figure 2). The calculated average number of sampling locations that was required for an allowable error of 15 per cent was 6.5 ± 3.3 (mean ± SD). A 20 per cent allowable error resulted in an average of 4.2 ± 2.3 measurements. Figure 2 View largeDownload slide Density distribution of required number of measurement locations per tree grouped by tree height classes. The number of trees per class is shown at the top. The mean per group is indicated by bold horizontal lines. Figure 2 View largeDownload slide Density distribution of required number of measurement locations per tree grouped by tree height classes. The number of trees per class is shown at the top. The mean per group is indicated by bold horizontal lines. Fit statistics of the linear model (1) revealed larger errors and higher variances for the smaller dataset 2 (Table 2). A strong autocorrelation of bark thickness within plots and within trees could be shown with the between-tree ( σˆtree2) and within-tree ( σˆε2) variation each making up about one-third of total variation ( σˆtotal2). The fewer trees were used for model fitting of equation (1), the less runs resulted in equation fits leading to equivalent predictions as the equation fit on the full data (Figure 3). In each sample size simulation with the same number of plots, the number of trees hardly affected the precision of the equation coefficients as long as five or more trees were sampled per plot. Predictions were equivalent to predictions of the full model in more than 94 per cent of the iterations, as long as five trees were sampled in each of 35 plots. Table 2 Parameter estimates with their standard error in parentheses, summary of fit statistics, and correlation structure of both datasets using equation (1) Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Note: MAE is the mean absolute error, RMSE is the root mean square error, PB is the percentage bias, σˆ2total is the total residual variance, σˆ2plot is the variance between plots, σˆ2tree is the variance between trees, and σˆ2ε is the variance within trees. Table 2 Parameter estimates with their standard error in parentheses, summary of fit statistics, and correlation structure of both datasets using equation (1) Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Parameter Dataset 1 (1970s) Dataset 2 (2010s) Fitting β0 4.06 (0.055) 5.39 (0.222) β1 0.05 (0.0002) 0.05 (0.0006) MAE (mm) 2.85 3.30 RMSE (mm) 3.72 4.25 PB (%) 15.42 16.01 Validation (10-fold CV) MAE (mm) 2.79 3.37 RMSE (mm) 3.54 4.24 PB (%) 15.72 16.36 σˆ2total 12.88 17.93 σˆ2plot/σˆ2total (%) 33.01 31.40 σˆ2tree/σˆ2total (%) 27.43 33.31 σˆ2ε/σˆ2total (%) 39.56 35.29 Note: MAE is the mean absolute error, RMSE is the root mean square error, PB is the percentage bias, σˆ2total is the total residual variance, σˆ2plot is the variance between plots, σˆ2tree is the variance between trees, and σˆ2ε is the variance within trees. Figure 3 View largeDownload slide Results from the equivalence test of equation (1) for the 16 different sample size combinations. Those shown are proportion of runs (out of 1000) that resulted in equivalent (P < 0.05, region of equivalence = 1 mm) predictions from the full training dataset (larger values are better). Labels on the x-axis refer to the number of plots and the number of trees per plot, respectively. Figure 3 View largeDownload slide Results from the equivalence test of equation (1) for the 16 different sample size combinations. Those shown are proportion of runs (out of 1000) that resulted in equivalent (P < 0.05, region of equivalence = 1 mm) predictions from the full training dataset (larger values are better). Labels on the x-axis refer to the number of plots and the number of trees per plot, respectively. Model evaluation Testing the different model forms from the literature showed that equation (6) could best describe both datasets in terms of MAE, RMSE and PB (Table 3). Blocked cross-validation revealed that equations (5) and (9) were best for predictions for datasets 1 and 2, respectively. Large errors in the validation process showed that equation (7) was not useful for predicting double bark thickness of Silver fir. Parameter estimates of all tested models can be found in Table 4. The prediction bias (expressed by the residuals) of equation (5) was equivalent for the two datasets (P < 0.05, region of equivalence: 0.5 mm), showing that the model could describe both datasets equally well. However, bark thickness was larger in dataset 1. Predictions for the newly assessed data (dataset 2) were significantly higher (P < 0.05, region of equivalence: 0.5 mm) when the model had been fit on dataset 1 compared with predictions made with a model fit from dataset 2. Table 3 Summary of fit statistics for mixed-effects models (4)–(9) and the respective prediction bias using only the fixed terms Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Best values per dataset are printed in bold. Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; MAE = mean absolute error; RMSE = root mean square error; PB, percentage bias. Table 3 Summary of fit statistics for mixed-effects models (4)–(9) and the respective prediction bias using only the fixed terms Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Dataset Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970 s) Fitting MAE (mm) 1.90 1.90 1.29 1.77 1.89 1.48 RMSE (mm) 2.64 2.62 1.77 2.48 2.61 2.01 PB (%) 9.98 9.93 6.92 9.27 9.95 7.90 Validation (10-fold CV) MAE (mm) 3.00 2.51 3.47 4.74 3.00 2.60 RMSE (mm) 3.77 3.23 4.49 5.63 3.76 3.38 PB (%) 16.88 13.90 19.78 35.88 16.79 14.64 2 (2010 s) Fitting MAE (mm) 2.12 2.29 1.50 2.08 1.98 1.83 RMSE (mm) 2.87 3.06 1.98 2.80 2.64 2.36 PB (%) 10.05 10.94 7.42 9.85 9.36 9.04 Validation (10-fold CV) MAE (mm) 3.58 3.03 3.86 7.38 3.62 3.01 RMSE (mm) 4.38 3.83 4.74 8.15 4.44 3.77 PB (%) 17.48 14.76 18.99 60.21 17.51 14.91 Best values per dataset are printed in bold. Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; MAE = mean absolute error; RMSE = root mean square error; PB, percentage bias. Table 4 Parameter estimates with their standard error in parentheses of equations (4–9) for both datasets Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; δ = random effect variance; CAR1 = first-order continuous autoregressive term; σ = residual standard error. Table 4 Parameter estimates with their standard error in parentheses of equations (4–9) for both datasets Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Dataset Parameter Equation (4) Equation (5) Equation (6) Equation (7) Equation (8) Equation (9) LME LME LME LME LME NLME 1 (1970s) β0 7.198 (0.3779) 4.317 (0.2939) −0.419 (0.0470) 7.199 (0.3977) −3.980 (0.0415) β1 0.040 (2.5E−04) 0.009 (3.6E−04) −4.69E−04 (3.1E−05) 0.055 (0.0012) 0.037 (0.0016) 0.674 (0.0127) β2 0.038 (2.5E−04) 0.606 (0.0084) −0.010 (8.41E−04) 2.4E−05 (5.03E−06) 5.057 (0.0936) β3 0.060 (0.0015) 0.034 (0.0012) −3.8E−08 (4.68E−09) 1.423 (0.0224) β4 −2.55E−04 (3.4E−05) 0.016 (3.3E−04) β5 9.2E−04 (3.3E−04) β6 0.254 (0.0195) δ 0.692 0.698 −0.098 0.713 0.671 −0.091 CAR 1 0.914 0.910 0.857 0.903 0.913 0.865 σ 0.053 0.050 0.183 0.044 0.059 0.192 2 (2010s) β0 8.641 (0.4824) 2.403 (0.5834) −0.390 (0.1343) 8.410 (0.6684) −5.261 (0.1841) β1 0.035 (7.1E−04) 0.013 (0.0010) −3.7E−04 (8.4E−05) 0.032 (0.0045) 0.034 (0.0047) 1.280 (0.0832) β2 0.034 (7.1E−04) 0.606 (0.0274) −0.003 (0.0030) 1.9E−05 (1.3E−05) 3.021 (0.1699) β3 0.039 (0.0063) 0.013 (0.0042) −2.8E−08 (9.6E−09) 2.270 (0.1070) β4 1.6E−04 (1.3E−04) 0.016 (9.1E−04) β5 0.007 (0.0011) β6 0.511 (0.1260) δ 0.520 0.538 −0.277 0.526 0.478 −0.308 CAR 1 0.898 0.897 0.824 0.891 0.886 0.839 σ 0.148 0.132 0.548 0.140 0.178 0.751 Note: LME = linear mixed-effects model; NLME = nonlinear mixed-effects model; δ = random effect variance; CAR1 = first-order continuous autoregressive term; σ = residual standard error. Geospatial effects The GAMM fitting of equation (10) and its sub-models obtained the best fit (BIC = 12,109.4) when only the tree-inherent variables dob, dbhob, h/H and Age were included as predictors, which is why only these predictors were kept for further analysis. Removing h/H as well only slightly reduced fit quality (BIC = 12,119.1). Total tree height, geographic coordinates, elevation and the climatic factors of temperature and precipitation could not significantly contribute to explaining the variation of bark thickness. The estimated partial effects of the remaining predictors on double bark thickness are shown with their smooth curves in Figure 4. Double bark thickness showed an increasing trend with increasing dob (Figure 4a) and with dbhob (Figure 4b). The effect of h/H was small, but obvious at the lowest measurement point at breast height (Figure 4c). A significant increase of DBT with increasing age was observed (Figure 4d). Figure 4 View largeDownload slide Estimated smooth functions describing the partial effect of the four covariates on double bark thickness (DBT): (a) the smooth of diameter outside bark, (b) the smooth of dbh outside bark, (c) the smooth of relative tree height, and (d) the smooth of age. Shaded areas are 95%-confidence intervals. Figure 4 View largeDownload slide Estimated smooth functions describing the partial effect of the four covariates on double bark thickness (DBT): (a) the smooth of diameter outside bark, (b) the smooth of dbh outside bark, (c) the smooth of relative tree height, and (d) the smooth of age. Shaded areas are 95%-confidence intervals. Bark thickness equations for roundwood volume determination Relative bark thickness was quite constant in the lower three quarters of tree height and increased slightly further up in the tree. Depending on the bucking pattern with varying log length, midpoints of logs are positioned at different relative tree heights. Nevertheless, results show that bark thickness equations referring to mid-diameters of long logs and of 5-m segments differ only slightly from each other (Table 5, Figure 5). When the equation coefficients for long logs were used as opposed to coefficients for 5-m segments, the calculated log volume was 0.8 and 0.7 per cent higher for long logs (n = 192) and 5-m segments (n = 807), respectively. The volumes of the individual logs were equivalent (paired TOST, P < 0.05, region of equivalence 0.01 m3) using the two sets of coefficients for the 5-m segments, but not equivalent for the long logs. Predicted bark thickness for small ends of 5-m segments with a diameter of up to 40 cm was higher than for the other tested reference diameters (Figure 5). Table 5 Best fit model coefficients with their standard error in parentheses and coefficient of determination of the functions (1) or (11) predicting double bark thickness with different reference diameters for datasets 1 and 2. All sampling locations means every 2 m in each tree. Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Table 5 Best fit model coefficients with their standard error in parentheses and coefficient of determination of the functions (1) or (11) predicting double bark thickness with different reference diameters for datasets 1 and 2. All sampling locations means every 2 m in each tree. Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Dataset Reference diameter Equation β0 β1 β2 R2 1 long logs: midpoint (11) −0.836 (0.5118) 0.756 (0.0341) −0.002 (0.0005) 0.79 2 all sampling locations (1) 5.322 (0.2330) 0.451 (0.0064) – 0.68 2 long logs: midpoint (1) 0.946 (0.8830) 0.563 (0.0241) – 0.74 2 5-m segments midpoint (1) 3.848 (0.4656) 0.498 (0.0123) – 0.67 2 5-m segments small end (1) 5.691 (0.3365) 0.463 (0.0104) – 0.71 Figure 5 View largeDownload slide Regression lines describing double bark thickness against the rounded diameter using different reference diameters, which were gained from different bucking variants of virtual logs produced from dataset 2. For comparison, also all sampled locations (every 2 m along the tree bole) are used as reference diameters. Figure 5 View largeDownload slide Regression lines describing double bark thickness against the rounded diameter using different reference diameters, which were gained from different bucking variants of virtual logs produced from dataset 2. For comparison, also all sampled locations (every 2 m along the tree bole) are used as reference diameters. Bark thickness at midpoints of long logs was smaller for virtual logs of dataset 2 than those of dataset 1 (Table 5, Figure 6). In Germany, the currently suggested boundaries for the diameter classes referring to bark deduction values of 2, 3 and 4 cm have been developed using dataset 1 and are 23, 39 and 56 cm, respectively (Anonymous, 2015). The developed bark thickness equation from dataset 2 would lead to new class boundaries at 25, 43 and 61 cm. Assuming that harvesting intensity and assortments are similar to those of the last 5 years, the application of the new class boundaries would mean a 0.8 per cent annual increase in sales volume of manually measured Silver fir logs. Figure 6 View largeDownload slide Double bark thickness plotted against the rounded mid-diameter with 95% confidence intervals of the regression lines. Vertical lines indicate the class boundaries of integer bark deduction values (horizontal lines: class midpoints 10–20, 30–40 and 40–50 mm dbt). Shifts are resulting from more shallow-barked Silver firs in the more recent data (1970s vs 2010s). Figure 6 View largeDownload slide Double bark thickness plotted against the rounded mid-diameter with 95% confidence intervals of the regression lines. Vertical lines indicate the class boundaries of integer bark deduction values (horizontal lines: class midpoints 10–20, 30–40 and 40–50 mm dbt). Shifts are resulting from more shallow-barked Silver firs in the more recent data (1970s vs 2010s). Discussion The development of accurate bark thickness equations depends on how well within-tree and between-tree variability is captured by an appropriate sampling design, adequate sample sizes, and the predictor variables in the chosen model. This study showed that the sampling design with 2-m increments between the sampling locations starting from 1.3 m above ground was adequate to assess the variability within a stem in most cases. Thus, this design is recommended for further studies of species with similar bark morphology. The sampling design with regularly spaced sampling locations (every 2 m along the tree bole) will lead to more measurements with increasing tree height. In contrast, a fixed number of relative positions along the tree bole, as is suggested in other studies (e.g. Korell, 1972; Feduccia and Mann, 1976; Kozak and Yang, 1981; Gordon, 1983; Laasasenaho et al., 2005), could reduce the support for each data point and would require additional effort in the field to recompute relative positions for each single tree. As shown by Stängle et al. (2016b), the type of bark gauge used in this study overestimated bark thickness by 0.52 mm (SD = 1.59 mm) when compared with CT-derived measurements on Norway spruce logs. This could potentially lead to biased bark thickness predictions for Silver fir as well. If a more accurate measurement device was used, a smaller number of measurements per location would be sufficient to reach the same level of accuracy. The observed high intraclass correlation within the residual variance of equation (1) shows that the relationship between diameter and bark thickness varies strongly between individual trees and plots. This was also supported by the Monte-Carlo simulations representing different sample sizes. The number of sampled plots strongly determined the number of iterations with predicted bark thickness values that were equivalent to predictions from the full-data model (Figure 3). Model choice is an important factor influencing the response in predictive capacity of bark thickness equations on sample size. For Norway spruce it was shown that a more complex model with a higher fit quality required fewer sampled plots and trees per plot (Stängle et al., 2016b). The chosen method of drawing plots and trees with replacement from the training data was necessary to account for differing numbers of trees per sampled plot and differing diameter distributions between plots. The variation observed cannot, therefore, be attributed only to sample size but also to the data-derived uncertainty. For that reason, the suggested sample sizes might slightly overestimate those actually required to achieve the intended precision of predictions. Results showed that a rather complex nonlinear model (equation (9)) and quite a simple model (equation (5)) were most suited for predicting bark thickness of Silver fir. Model (6) described both datasets best, but it failed to predict bark thickness as precisely as models (5) and (9). The same results were reported for Norway spruce in the same region (Stängle et al., 2017). Depending on the purpose of application, different explanatory variables are available and different degrees of computational complexity are preferable. If bark thickness is to be modelled for inventory data, for example, a complex model with many explanatory variables, such as model (9) is suggested because taper curves outside bark are often available. If, in contrast, total and relative tree height are not known, as in the situation when a harvester measures a tree, the best possible model without these covariates should be selected (equation (5)). Effects of total tree height were generally absent and relative tree height only marginally explained variation in bark thickness. Equation (7), which is mainly based on relative and total tree height as predictors, cannot be recommended for Silver fir, although it was proposed for different Pinus species (Cao and Pepper, 1986). The high influence of absolute and relative tree height on bark thickness variation for pines was also confirmed by Murphy and Cown (2015), who reported bark volume variation within trees, between trees, and between sites for Pinus radiata (D. Don). For Norway spruce, effects of tree height and relative tree height on bark thickness variation were shown as well (Stängle et al., 2017). As Norway spruce and most Pinus species are less shade tolerant than Silver fir, a higher correlation between total tree height and age can be expected for those species. Therefore, the effect of total tree height for those species could, at least partially, be caused by collinearity. The variation observed between the sample plots could not be linked to large-scale spatial patterns, such as climatic zones and effects of geographic coordinates were generally absent. Moreover, the climatic factors of temperature and precipitation as well as elevation could not explain differences between plots. Tree age, however, was found to have a significant positive effect on bark thickness. These findings are in accordance to those of Laasasenaho et al. (2005), Sonmez et al. (2007) and Stängle et al. (2017) for Norway spruce. This leads to the conclusion that higher growth rates result in smaller relative bark thickness for those species. Other studies also show that increased site quality, which can usually be linked to faster growth of individual trees, reduced bark thickness of Silver fir (Božić et al., 2007). It seems that no single environmental predictor can explain differences in bark thickness. Rather, the combination of many factors influencing the productivity of a stand and the growth of individual trees is required. For future studies of bark thickness, the assessment of tree age or site index is, therefore, strongly recommended. The observation that faster growing trees have thinner bark might explain the observed thinner bark in the more recent measurements of dataset 2. Since growth of Silver fir in southwestern Germany has slightly accelerated in the twentieth century (Kohnle et al., 2014), it can be expected that the trees of the new dataset have grown faster than those trees measured 40 years earlier. The small difference in relative bark thickness along the stem was apparent when bark thickness equations were fitted on mid-diameters of long logs and of 5-m segments. As there were hardly any differences in volume calculation, the same equation coefficients could be used for different log assortments. However, the small ends of 5-m segments showed larger relative bark thickness as they are positioned further up in the tree. This difference, especially for small diameters, suggests that different equation coefficients should be used for bucking optimization that is based on minimum small end diameters. Alternatively, a model, which includes information on the relative height of the sampling location, could be chosen. Thinner bark in the more recent measurements resulted in different diameter classes for integer bark deduction for the two datasets. The values that are currently in use in Germany (Anonymous, 2015) are based on the older measurements. Therefore, they tend to overestimate bark thickness and lead to biased inside-bark volume estimations. Conclusion To parameterize bark thickness equations, within-tree and between-tree variation in Silver fir bark thickness requires several sampling locations per tree and a minimum of five sampled trees in at least 35 plots in the study region. Bark thickness of Silver fir can be described well using tree diameter at the sampling location and the breast height diameter as predictor variables. Faster growing trees were shown to have relatively thinner bark. From the demonstrated influence of tree age – and therefore growth rate – on bark thickness, it can be inferred that changing growth conditions can alter bark allometry, which is why bark thickness equations should be regularly reviewed for their validity. For measuring roundwood, the same bark deduction values can be suggested for long logs (6–21 m) and 5-m segments. For bucking optimization based on minimum small end diameters, different equation coefficients that reflect higher relative bark thickness at larger relative tree height are recommended. The findings show the need to update the current bark deduction values for manual log measurement in the state of Baden-Württemberg in order to estimate log volume more accurately. As the same values are currently used in large parts of Central Europe, further validation measurements in other federal states in Germany and in neighbouring countries are suggested. Acknowledgements The methods described in this manuscript were partially developed with support from Professor Aaron Weiskittel (University of Maine) and additional statistical advice from Dr Gerald Kändler (FVA Baden-Württemberg) was appreciated. We would like to thank two anonymous reviewers, the Editor-in-Chief Dr Gary Kerr, and especially the handling Editor Professor Alexis Achim for their comments and suggestions that helped to improve previous versions of this manuscript. Conflict of interest statement None declared. Funding This study was funded by the Ministry for Rural Affairs and Consumer Protection Baden-Württemberg, Germany (MLR) and was conducted at the Forest Research Institute Baden Württemberg (FVA). References Altherr , E. , Unfried , P. , Hradetzky , J. and Hradetzky , V. 1974 Statistische Rindenbeziehungen als Hilfsmittel zur Ausformung und Aufmessung unentrindeten Stammholzes: Teil I: Kiefer, Buche, Hainbuche, Esche und Roterle. In Mitteilungen der Forstlichen Versuchs- und Forschungsanstalt, Forstliche Versuchs- und Forschungsanstalt Baden-Württemberg, Freiburg im Breisgau, Germany, pp. 137. Altherr , E. , Unfried , P. , Hradetzky , J. and Hradetzky , V. 1978 Statistische Rindenbeziehungen als Hilfsmittel zur Ausformung und Aufmessung unentrindeten Stammholzes: Teil IV: Fichte, Tanne, Douglasie und Sitka-Fichte. In Mitteilungen der Forstlichen Versuchs- und Forschungsanstalt, Forstliche Versuchs- und Forschungsanstalt Baden-Württemberg, Freiburg im Breisgau, Germany, pp. 294. Anonymous . 2015 Rahmenvereinbarung für den Rohholzhandel in Deutschland (RVR). Deutscher Forstwirtschaftsrat e.V.; Deutscher Holzwirtschaftsrat e.V. http://www.saegeindustrie.de/rvr/docs/dynamisch/6205/rvr_gesamtdokument_2.auflage_stand_oktober_2015.pdf (accessed on 14 July, 2016) Antony , F. , Schimleck , L.R. , Daniels , R.F. , Clark , A. , Borders , B.E. , Kane , M.B. , et al. 2015 Whole-tree bark and wood properties of loblolly pine from intensively managed plantations . For. Sci. 61 , 55 – 66 . Barton , K. 2016 MuMIn: Multi-model inference: R package version 1.15.6. Baskerville , G.L. 1972 Use of logarithmic regression in the estimation of plant biomass . Can. J. For. Res. 2 , 49 – 53 . Google Scholar CrossRef Search ADS Božić , M. , Čavlović , J. , Vedriš , M. and Jazbec , M. 2007 Modeling bark thickness of silver fir trees (Abies alba Mill.) . Šumarski list 131 , 3 – 12 . Cao , Q.V. and Pepper , W.D. 1986 Predicting inside bark diameter for Shortleaf, Loblolly, and Longleaf pines . South. J. Appl. For. 10 , 220 – 224 . Cellini , J.M. , Galarza , M. , Burns , S.L. , Martinez-Pastur , G.J. and Lencinas , M.V. 2012 Equations of bark thickness and volume profiles at different heights with easy-measurement variables . Forest Sys. 21 , 23 – 30 . Google Scholar CrossRef Search ADS Dimitrov , E.T. 1976 Mathematical models for determining the bark volume of spruce in relation to certain mensurational characteristics (translated from Hungarian) . Gorstoskopanska Nauka (Forest Science) 13 , 52 – 63 . Doruska , P.F. , Patterson , D. , Hartley , J. , Hurd , M. and Hart , T. 2009 Newer technologies and bioenergy bring focus back to bark factor equations . J. For. 107 , 38 – 43 . EUFORGEN . 2009 Distribution map of Silver fir (Abies alba). http://www.euforgen.org (accessed on 14 July, 2016). Feduccia , D.P. and Mann , W.F. 1976 Bark thickness of 17-year-old loblolly pine planted at different spacings. USDA For. Serv., p. 2. Fonseca , M.A. 2005 The Measurement of Roundwood: Methodologies and Conversion Ratios . CABI , p. 269 . Google Scholar CrossRef Search ADS Freund , H. and Grehn , J. 1951 Handbuch der Mikroskopie in der Technik. 1 - Bd. 5, Mikroskopie des Rohholzes und der Rinden . Umschau-Verlag , p. 456 . Gordon , A. 1983 Estimating bark thickness of Pinus radiata . New Zeal. J. For. Sci 13 , 340 – 348 . Hannrup , B. 2004 Funktioner för skattning av barkens tjocklek hos tall och gran vid avverkning med skördare. (Functions for prediction of bark thickness of Norway spruce and Scots pine at CTL-harvesting). Skogforsk, Uppsala, p. 33. Hoffmann , J. 1958 Untersuchungen über die Rindenstärke der Fichte auf verschiedenen Standorten im südöstlichen Thüringer Wald . Wiss Z. TU Dresden 7 , 361 – 368 . Kohnle , U. , Albrecht , A. , Lenk , E. , Ohnemus , K. and Yue , C. 2014 Growth trends driven by environmental factors extracted from long term experimental data in southwest Germany . AFJZ 185 , 97 – 117 . Kohnle , U. , Hein , S. , Sorensen , F.C. and Weiskittel , A.R. 2012 Effects of seed source origin on bark thickness of Douglas-fir (Pseudotsuga menziesii) growing in southwestern Germany . Can. J. For. Res. 42 , 382 – 399 . Google Scholar CrossRef Search ADS Korell , U. 1972 Über Rindendicken der Fichte . Beitr. Forstwirtsch 2 , 54 – 55 . Kozak , A. and Yang , R.C. 1981 Equations for estimating bark volume and thickness of commercial trees in British Columbia . For. Chron. 57 , 112 – 115 . Google Scholar CrossRef Search ADS Kublin , E. and Breidenbach , J. 2015 TapeR: Flexible tree taper curves based on Semiparametric Mixed Models: R package version 0.3.3. Laasasenaho , J. , Melkas , T. and Aldén , S. 2005 Modelling bark thickness of Picea abies with taper curves . For. Ecol. Manage. 206 , 35 – 47 . Google Scholar CrossRef Search ADS Li , R. and Weiskittel , A.R. 2011 Estimating and predicting bark thickness for seven conifer species in the Acadian Region of North America using a mixed-effects modeling approach: comparison of model forms and subsampling strategies . Eur. J. Forest. Res. 130 , 219 – 233 . Google Scholar CrossRef Search ADS Loetsch , F. , Zöhrer , F. and Haller , K.E. 1973 Forest Inventory Vol. II . BLV Verl. Ges , p. 436 . Marshall , H.D. , Murphy , G.E. and Lachenbruch , B. 2006 Effects of bark thickness estimates on opimal log merchandising . For. Prod. J 56 , 87 – 92 . Martin , R.E. and Crist , J.B. 1970 Elements of bark structure and terminology . Wood Fiber 2 , 269 – 279 . Murphy , G. and Cown , D. 2015 Within-tree, between-tree, and geospatial variation in estimated Pinus radiata bark volume and weight in New Zealand . New Zeal. J. For. Sci. 45 , 18 . Google Scholar CrossRef Search ADS Niklas , K.J. 1999 The mechanical role of bark . Am. J. Bot. 86 , 465 – 469 . Google Scholar CrossRef Search ADS PubMed Persson , B. and Downie , B. 1992 Variation in bark thickness of young Pinus contorta var. latifolia Engelm. in Sweden . Scand. J. For. Res. 7 , 99 – 106 . Google Scholar CrossRef Search ADS Pinheiro , J. , Bates , D. , DebRoy , S. , Sarkar , D. and Team , R.C. 2016 nlme: Linear and Nonlinear Mixed Effects Models: R package version 3.1–124. Pinheiro , J.C. and Bates , D.M. 2000 Mixed-effects Models in S and S-PLUS . Springer , p. 528 . Google Scholar CrossRef Search ADS Pizzi , A. 2016 Wood products and green chemistry . Ann. For. Sci. 73 , 185 – 203 . Google Scholar CrossRef Search ADS R Core Team . 2015 R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing . Robinson , A.P. 2016 equivalence: Provides tests and graphics for assessing tests of equivalence: R package version 0.7.1. Schielzeth , H. and Nakagawa , S. 2013 Nested by design: model fitting and interpretation in a mixed model era . Methods Ecol. Evol. 4 , 14 – 24 . Google Scholar CrossRef Search ADS Schmidt-Vogt , H. 1986 Die Fichte–Ein Handbuch in zwei Bänden, Band II/1: Wachstum, Züchtung, Boden, Umwelt, Holz . Parey , p. 563 . Skogforsk . 2012 Standard for Forest Data and Communication Appendix: Definitions of variables - General and country specific, http://www.skogforsk.se/contentassets/b063db555a664ff8b515ce121f4a42d1/appendix1_eng_120418.pdf (accessed on 27 June, 2016). Sonmez , T. , Keles , S. and Tilki , F. 2007 Effect of aspect, tree age and tree diameter on bark thickness of Picea orientalis . Scand. J. For. Res. 22 , 193 – 197 . Google Scholar CrossRef Search ADS Sprugel , D. 1983 Correcting for bias in log-transformed allometric equations . Ecology 64 , 209 – 210 . Google Scholar CrossRef Search ADS Stängle , S.M. , Sauter , U.H. , Brüchert , F. and Kändler , G. 2016 a Überprüfung der Rindenabzugswerte für Fichten-Stammholz in Baden-Württemberg (A review of bark deduction values for Norway spruce logs in Baden-Württemberg) . Forstarchiv 87 , 162 – 169 . Stängle , S.M. , Sauter , U.H. and Dormann , C.F. 2017 Comparison of models for estimating bark thickness of Picea abies in southwest Germany: the role of tree, stand, and environmental factors . Ann. For. Sci. 74 , 16 . Google Scholar CrossRef Search ADS Stängle , S.M. , Weiskittel , A.R. , Dormann , C.F. and Brüchert , F. 2016 b Measurement and prediction of bark thickness in Picea abies: assessment of accuracy, precision, and sample size requirements . Can. J. For. Res. 46 , 39 – 47 . Google Scholar CrossRef Search ADS Temesgen , H. , Affleck , D. , Poudel , K. , Gray , A. and Sessions , J. 2015 A review of the challenges and opportunities in estimating above ground forest biomass using tree-level models . Scand. J. For. Res. 30 , 326 – 335 . Weiskittel , A.R. , Hann , D.W. , Kershaw , J.A. and Vanclay , J.K. 2011 Forest Growth and Yield Modeling , 415 . Wiley . Google Scholar CrossRef Search ADS Wilhelmsson , L. , Arlinger , J. , Spångberg , K. , Lundqvist , S.-O. , Grahn , T. , Hedenberg , Ö. and Olsson , L. 2002 Models for Predicting Wood Properties in Stems of Picea abies and Pinus sylvestris in Sweden . Scand. J. For. Res. 17 , 330 – 350 . Google Scholar CrossRef Search ADS Wood , S.N. 2006 Generalized Additive Models: An Introduction With R . Chapman & Hall/CRC , p. 392 . Wood , S.N. 2011 Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models . J. Roy. Stat. Soc. Ser. B (Stat. Method.) 73 , 3 – 36 . Google Scholar CrossRef Search ADS Zacco , P. 1974 Barktjockleken hos sågtimmer (The bark thickness of saw logs). In Rapport . The Swedish University of Agricultural Sciences, Department of Forest Products. Stockholm , p. 53 . © Institute of Chartered Foresters, 2017. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Forestry: An International Journal Of Forest ResearchOxford University Press

Published: Oct 27, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off