Modelling bark volume for six commercially important tree species in France: assessment of models and application at regional scale

Abstract

Key message

A set of models of bark thickness at breast height and bark volume are now available for six species in France. A common model suitable for predicting bark volume was proposed for all species. A small but significant altitude effect on bark thickness at breast height was detected for three species.

Context

The growing demand for wood energy and bio-molecules requires a thorough evaluation of forest biomass, particularly bark.

Aims

The objective of this study is to have statistical models of bark volumes for the six main forest species present in North-Eastern France and to be able to estimate regional bark biomasses and quantities of chemical extractives at regional scale.

Methods

A large databank gathering bark thickness measured at different heights in France was used for selecting literature or new alternative models of tree bark volume. These models were applied to the available forest inventory data from North-Eastern France to estimate the regional bark volume. Secondly, by multiplying these volumes by basic density data and extractive content recently obtained, bark biomasses and extractives quantities were deduced.

Results

The first results consist in a set of species-specific models of bark thickness at breast height with R² around 0.70 and a relative RMSE around 30% which is an improvement of 0.1 for R² and of 1–2% for relative RMSE depending on the species compared to the best models from the literature. The second results consist in a set of species-specific models of tree bark volumes with R² of 0.90 and a relative RMSE which varies between 22% when bark thickness at breast height is included and 40% when it is predicted. A significant relationship between bark thickness at breast height and altitude was also observed. The bark resources of Grand Est and Bourgogne-Franche-Comté regions were estimated at 558 000 m³/year and 611 000 m³/year respectively representing between 5.5% and 15% of the stem volume depending on the species. The propagation of the measurement error of bark gauge was estimated at 5% for model of bark thickness at breast height and 24% for bark volume model.

Conclusion

These results constitute an important contribution for a better knowledge of the bark resource at a regional scale and may help to optimise bark valuation by the forest-wood sector.

1 Introduction

Bark is a multifunctional structure absolutely necessary for the tree life (Rosell 2019). It protects living tissues, i.e., sapwood, phloem, dormant buds (Charles-Dominique et al. 2015), from high temperatures (Pausas 2015; Rosell 2016) or very low ones (De Antonio et al. 2020), superficial injuries due to rock falls, harvesting, herbivores and insects (Theander 1985; Harun and Labosky 1985), or pathogenic micro-organisms (Franceschi et al. 2005). The bark also has less known functions such as storage of water (Levia and Herwitz 2005) or mineral reserves (Schowalter and Morrell 2002), and is involved in tree biomechanics (Clair et al. 2019). Consequently, bark participates through its multiple functions to the adaptation of trees to different ecosystems (Rosell et al. 2014).

To fulfil these different functions, bark can be described by different variables: Its thickness at different heights (Gordon 1983), its total volume in a tree (Wehenkel et al. 2012; Rosell et al. 2017), its proportion of the volume of the stem (Cellini et al. 2012), its density (Miles and Smith 2009), its structure, particularly different in young or old bark (Dedrie et al. 2015), its yields and quantities of different chemical compounds (Jyske et al. 2014; Trivelato et al. 2016; Feng et al. 2013; Brennan et al. 2020), especially carbon (Jones and O’Hara 2018; Castaño-Santamaría and Bravo 2012) and minerals (Buamscha et al. 2007). All these variables not only participate in relating structures and functions in living trees but also are critical in determining the value of harvested timber.

To fulfil these different functions, bark can be described by different variables: Its thickness at different heights, its total volume in a tree, its proportion of the volume of the stem, its density, its structure, particularly different in young or old bark, its yields and quantities of different chemical compounds, especially carbon and minerals. All these variables not only participate in relating structures and functions in living trees but also are critical in determining the value of harvested timber.

To fulfil these different functions, bark can be described by different variables: Its thickness at different heights, its total volume in a tree, its proportion of the volume of the stem, its density, its structure, particularly different in young or old bark, its yields and quantities of different chemical compounds, especially carbon and minerals. All these variables not only participate in relating structures and functions in living trees but also are critical in determining the value of harvested timber.

The bark is thus a rich and very accessible raw material and human has exploited it for a very long time (Pasztory et al. 2016; Harkin and Rowe 1971) to produce for instance medicines (Turner and Hebda 1990; Anderson 1955; Rastogi et al. 2015), perfume (e.g., cinnamon), latex corks (Thomas et al. 1995) or insulating material (Gil 2014). In the current economy of the forest-wood sectors, bark is a by-product of the primary wood processing industries which contributes very little to their total turnover. Its economic value, mostly coming as a horticultural substrate or as a fuel (Lu et al. 2006), may reach 15% for all by-products pooled together (Chalayer 2015).

Not all softwood bark is suitable for use as horticultural substrate and, compared to wood, bark is a second choice fuel due to its high humidity and high mineral composition (Adler 2007) despite its high extractive content which improves its calorific value (Telmo and Lousada 2011; Tenorio and Moya 2013; Fuwape 1989). Only boilers specifically designed to burn biomass with a high content of water, impurities (soil, gravel, etc.), and mineral matter (and therefore ash) can burn it efficiently (Martin 2015).

It turns out that society is increasingly demanding biomolecules that are considered less harmful to health than petro-sourced products and are produced in a more environmentally friendly manner. Industries want to take advantage of this trend. They need consequently to answer questions about the volume of the different markets, the availability of raw materials rich in these biomolecules and the industrial processing of the raw material.

To estimate the available raw material, it is possible to make advantage of data and information already available. Bouvet and Deleuze (2013) have provided values of bark percentage to calculate stem bark volume from stem volume. There are also data on basic density to convert fresh wood volumes into dry biomasses (Billard et al. 2020). Moreover, the rate of extractable chemical compounds (Brennan et al. 2020) expressed in grams of extractives/gram of dry matter can be used to obtain the mass of extractives.

Depending on the intended purpose, bark quantity has been estimated either from bark thickness (Stängle et al. 2017; Muhairwe 2000; Laasasenaho et al. 2005), bark area (or proportion in bark area) or from whole stem bark biomass (e.g., Zianis et al. 2005). There are numerous studies that have developed bark thickness models (e.g., Gordon 1983; Muhairwe 2000; Hannrup 2004; Van Laar 2007). However, few studies have provided models to directly predict the volume of bark from tree variables such as diameter at breast height, tree height, and bark thickness at breast height when it is available (e.g., Kozak and Yang (1981)). Often cited, Meyer (1946) computed the bark volume of a tree on the basis of a constant proportionality between over-bark and under-bark stem diameters, generally measured at breast height. More recently, Wehenkel et al. (2012) and Liepiņš et al. (2015) have modelled the proportion of bark volume.

Modelling bark volume (B_v) and bark thickness at breast height (BTBH) for France is the central objective of this article. On the one hand, when models exist but are not adapted to the French resource they must be readjusted as advised by Jenkins et al. (2003) and Stängle et al. (2017). On the other hand, if the models do not exist, new ones must be developed. In addition, models of BTBH will also be designed since BTBH can be used as an input variable in bark volume models, yet they are not always measured in the field. Numerous bark data have been collected by French research and development organisations: INRAE (Institut National de Recherche pour l’Agriculture, l’Alimentation et l’Environnement), FCBA (institut technologique Forêt Cellulose Bois-construction, Ameublement), ONF (Office National des Forêts), IGN (Institut National de l’Information Géographique et Forestière). This database is managed by the FCBA (Bouvet and Deleuze 2013).

As a secondary objective, we applied the best models to forest inventory data to provide an initial estimate of the available resource, e.g. for the chemical industry. The resource must be assessed in terms of bark volume, bark biomass, and quantity of extractives. We provide this estimate at the scale of two French regions: Grand-Est and Bourgogne Franche-Comté. For the sake of clarity, we will consider only the total content of extractives and not the content of specific chemical families.

This study concerns six temperate tree species that are present and industrially processed in Eastern France: silver fir (Abies alba Mill.), Norway spruce (Picea abies L.), Douglas fir (Pseudotsuga menziesii Mirb.), European beech (Fagus sylvatica L.), sessile oak (Quercus petraea Matt.), and pedunculate oak (Quercus robur L.).

2 Materials and methods

BTBH models were developed from the data recorded during several campaigns of the IGN. Since 2008, these measurements are no longer carried out as the database has been considered complete. Models for total bark volume (B_v) were developed from a “research” dataset called EMERGE. Both datasets are described in detail in Section 2.1.

Table 1 summarises the abbreviations used in this paper.

Table 1 List of abbreviations

2.1 Datasets

In this work, three different databases were used namely EMERGE, IGN, and French NFI.

2.1.1 EMERGE dataset

EMERGE was a project led by the French National Forest Office (ONF) and supported by the French National Research Agency (ANR) (Deleuze et al. 2013). Its purpose was to estimate the available biomass in French forests. Eight French research and development organisations worked together on this project. The EMERGE dataset combined several subsets of data collected by the various partners. This dataset contained bark thickness (B_t) measurements performed in several French regions (Fig. 1), statistically representative of the French resource. The measurements were made using a Swedish bark gauge at several heights along the stem but often not at breast height. In the latter case, a linear interpolation between the two values collected from closest values to 1.3 m allowed for an estimation at this height.

Fig. 1

Location of the sites where the EMERGE data were recorded

Usual tree variables (total tree height, H_tot and diameter at breast height, DBH) were measured and additional information about the corresponding plots such as geographic coordinates and altitude were recorded. The tree DBH distribution is given in Appendix Fig. 5.

The composition of EMERGE database is given in Table 2.

Table 2 Composition of EMERGE dataset

2.1.2 IGN dataset

The second dataset, called IGN is a “resource” dataset. BTBH was measured on a large number of trees everywhere in France. This dataset includes only one measurement at breast height made with a Swedish bark gauge as well as DBH, H_tot and location of the trees (latitude, longitude and altitude). Table 3 shows the composition of the IGN dataset.

Table 3 Composition of IGN dataset

2.1.3 French NFI data in grand Est and Bourgogne-Franche-Comté

The National Forest Inventory (NFI) is a continuous statistical survey of French metropolitan forests, undertaken by the IGN. The NFI is carried out in public and private forests, regardless of whether they are available for wood supply. The NFI design features a systematic sampling grid with squared cells of 1 km side (Colin et al. 2017; Hervé 2016). Each year, 10% of the cells are sampled according to two phases. In the first phase, approximately 80 000 photo plots are interpreted to assess land cover and land use. In the second phase, approximately 7 000 temporary ground plots are established at a sub-sample of first phase photo plot locations that have forest land use. On each plot, every tree with a DBH of at least 7.5 cm is measured. Measurements include DBH and H_tot. Typically, five annual NFI samples are combined to calculate statistics for the forest resource. Moreover, since 2010, the French NFI provides direct measurements of logging based on a re-inventory of temporary plots placed in the inventory 5 years earlier. For example, in 2014, the NFI returned data on plots for the inventory in 2009, which can serve as an objective for giving an estimate of the mean annual harvest between 2009 and 2014.

French NFI visits approximately 1 650 new inventory plots each year in the Grand Est and Bourgogne-Franche-Comté regions (870 plots and 780 plots, respectively). In this paper, we used the observations performed on cut trees during the 2014–2018 periods to infer the amount of bark volumes harvested each year in the two regions between 2009 and 2018. For that, the best BTBH models developed for each species were applied on all felling trees recorded by the NFI to estimate BTBH from tree DBH and altitude whenever significant. Then, the best B_v models were applied to estimate stem bark volume of felling trees based on estimated BTBH, tree DBH and tree H_tot. Finally, the five annual NFI samples (2014–2018) were combined to estimate the volume of bark harvested each year, between 2009 and 2018, from the forest resource in the two regions.

2.2 Calculating the bark volume

First, all B_t measurements from the EMERGE dataset were used to compute over- and under-bark cross-sectional areas at different heights. Over- and under-bark stem area, respectively A_ob and A_ub at a given height were calculated with Eqs. 1 and 2.

$$ A_{ob} = \frac{D^{2}_{ob}}{4}\cdot{}\pi $$

(1)

$$ A_{ub} = \frac{(D_{ob}-2B_{t})^{2}}{4}\cdot{}\pi $$

(2)

where D_ob is the over-bark stem diameter and B_t the bark thickness.

To calculate the bark volume of each tree, over- or under-trunk portions were considered as a stack of truncated cones with a volume given by Eq. 3 with the tree top portion considered as a cone (thus with A₂ = 0) and the stump portion calculated by extrapolation from the areas calculated from the two lowest measures:

$$ V_{truncated\ cone} = \frac{H_{c}}{3}\cdot{}(A_{1}+\sqrt{A_{1}\cdot{}A_{2}}+A_{2}) $$

(3)

where H_c is the truncated cone height and A₁ and A₂ the area of, respectively, its lower section and upper sections.

The sum of over- and under-bark volumes of all truncated cones gives the total over- and under-bark stem volumes for the tree (Fig. 2a). And as a consequence, their difference gives the bark volume.

Fig. 2

a Sequence of truncated cones modelling a stem. b Simplified geometric calculation of bark volume. BTBH: Bark thickness at breast height; DBH: Diameter at breast height; D0_ob and D0_ub: Diameters at ground level over and under bark; H_tot: total tree height

2.3 BTBH modelling

We first adjusted literature models on the FCBA database and then design alternative models. Finally, we selected the most relevant models based on the lowest Akaike Information Criterion (AIC) and relative root mean square error (RMSErel).

2.3.1 BTBH models from the literature

We selected three models from Wilhelmsson et al. (2002) for Norway spruce and (Pinus sylvestris L.) (Eq. 4); Cao and Pepper (1986) for North-American species (Pinus echinata Mill., Pinus taeda L. and Pinus palustris Mill.) (Eq. 5) and Gordon (1983) for Pinus radiata D. Don in New-Zealand (Eq. 6).

$$ BTBH=exp(a+b\cdot{}DBH+c\cdot{}ln(DBH)) $$

(4)

$$ BTBH=DBH\left( a + b\cdot{}\frac{1.3}{H_{tot}}+c\cdot{}\left( \frac{1.3}{H_{tot}}\right)^{2}+d\cdot{}H_{tot}\right) $$

(5)

$$ BTBH=a+b\cdot{}DBH+c\cdot{}DBH^{2}+d\cdot{}DBH^{3} $$

(6)

2.3.2 Alternative BTBH models

After carefully analysing several equations we selected Eq. 7 as a base for our modelling, reflecting the strong relationship with DBH.

$$ \begin{array}{@{}rcl@{}} BTBH = a\cdot{}DBH^{b} \end{array} $$

(7)

The addition of an intercept to Eq. 7 was studied but we preferred to remove it, even if this parameter was significant, because we considered that a tree with DBH = 0 should logically have BTBH = 0.

For assessing whether a relation between altitude and BTBH could be introduced, quantitative altitude values were transformed into ten increasing altitude classes and their mean was assigned to each tree belonging to their corresponding class. To make sure the classes were approximatively of comparable size, we merged classes representing less than 5% of our data with their closest class. The parameters a and b in Eq. 7 were then adjusted according to altitude class. We thus obtained around ten values of a and b, for each species, that were analysed in relation with the mean altitude of trees from each class. Sometimes and only on a, we observed a significant altitude effect. In this case Eq. 7 becomes Eq. 8.

$$ BTBH = (c\cdot{}alt+d)\cdot{}DBH^{b} $$

(8)

with alt as the altitude

2.4 Bark volume models

2.4.1 Bark volume models from the literature

Meyer (1946) proposed a method to calculate B_v using the a coefficient according to Eq. 9.

$$ a= \frac{\sum{DBH}-2*\sum{BTBH}}{\sum{DBH}} $$

(9)

where the sum is calculated on data from a pool of trees of the same species.

Bark volume for individual tree may therefore be computed using Eq. 10.

$$ B_{v}= V_{ob}\cdot(1-a^{2}) $$

(10)

We determined the a coefficient from the IGN dataset for each species. We then calculated the bark volume of the trees from the EMERGE dataset using the a coefficient for this species and Eq. 10. This method was developed by Meyer (1946) for Cinchona but he also proposed an application for other species such as Tsuga canadensis, Pinus strobus, Quercus alba and Acer rubrum.

Kozak and Yang (1981) proposed a model of bark volume following Eq. 11. This model was applied on 32 species including both hardwood species and softwood species.

$$ B_{v}=a\cdot{}DBH^{c}\cdot{}BTBH^{d}\cdot{}H_{tot}^{e} $$

(11)

2.4.2 Alternative model of bark volume

We built a new model of B_v by considering the difference between two cones V_ob and V_ub which are the over-bark and under-bark stem volumes assuming the trunk is cone-shaped (Eq. 12 and Fig. 2 b).

$$ B_{v} = V_{ob}-V_{ub} = \frac{\pi}{12}\cdot{}H_{tot}\cdot{}(D0^{2}_{ob}-D0^{2}_{ub}) $$

(12)

where D0_ob and D0_ub are the over bark and under bark diameters at the ground level.

D0_ob and D0_ub can be calculated from DBH and BTBH with Eqs. 13 and 14, respectively, on the basis of the Thales’ theorem.

$$ D0_{ob}=DBH\cdot{}\frac{H_{tot}}{H_{tot}-1.3} $$

(13)

$$ D0_{ub}=(DBH-2BTBH)\cdot{}\frac{H_{tot}}{H_{tot}-1.3} $$

(14)

We obtained the theoretical Eq. 15.

$$ B_{v} = \frac{\pi}{3} \cdot \frac{H_{tot}^{3}}{(H_{tot}-1.3)^{2}} \cdot BTBH \cdot \left[ DBH-BTBH \right] $$

(15)

For mature trees, by assuming that H_tot ≫ 1.3 m and DBH ≫ BTBH, Eq. 15 can be simplified to give Eq. 16.

$$ B_{v} = \frac{\pi}{3} \cdot H_{tot} \cdot BTBH \cdot DBH $$

(16)

Models following Eqs. 17 and 18 were designed from Eqs. 15 and 16, respectively.

Nevertheless, Eqs. 15 and 16 are only exact for a perfectly cone shaped stem, both over bark and under bark. To account for the geometric difference with a real stem the two models following Eqs. 17 and 18, in which parameters a and b can be adjusted statistically, were designed.

$$ B_{v} = a\cdot \frac{H_{tot}^{3}}{(H_{tot}-1.3)^{2}} \cdot BTBH \cdot \left[ DBH-BTBH \right] + b $$

(17)

$$ B_{v} = a \cdot H_{tot} \cdot BTBH \cdot DBH+b $$

(18)

Finally, we obtained the model following Eq. 19 by replacing BTBH by its prediction $\widehat {BTBH}$.

$$ B_{v} = a \cdot H_{tot} \cdot \widehat{BTBH} \cdot DBH+b $$

(19)

$\widehat {BTBH}$ is the predicted value obtained with the best BTBH models depending on the species.

2.5 Statistical methods

2.5.1 Modelling

All the regressions were carried out using the R software (R Core Team 2018). Two methods have been applied, depending on the modelled variable.

For BTBH, the R function nls was used to fit the model. A total of 13 outliers were removed from the datapool (0 observation for silver fir and for Norway spruce, 7 for Douglas fir, 2 for European beech, 2 for sessile oak, and 2 for pedunculate oak). These observations were removed based on graphical analysis.

For B_v, bark volume measurements showed strong heteroscedasticity with variance increasing with tree size. For handling heteroscedasticity in the non linear regression analyses, we used the gnls function, provided by the nlme package (Pinheiro et al. 2019) with a variance structure described by the function varPower (Eq. 20):

$$ Var(\epsilon)=\sigma^{2} \cdot (DBH^{2} \cdot H_{tot})^{2\theta} $$

(20)

The Breusch-Pagan test (BP-test) was used to verify the efficiency of this formulation by analysing the heteroscedasticity of the model residuals.

The model has been considered without heteroscedasticity if the obtained p-value is above 0.1.

In order to assess if species-specific models were required especially for distinguishing the two oak species, the R function nlsList provided by the nlme package (Pinheiro et al. 2019) was used.

This function partitions the data according to the levels of a grouping factor, in our case, the species, and gives specific fits for each data partition, using same models for different factors.

2.5.2 Cross-validation

All models, either from the literature or alternative ones, were validated through a cross-validation method as follows: The data were separated according to their plot of origin. Then these parts were merged randomly to form ten groups of approximatively the same size. Each group was used as validation set, the remaining observations being used as training set. A fitted value was thus obtained from a validation set for each measure made. These two datasets were used to calculate RMSE and relative RMSE (RMSErel) in%, as recommended by Mayer and Butler (1993), using Eqs. 21 and 22, respectively.

$$ RMSE=\sqrt{\frac{\sum(y-\widehat{y})^{2}}{n}} $$

(21)

$$ RMSErel=\frac{RMSE}{\bar{y}}\cdot{}100 $$

(22)

where y is the measured value, $\widehat {y}$ the fitted value, $\bar {y}$ the mean of observed values and n the number of observations.

The R² value of the model was obtained by fitting a linear regression between observed and fitted value, as recommended by Piñeiro et al. (2008). This R² was taken to characterize the quality of this modelling.

Then our data were split again, calibration and validation performed twice more. Three values of the indicators RMSE, RMSErel, and R² were thus obtained. Finally, these three values were used to calculate their mean value, which are presented in Section 3. However, the parameter values are issued from the model adjustment made on all datasets.

We also calculated the AIC (Akaike 1973) for each model. We used for this a model adjusted on all data, without cross-validation.

2.6 Propagation of measurement error (PME) for Norway spruce

According to the analysis done by Stängle et al. (2016), measuring B_t using a Swedish bark gauge is subject to an error. This measurement error was identified for Norway spruce species in the same analysis as an overestimation of the B_t by 13.6% ± 28.4% (mean ± standard deviation) relative to its true value. The relevance of this systematic error calls for considering it within the fitted model. To this end, we use a “toy” Monte Carlo method based approach. Its design is based on the rationale found in different references (Rocha and Nogueira 2012; Clarkson 2014; Mahmoud and Hegazy 2017). The method can be summarized in four steps.

The first step is to define a bound for each measurement we have of the B_t, based on the identified measurement error. The second step is to take random values of errors (according to a Gaussian generation) at each measurement we have of B_t within the defined bounds. The third step is to fit the model according to these random values. Then, we repeat for a sufficient number of times (> 1000), where at each simulation we get different values for the regression parameters, and eventually determine their mean and standard error. Finally, the fourth step is to use the formula (Eq. 23) of error propagation (Ku 1966) to calculate the propagation from each regression parameter onto the model.

$$ E_{f}=\sqrt{ \left( \frac{\partial f}{\partial x}\right)^{2} {s_{x}^{2}} + \left( \frac{\partial f}{\partial y} \right)^{2} {s_{y}^{2}} + \left( \frac{\partial f}{\partial z} \right)^{2} {s_{z}^{2}} + \cdots} $$

(23)

where f is a function, which in our case, is the model we are interested in.

2.7 Consistency of our model predictions with proportions of bark on over-bark trunk volume

In order to check the relevance of our modelling results, we translated results from the IGN dataset into bark proportions (B_p, calculated with Eq. 24) using bark volume estimated by our selected models and over-bark stem volume estimated following (Tran-Ha et al. 2007) (Eq. 25). We translated the results from Emerge dataset into B_p using either the bark volume and over-bark trunk volume calculated according to the procedure found in Section 2.2 and Fig. 2 (Method 1), or the bark volume estimated by our selected models with BTBH measured and over-bark stem volume calculated as previously (Method 2). We then compared the results obtained according to these procedures with the results published in (FCBA 2019) and adapted them to our dataset (Emerge or IGN) according to the repartition of tree’s DBH and obtained by Meyer (1946) where, according to Eq. 10, B_p corresponds to (1 − a²). To calculate the B_p, (FCBA 2019) used the Emerge dataset to build a relation between the ratio of bark thickness to log radius at 1.3 m and B_p on the tree, then applied it to the IGN dataset. These relations were built for different DBH classes for each species.

$$ B_{p}=\frac{B_{v}}{V_{ob}} $$

(24)

$$ V_{ob}=a\cdot{}DBH+b\cdot{}\pi \left (\frac{DBH}{2} \right )^{2}\cdot{}H_{tot} $$

(25)

The a and b coefficients were estimated by Tran-Ha et al. (2007) for several species including our six targeted species.

For all these comparisons, we used Eq. 26 in order to calculate a mean B_p of all our trees.

$$ B_{p}=\frac{\sum{\frac{B_{v}}{V_{ob}}}}{n} $$

(26)

where n is the number of observations.

3 Results

3.1 BTBH modelling

The detailed plots showing data, trend lines and models for each species are given in Appendix Figs. 6 and 7.

3.1.1 Fit of literature models

Table 4 summarises the results obtained for models following Eqs. 4, 5, 6, 7 and 8.

Table 4 Attributes of the different models of bark thickness at breast height (BTBH). Model following Eq. 8 is shown when altitude effect is significant. The statistical significance is indicated by: NS: p ≥ 0.1; *:p < 0.05; **: p < 0.01; ***: p < 0.001. R² is the coefficient of determination. RMSErel is the relative Root Mean Square Error. AIC is the Akaike Information Criterion. PME is the propagation of measurement error for Norway spruce only. The models in bold correspond to the model finally retained in this study. It must be remembered that the different parameters a to d are model-dependent. Standard errors are given in Appendix Table 9

With regard to the model proposed by Cao and Pepper (1986) for the sessile oak species, d parameter was found insignificant. The model was thus recalibrated without this parameter. This new model is the one presented here. In terms of R², RMSErel and AIC, all models presented in Table 4 are quite similar. It can be observed that except for pedunculate oak and Norway spruce, the model proposed by Cao and Pepper (1986) is slightly better than others. For pedunculate oak, the three models are really equivalent. The PME calculation, which is possible in the case of the Norway spruce species, showed a 5.12%, 4.49%, and 5.67% error propagation in the case of Cao and Pepper (1986), Wilhelmsson et al. (2002), and Gordon (1983), respectively. For models following Eqs. 7 and 8, PME were 0.84% and 0.19%.

Models following Eq. 7 appeared to be rather good, with R² of 0.65 for European beech, 0.67 for sessile oak, 0.68 for pedunculate oak, 0.69 for Norway spruce, 0.71 for silver fir, and 0.77 for Douglas fir. However, the RMSErel around 30% (up to 37.3% for European beech) shows that the natural variation of B_t cannot be perfectly described by a simple model. For all species, the model according to Eq. 7 was thus not better than the best models found in the literature except for Douglas fir, pedunculate oak and European beech. The new models had an accuracy close to that of the literature models but had the advantage of having only two parameters and better results in terms of propagation of measurement error for the Norway spruce species.

After adjusting the model on altitude classes, we observed a significant relation between a and altitude only for silver fir, Norway spruce and European beech. This relations are given in Appendix Fig. 8. Furthermore, it can be seen that models following Eq. 8 are better than models of the literature. At the end we selected the model following Eq. 8 for silver fir, Norway spruce, and European beech, models following Eq. 7 for Douglas fir and pedunculate oak, and model following Eq. 5 for sessile oak (Table 4). Figure 3 shows the form of the selected models predicting BTBH for each species.

Fig. 3

Prediction of bark thickness at breast height (BTBH) from diameter at breast height (DBH) according to models selected in Table 4. a Models depending on altitude for silver fir, Norway spruce and European beech. b Models without altitude for Douglas fir, pedunculate oak and sessile oak. The dotted lines show the extrapolated parts of the BTBH

3.2 Models of B _v

Table 5 shows the modelling results for all species.

Table 5 Models of bark volume (B_v) with bark thickness at breast height (BTBH) and stem volume over-bark (V_ob) for Eq. 10 (Meyer 1946). The statistical significance is indicated by: NS: p ≥ 0.1; *: p < 0.05; **: p < 0.01; ***: p < 0.001. R² is the coefficient of determination. RMSErel is the relative Root Mean Square Error. AIC is the Akaike Information Criterion. PME is the propagation of measurement error for Norway spruce only. BP-test is the p-value of the Breusch-Pagan test without the structuration of the variance and with the structuration of the variance. The models in bold correspond to the models finally retained in this study. Standard errors are given in Appendix Table 10

First, by considering Breusch-Pagan test (BP-test) results, we observed that weights argument of gnls function handled correctly heteroscedasticity for all species except for the European beech models.

It can be observed that there is only a small difference between model following Eq. 17 and model following Eq. 18, in terms of AIC, RMSErel and R². Consequently, the approximation we made by assuming, DBH >> BTBH and H_tot >> 1.3 seems to be proper.

It can also be observed that results obtained with Eq. 18 are much better than the ones obtained by Meyers’ method (Eq. 10). RMSErel is decreased by 13% for silver fir, 12% for Norway spruce, over 9% for Douglas fir, 5% for sessile oak, 1.5% for pedunculate oak and 15% for European beech.

By comparing models following Eqs. 18 and 11 (Kozak and Yang 1981), it can be observed that, even if in terms of R² both models seem rather good, in terms of relative RMSE and AIC, the model proposed by Kozak and Yang (1981) is significantly better. We select this model even if the one following Eq. 18 is simpler.

The PME is, in the case of Norway spruce, approximately 24.5% for model following Eq. 18, while 65.73% in the case of Kozak and Yang (1981). As Meyer’s is not a model, PME can not be calculated.

Table 6 presents the results obtained with models following Eqs. 19 and 11 with parameters summarised in Table 5 and BTBH estimated with selected models presented in Table 4.

Table 6 Models predicting bark volume (B_v) with an estimate of bark thickness at breast height ($\widehat {BTBH}$) provided by the models in bold of Table 4. Coefficient estimates are provided in Tables 4 and 5. The models in bold correspond to the models finally retained in this study

By comparing models following Eq. 19, in Tables 6 and 5, several differences can be observed. For silver fir, Douglas fir, and European beech, a large increase of RMSErel can be detected when instead of using the model with BTBH (Eq. 18) we use the model with $\widehat {BTBH}$. For Norway spruce, pedunculate oak, and sessile oak, a small increase of RMSErel can also be observed between the two models respectively of 6%, 4%, and 6%.

The same accuracy loss can be observed for the model proposed by Kozak and Yang (1981). It can be nevertheless pointed out that for sessile oak, model following Eq. 19 became better than models proposed by Kozak and Yang (1981).

3.3 Bark proportion

In Table 7, column Emerge Method 1 represents the measured B_p and column Emerge Method 2 the B_p predicted with the bark volume model of Eq. 11 (Kozak and Yang 1981, best fit model see Table 5). Column (FCBA 2019) adapted for IGN, (FCBA 2019) adapted for Emerge and Meyers’Coefficient, (1 − a²) are added to simplified the comparison with the results and will be discuss in the discussion section.

Table 7 Bark proportions expressed in percent calculated for IGN and Emerge datasets with methods 1 and 2 (see Section 2.7), published in (FCBA 2019) and quantified by the Meyer’s coefficient (1 − a²)

Between methods 1 and 2, the only change is the prediction of bark volume. The prediction over-estimates the bark proportion for all species.

3.4 Application of BTBH and B _v models to NFI data to estimate bark resources at regional scale

Applying the best BTBH and B_v models to NFI data, as well as specific bark densities (Billard et al. 2020) and extractive content (Brennan et al. 2020), we were able to provide estimates of the volume, biomass and amount of extractives in the bark of all trees cut annually in Grand Est and Bourgogne-Franche-Comté regions. We will consider the whole stem, stump and tree top included. In total for the two regions, about 1.2 million m³ of bark were cut each year between 2009 and 2018. This represents a total bark biomass of 600 000 t/year and about 160 000 t/year of bark extractives. For each metric (volume, biomass, amount of extractives), the annual cutting was slightly higher in Bourgogne-Franche-Comté than in Grand Est (52% vs 48% of the total for the two regions, respectively). In addition, they were slight differences in harvest profiles in terms of species composition between the two regions. In Bourgogne-Franche-Comté, the annual cutting was dominated by pedunculate oak (22% of the regional bark volume cut each year) followed by Norway spruce (21%), Douglas fir and silver fir (18% each). By contrast in Grand Est, the annual cutting was dominated by Norway spruce (27% of the regional bark volume cut each year) followed by silver fir (20%), pedunculate oak (18%) and European beech (15%).

Table 8 and Fig. 4 show the ressources calculated for the six species in the two regions, Grand-Est and Bourgogne-Franche-Comté.

Fig. 4

a Repartition of bark volume removals between species. b Removals depending on species for two regions. Error bars represent the 95% confidence interval around means

Table 8 Summary of bark resources calculated for each species in Bourgogne-Franche-Comté (BFC) and Grand-Est (GE) regions. Values correspond to ± 95% confidence interval around mean

4 Discussion

4.1 Accuracy of measurements

Since a Swedish bark gauge has a high degree of inaccuracy (Stängle et al. 2016; Althen 1964), it is noteworthy to identify the different sources of errors that are propagated to the target models. To achieve this, we have used a method based on a Monte Carlo technique (see Section 2.6). At this stage, it is appropriate to mention that the amount of error propagated onto the model does not speak for the goodness of the model, but rather determines its accuracy and prediction power. Determining whether the model is good or not would require an instrumental variable method in the case of non-linear models (Pearl 2009), or a dedicated Monte Carlo simulations with pseudo-data and a detailed study of the bias and variance that could be generated on model parameters; this is outside the scope of this paper. In this article, we have applied this method on two kinds of selected models, BTBH and B_v.

For the BTBH model, we have seen that the error propagation (when altitude is unknown) for Eq. 5 is fairly small (5.12%). However, when altitude is known (Eq. 8) the error propagation is relatively smaller (0.19%). For the B_v models, Kozak and Yang’s 1981 model incurred a 64.8% error propagation. This can only be explained by the number of B_t measurements that are used to calculate B_v in the first place (several truncated cones with different B_t measurements), each measurement possessing an error on its own, eventually affecting its propagation within the model. It is expected that a better measurement technique could induce a higher inherent model accuracy.

4.2 Model application

In this work we applied models found in the literature to our targeted species. The B_v models following Eqs. 9 and 11 were initially designed for softwood and hardwood species. The differences in stem shape between both types of species were not taken into account. However, the models seem generic enough to fit both types of species in an equivalent way.

It can nevertheless be pointed out that for hardwood species, total bark volume may have been seriously underestimated as we do not consider the bark of branches. Ver Planck and MacFarlane (2014) estimated that, on average, branches of hardwoods account for 41% of the total wood volume in a tree.

4.3 Comparison of BTBH models

In spite of their variable forms, none of the tested models in this study provided RMSErel lower than 20%. This may be due to the natural variability of the bark, the measurement error described above, and/or the fact that no independant variables complementary to DBH and altitude were identified. Nevertheless, an effect of altitude was evidenced. The models that we recommend to use are shown in bold in Table 4.

4.4 Altitude effect

An effect of altitude was observed for three species (silver fir, Norway spruce, European beech). This effect is positive: the higher the altitude, the larger the parameter a and the more B_t increases with DBH. One of the reasons could be a different allocation of biomass between bark and the other parts of the tree, with the increasing altitude. No effect was observed for the other species (Douglas fir, sessile oak and pedunculate oak).

The difference between species can be first explained by a statistical reason: For some species data are available over a wide range of altitudes, while for others this range is more restricted. The altitude effect can be tested only in the first case. This hypothesis is confirmed because it turns out that for species for which the altitude effect is significant, i.e. silver fir, Norway spruce and European beech, the altitude ranges are respectively 0–1800 m, 0–2000 m and 0–1600 m. On the other hand, the range of altitudes is much more restricted for the other three species, respectively 0–1400 m, 0–1100 m and 0–1000 m for Douglas fir, sessile oak and pedunculate oak. To be exact it should be verified if in all cases, and particularly for Douglas fir, sessile oak, and pedunculate oak, this altitude range reflects the real distribution of trees elevation or represents only a part of it.

With regard to the effect of environmental factors, many studies have showed relationships between B_t and fire (Pausas 2015; Schafer et al. 2015; Bauer et al. 2010). Although fire is not a common issue for these species in France, bark thermal insulation properties are worth considering. Indeed, De Antonio et al. (2020) showed that bark properties, especially low bark density, protected buds against frost while Molina et al. (2016) showed a link between resistance to frost and B_t. They have worked on species of Brazilian savannah and in Patagonia. It can be assumed that this protection strategy also applied to our species. Further investigations may be required to analyse the actual effect of frost on increasing B_t. Moreover, it is questionable whether the protection strategy of Douglas fir, sessile oak, and pedunculate oak against frost may be based more on decreasing bark density than on increasing B_t.

Another factor influencing relative bark thickness is growth rate (Stängle et al. 2017; Stängle and Dormann 2018; Laasasenaho et al. 2005). Growth rate has a negative influence on relative bark thickness, meaning that the slower a tree grows, the bigger B_t will be compared to DBH. The fact that trees grow slower at higher altitudes, given the more difficult environmental conditions, may explain the observed relationship. However, we are not able to separate the influence of growth rate and the influence of altitude since the age of the trees was not recorded.

Rosell et al. (2014) showed that bark also has a function of water storage. Water storage mainly increases with B_t. As studied by Antoni et al. (2011), French mountains have a smaller usable reserve of water than lowland areas which could influence B_t. However, usable reserve is not only dependent on altitude, and thus further research is needed to validate this hypothesis.

If one of these reasons can be validated, it is reasonable to think that using these values (temperature or usable reserve of water) will enable a better model to be built, that would be better able to represent the actual effect that is linked to the variation of B_t.

Other variables such as latitude, longitude, and type of forest (coppice or high forest), and H_tot were tested but no significant influence was found.

4.5 Comparison of B _v models

In the Results section, we have selected the model following Eq. 11 taking into account RMSErel and R² values. However, it can be reasonably argued that the model following Eq. 18 is preferable since it has fewer parameters and thus is more robust and simpler to use. Moreover, the error associated to the parameters is smaller for the model following Eq. 18 than for the one following Eq. 11 (Appendix Table 10). Thus the model following Eq. 18, although slightly less accurate than model following Eq. 11, could be a good alternative.

Comparing Tables 5 and 6 it appears that, for silver fir, Douglas fir, and European beech, the accuracy losses by replacing BTBH by $\widehat {BTBH}$ are respectively of 23%, 14%, and 25%. Although BTBH measurements are found to be helpful when predicting B_v for these species, the models predicting BTBH for these species are not enough accurate to predict correctly B_v.

4.6 Comparison of B _p values predicted by the models

As shown by Table 7, compared to Emerge Method 1, the values predicted using the Meyer’s coefficients all over-estimated B_p except for European beech.

The difference with FCBA (2019) may be due to the method applied to the provided data. Indeed, FCBA (2019) provided B_p for different ranges of DBH starting from 25 cm. However, for computing the mean proportion we averaged these values weighed by the number of trees in each class and we applied the value of the lower class (25–30 cm) to all trees smaller than 25 cm.

The IGN dataset (column IGN) is closer to (FCBA 2019) (column (FCBA 2019) adapted for IGN) except for sessile oak, pedunculate oak, and Norway spruce. The high proportion of trees with a DBH smaller than 20cm can explain this difference, especially for Norway spruce. One may wonder if, for pedunculate oak and sessile oak, there is an important variation of B_p with respect to DBH for small trees. It can also be observed that the assumption made by Meyer (1946) of a constant B_p along the stem is rightful for silver fir, Norway spruce, and European beech. Indeed, the Meyers’ coefficient is close to the bark ratio measured for our trees except for Douglas fir and the two oak species.

5 Conclusion

To assess regional bark availability in terms of volume, biomass and quantities of extractive, we built several models to predict bark thickness at breast height and tree bark volume, from usual tree measurements such as total height and diameter at breast height and a site variable, altitude. This modelling was achieved for six temperate species: Silver fir, Norway spruce, Douglas fir, sessile oak, pedunculate oak, and European beech. We observed a statistical influence of altitude on bark thickness at breast height for three species, silver fir, Norway spruce, and European beech which opens the door to more ecological studies on bark. These models have been fitted on a pool of data collected in France. Considering the number of trees studied, the diversity of measurements made, and the number of bark thickness measurements made, these data are particularly rich and unique. The model set includes models developed elsewhere for other species or new ones created in this study. In this paper we have been able to valorise national forest inventory data, newly collected basic density and extractive rate data.

Appendix

Fig. 5

Diameter at breast height (DBH) distribution of the trees from Emerge and IGN datasets, for silver fir, Norway spruce, Douglas fir, pedunculate oak, sessile oak and European beech

Fig. 6

Relationship between standardized residuals and fitted values for silver fir, Norway spruce, Douglas fir, pedunculate oak, sessile oak and European beech for the models selected in Table 4. The red lines correspond to the horizontal axes, the green lines give the trend of the residuals

Fig. 7

Relationship between standardized residuals and fitted values for silver fir, Norway spruce, Douglas fir, pedunculate oak, sessile oak and European beech for the models selected in Table 5. The red lines correspond to the horizontal axes, the green lines give the trend of the residuals

Fig. 8

Relationship between altitude and the value of a parameter in model following Eq. 7 for silver fir, Norway spruce and European beech. Red line corresponds to the regression done

Table 9 Standard error associated to the parameters estimated for the models of bark thickness at breast height (BTBH). The model following Eq. 8 is shown when altitude effect is significant. The parameter values are given in Table 4

Table 10 Standard error associated to the parameters estimated for the model of bark volume (B_v). The parameter values are given in Table 5