Methods overview
A growth and harvest simulator was developed for the purpose of this study (see “Methods”—Growth and harvest simulator). The initial conditions of the forest and the growth equations of the simulator were parameterized by making use of the French National Forest Inventory Data (hereafter called NFI) (Fig. 5) for the years 2008–2012 [50].
Each plot of the IFN dataset was assumed to be homogeneous in species and stand structure, be it high stand or coppice. Density, fertility and exploitability indices were estimated from the IFN data (see “Characterizing stands using the French National Inventory data” section), used as proxies for the current management state (see “Management approaches” section), and used to build the growth and harvest simulator based on empirical growth relationships derived for volume and diameter (see “Methods”—Growth and harvest simulator).
The simulator assumes that the wood market is supply-driven: stands are harvested when they are mature, irrespective of the demand and wood price. Such a simplification is justified by the relative inelasticity of wood supply [75] and obviates the need to include wood demand or wood price scenarios. The simulated harvest at the first model time-step, 2015, was calibrated to match the nationwide total harvest per species [50], by adjusting the thresholds used to define the different management approaches, i.e., harvest diameter and thinning ratio (see “Methods” Calibration). By combining the descriptions of management and growth, and assuming constant environmental conditions, the evolution of French forest was simulated for a set of diverse but realistic management scenarios (Fig. 5). Furthermore, the wood-use chain was reconstructed from wood industry data. Forest growth, wood use, lifetime, and substitution potential of all wood products were accounted for in the carbon balance of the set of management scenarios which represented different intensification approaches.
The growth and harvest simulator was ran for a Business as Usual (BaU) scenario for all of France between 2010 and 2115 with 2040–2115 results solely used to determine the carbon parity time (see “Discussion”—Carbon parity time of intensified management scenarios). Subsequently, simulated annual harvests were fed into the French wood-use chain to project the stock of wood products as well as the energy production and emissions from wood use.
In this study, BaU represents a scenario in which no political action is taken, and no economic incentive exists, to increase wood harvest above present-day levels and the transformation industry is able to adapt its capacity to absorb the variation in wood supply. The BaU scenario thus assumed that no change in practice occurs, either in harvesting levels or in wood use. In the future, forests are thus assumed to be managed as they are today (see “Discussion”—Different management, different policy) and the harvested wood is assumed to be transformed in the same proportions into timber, pulp, and paper and energy with the same transformation efficiencies as today—whatever the simulated wood supply (see “Methods”—Wood use modelling). Changes to this assumption would require socioeconomic hypotheses that would impair the understanding of the processes driving the carbon balance and fall beyond the scope of this study.
Characterizing stands using the French National Inventory data
Every year a tenth of a 1 × 2 km grid over France is visited by IFN, resulting in the inventory of around 6500 to 7000 forest plots a year. On each inventory point, four circular plots of decreasing radius are established to measure and observe different variables with an appropriate level of detail [49]. Each stand is described in terms of the physical properties of the terrain, e.g., topography and soil properties, and its vegetation characteristics. Mainly tree-level and stand-level measurements were used in this study.
This study explicitly analysed the data of the 14 tree species which comprise 85% of the total area of forest in France (Table 1). Inventory plots for which the main species was not one of these 14 selected species, were pooled in the categories ‘Other hardwood’ or ‘Other softwood’. Some further simplification was made when upscaling the IFN data from the plot to the regional or national level. The IFN methodology under-samples a few homogeneous regions [49]. When upscaling the data, the effect of this sampling bias needs to be corrected by weighting plots by a statistical weight that describes their representativeness in terms of forest area. These statistical weights, however, are not publicly available [76]. Hence, an equal weight for all plots was used when averaging the variables within each species-region combination. The NFI-reported areas of each species-region combination is then used to upscale the data to the national scale.
Site index
Site index is an integrated metric commonly used in forestry to evaluate the quality of a site for tree growth [77]. It usually takes the form of a height at a given reference age, and is usually set to 100 years, as is the case in this study. From the IFN measurements of tree-level age and height, we reconstructed the stand-level variables of average stand age (average of the measured ages) and dominant height (h0; average height of the 100 trees with largest diameters) for each species. A Hossfeld II type equation [78]—a sigmoidal function commonly used to model biological growth—was fitted to the dominant height-age data, and used as a reference to describe the dominant height (h0; m) as a function of age [79]:
$$\begin{array}{*{20}c} {Guidecurve{:}\, h_{0\,guide} \left( {age} \right) = \frac{{a_{guide} }}{{1 + \left( {\frac{b}{age}} \right)^{c} }} = a_{guide} *f\left( {age} \right)} \\ {h_{0\,guide} \left( {100} \right) = H100_{guide} = a_{guide} *f\left( {100} \right)} \\ {\frac{{H100_{guide} }}{{h_{0\,guide} \left( {age} \right)}} = \frac{{f\left( {100} \right)}}{{f\left( {age} \right)}}} \\ \end{array}$$
(1)
where a, b, and c are fitted coefficients, h0 (m) and age (years) are the plot dominant height and age, respectively (Additional file 6: Table S3). Assuming a constant shape across fertility classes (constant parameters b and c, hence same function f) the ratio between the site index and the height at age t (years) are constant across fertility indices [80]:
$$\frac{{H100_{plot} }}{{h_{0\,plot} \left( {age} \right)}} = \frac{{H100_{guide} }}{{h_{0\,guide} \left( {age} \right)}}$$
(2)
where H100plot (m) is the dominant height at an age of 100 years sought for a given plot, h0plot(age) (m) is the observed height for the current age of the plot, H100guide (m) is the dominant height calculated from the guide curve for an age of 100 years, and h0guide(age) (m) is the dominant height calculated from the reference curve for the age of the considered plot. H100 was used as the site index and the fertility indices for each species were aggregated into fertility classes for which the limits were defined by the observed quantiles of the distribution of site indices for each species (Additional file 7: Figure S5).
Density index
The relative density index (DI) is commonly used to quantify management intensity [81]. It has been defined as the ratio of the number of stems in a stand and the expected number of stems for a stand of the same mean diameter [82] under self-thinning [83]. For each species, the boundary of the log–log relationship between the observed quadratic mean diameter and the stand’s number of stems for stands with age between the 10th and 90th percentiles was fitted to quantify the self-thinning relationship [84]. Owing to the presence of unmanaged stands in French forests, some of the observations can be expected to be close to self-thinning, which is required for this approach to be valid.
The DI is then calculated as:
$$DI = \frac{N}{{N_{max} }} = \frac{N}{{e^{g} Dg^{h} }}$$
(3)
where DI is the density index (unitless) for a stand characterized by its number of stems N (trees ha−1) and quadratic mean diameter Dg (cm) as calculated from the inventory’s circumference data. Nmax (stems ha−1) is the maximum number of trees with quadratic mean diameter Dg according to the fitted self-thinning relationship. Parameters g and h are fitted to define the self-thinning function as log(N) = g + h log(Dg). An example of a fitted species-specific self-thinning relationship is shown in Additional file 8: Figure S6 and parameters of the self-thinning relationships are shown in Additional file 6: Table S3.
Exploitability index
Following the methodology established by IGN [49], plot-level records of terrain slope, ruggedness, distance to forest ride, carrying capacity of the ground for machinery were used to define four levels of exploitability: very easy, easy, difficult, impossible. The details of the categories are shown in Additional file 9: Table S1.
Management approaches
Plots were classified as one of four management approaches based on their physical and biological characteristics: quadratic mean diameter, fertility index, relative density index, and exploitability index. The use of the exploitability index as a physical constraint builds on the assumption that the stands which were easier to harvest and thin, are more likely to have been harvested. However, depending on the wood price, targeting stands with a lower exploitability index could be justifiable. Price dependency may have contributed to the observation that P. pinaster stands were managed irrespective of their exploitability indices. Harvest and thinning of all other species started with the exploitability index set as ‘difficult to access’.
Discontinuity in the statistical distribution of the diameter was used to determine the expected clearcut diameter for each species and each fertility class. The discontinuity was identified using a segmented regression approach. The rationale for this is the assumption that for a given species, under idealized management and an even age-class distribution, all diameter classes below the harvest diameter would have a similar number of stands whereas no stand would have a mean diameter above the harvest diameter. Subsequently, the clearcut diameters calculated with this method served as a prior and were adjusted during calibration of the harvested wood volume per species, as reported by IFN for the year 2010 [50]. Also, for each species, a density threshold was set based on DI distributions with a segmented regression approach similar to the one used for the diameters. Subsequently, the final density threshold was determined through calibration of the harvested wood volume per species, as reported by IFN for the year 2010 [50].
Each plot is assigned to only one of the four management approaches by following an assignment order. First, plots with too low an exploitability index were classified as unexploitable. Then, the remaining plots with diameters above the calibrated expected quadratic mean diameter were referred to as ‘harvest-delayed’ and assumed not to be harvested for ownership reasons, for example, fragmentation of properties or disinterest of the owner or for conservation. Under present management, or Business as Usual, these plots are neither thinned nor harvested. Then, plots with densities above the final DI threshold were referred to as overstocked. In a Business as Usual scenario, these plots will not be thinned but will be harvested. Finally, the remaining forest with diameters, density and exploitability not matching any of the above criteria, were classified as actively managed. For P. pinaster and Picea abies stands, the harvest statistics were best matched when no threshold was set either on stocking density or on exploitability index.
Growth and harvest simulator
Forest growth
The growth and harvest simulator calculates the evolution of each plot in the inventory, based on its species, density and fertility characteristics. The 5-year radius increment and the 5-year volume increment at the stand-level are the two key variables of the growth and harvest simulator. The incremental variables are modelled, rather than absolute variables, because the initial distributions of the volume and diameter are preserved. Volume and radius increment are independently simulated (Eqs. 4 and 5, respectively) as a function of stand age, its growing stock and stand fertility. The structure of the equations used in the growth and harvest simulator are commonly used in the forest modelling community [80].
$${ \log }\left( {IV_{stand} } \right) = a_{H100} + {\text{a}}_{DI} \cdot { \log }\left( {DI} \right) + {\text{a}}_{age} \cdot { \log }\left( {age} \right) + \varepsilon$$
(4)
$$log\left( {Ir_{stand} } \right) = b_{H100} + b_{DI} \cdot log\left( {DI} \right) + b_{age} \cdot log\left( {age} \right) + \varepsilon$$
(5)
where IV is the 5-year volume increment (m3 ha−1), age is the plot age (years), DI is the density index (unitless), ε is the residual error term, a and b the regression coefficients for the log-transformed variables [80].
Ultimately, Eqs. 4 and 5 were to be used to simulate the future productivity of French forests under different intensification scenarios and BaU. This, however, required that the parameters of this set of equations, i.e., aH100, aDI, aage, bH100, bDI, bage were fitted first on inventory-based plot observations. As the stand-level volume and above-bark diameter increments are not reported in the inventory, IVstand and Irstand were estimated from the plot measurements. The subsequent paragraphs of this section detail how IVstand and Irstand were calculated.
Derivation of stand-level volume and radius increment
The modelling at the stand level required making an assumption on stand composition. Each stand in the national inventory data was assumed to be homogeneous in terms of its species composition. The species for which the age was recorded in the NFI, was considered the main tree species of the stand and was used in the species-specific equations of the growth and harvest simulator.
The stand-level volume increment (IVstand; m3) is calculated as the weighted sum (wtree) of the volume increments of all trees (ivtree; m3) in the stand irrespective of their species.
$$IV_{stand} = \mathop \sum \limits_{stand} iv_{tree} \cdot w_{tree}$$
(6)
where, wtree (unitless) is the statistical weight of each tree within a plot as a result of the tree’s selection area and proximity to its border, and ivtree(m3 per 5-years) is the 5-year tree-level volume increment. The variable ivtree is not provided by the NFI and was, therefore, reconstructed from data on the tree volume (vtree; m3), tree radius (rtree; m) and tree radius increment (irtree; m per 5-years) all three reported by IGN, combined with bark ratios (β; unitless) from a French measurement campaign [85]. For this, relations (f) between reported tree volume (vtree; m3) and tree radius (rtree; cm) were fitted independently for each species (Eq. 8) yielding current and previous time-steps tree volume estimates (\(\tilde{v}_{tree\left( t \right)} ,\tilde{v}_{{tree\left( {t - 5} \right)}}\); m3) from reported current radius (rtree(t); m) and reconstructed previous time-step radius (\(\tilde{r}_{{tree\left( {t - 5} \right)}}\); m) all considered above-bark (Eqs. 8 and 9). Previous time-step above-bark radius was calculated from inventory-reported current radius (rtree(t); m), below-bark radius increment (\(\widetilde{ir}_{tree\left( t \right)}\); m per 5-years) and bark ratios (β) (Eq. 10).
$$iv_{tree} = v_{tree} \left( t \right)\left( {1 - \frac{{\tilde{v}_{{tree\left( {t - 5} \right)}} }}{{\tilde{v}_{tree\left( t \right)} }}} \right)$$
(7)
$$v_{tree\left( t \right)} = f\left( {r_{tree} } \right) + \varepsilon = \tilde{v}_{tree\left( t \right)} + \varepsilon$$
(8)
$$\tilde{v}_{{tree\left( {t - 5} \right)}} = f\left( {\tilde{r}_{{tree\left( {t - 5} \right)}} } \right)$$
(9)
$$r_{{tree\left( {t - 5} \right)}} = \frac{{\left( {1 - \beta } \right)r_{tree\left( t \right)} - ir_{tree\left( t \right)} }}{1 - \beta }$$
(10)
For the simulation of diameter–related variables (quadratic mean diameter, diameter growth and stand relative density), it was assumed that each stand could be simulated by the replication of an ‘idealized average tree’ of the main species. The characteristics of the idealized average tree of the stand were based on the observed characteristics of all the trees of this species in the stand. The stand-level radius increment (Irstand) was hence calculated as the average radius increment of all trees of the main species only (irtreeSp), as reported by the inventory, to account for diversity of growth strategies between different species in mixed stands.
$$Ir_{stand} = \frac{{\mathop \sum \nolimits_{stand} ir_{treeSp} \cdot w_{treeSp} }}{{N_{treeSp} }} = \frac{{\mathop \sum \nolimits_{stand} ir_{treeSp} \cdot w_{treeSp} }}{{\mathop \sum \nolimits_{stand} w_{treeSp} }}$$
(11)
where Irstand is the 5-year radius increment of the idealized average tree replicated in the stand (m), irtreeSp and wtreeSp are individual tree radius increment (m) and weighting coefficient (unitless), respectively, as given in the NFI data.
The quadratic mean diameter on which harvest decision is based is also calculated from the average basal area for trees of the main species with gtreeSp (m2) the individual basal area of trees of the main species and wtreeSp their statistical weights (unitless).
$$Dg_{avgtree} = 100\sqrt {\frac{4}{\pi }\frac{{G_{Sp} }}{{N_{Sp} }}} = 100\sqrt {\frac{4}{\pi }\frac{{\mathop \sum \nolimits_{stand} g_{treeSp} }}{{\mathop \sum \nolimits_{Stand} w_{treeSp} }}}$$
(12)
where the subscript Sp refers to the subset of the stand trees being of the main species. Dgavgtree is the quadratic mean diameter of the idealized average tree of the stand, GSp and NSp are the basal area (m2) and number of stems of the trees of the stand’s main species (unitless).
The idealized average tree is then replicated to make up an idealized stand for which the volume is identical to the volume reported in the NFI, providing a reconstructed number of stems in the stand Navgtree used for the calculation of the relative density index.
$$v_{avgtree} = \frac{{V_{sp} }}{{N_{sp} }} = \frac{{\mathop \sum \nolimits_{stand} v_{treeSp} \cdot w_{tree} sp}}{{\mathop \sum \nolimits_{stand} w_{treeSp} }}$$
(13)
$$N_{avgtree} = \frac{{V_{stand} }}{{v_{avgtree} }}$$
(14)
where vavgtree is the volume of the idealized average tree, VSp and NSp are the total volume (m3) and number of stems of the stand’s main species (unitless), respectively, and vtreeSp and wtreeSp are individual tree volume (m3) and statistical weighting coefficient (unitless), respectively. The latter variable describes the representativeness of the tree within the plot, as given in the NFI data.
Mortality
No relationship was found in the national inventory data between mortality on the one hand, and age, standing volume and density index on the other hand. Hence, species-specific average mortality rates derived from the national inventory data were applied to all plots as a yearly percentage of standing volume irrespective of plot age. The species-specific mortality rate was defined as the average dead to living volume recorded in the inventory data. Hence, the growth of a plot is the result of the balance between 5-year volume increment and 5-year mortality. After a given age, if a plot’s mortality volume is above its volume increment, then the stand’s volume decreases accordingly. Dead organic matter and soil carbon pools are not included in the model due to the large uncertainties on the effect of management on these carbon stocks and expected minimal differences they would bring between our scenarios.
Stand replacement
According to the inventory methodology, trees with diameters below 7.5 cm in diameter are not inventoried [49]. The stands where such trees dominate have not been included for deriving the equations of the growth and harvest simulator, which are hence not valid for very young stands. The first 5 years of post-harvest growth are therefore not simulated by the growth simulator. When stands of more than 200 trees were classified, and trees aged less than 10-years old were reported, but identified as being below the inventory threshold, their per species characteristics were averaged and assumed for the first post-harvest time-step. Data for stands classified as below the inventory threshold and containing less than 200 trees were assumed representative of the recently harvested stands rather than the regrowth and hence were not considered [49]. For species where no data was available according to this constraint, post-harvest stands were replaced by a hypothetical regrowth stand with characteristics of the across species average diameter and volume. Following a harvest, the density and fertility indices of the newly established stand were considered constant as proxies of fertility and management.
Thinning and clearcutting
For each species a thinning ratio was defined based on yield tables [86]. The standing volume given by this ratio was subtracted each year from all plots which were subject to thinning, i.e., the actively managed stands under BaU. Thinning affects only stand volume because in our growth and harvest simulator the diameter is not updated. A diameter criterion determined by a species-specific segmented regression was used to decide the final harvest. In the scenarios mobilizing biomass from harvest-delayed stands, the mobilization objective was gradually implemented over the length of the simulation (30 years) following a linear trend.
Calibration
The harvested volume of a given species largely depends on the distribution of its stands across the four management approaches and of the yearly thinning ratio. Hence, an initial species-specific mismatch between observed and simulated harvest volumes could be decreased by calibrating the prior thresholds defining the different management approaches. Management approaches are defined by thresholds for exploitability, density index, and harvest diameter.
The threshold on exploitability index was the first to be calibrated as it makes the strongest and most direct impact on harvested volume. Indeed a stand in which exploitability is below the exploitability threshold is simply not managed. As a prior, all stands with exploitability index ‘impossible’ were classified as unexploitable. After iterative comparison of simulated and reported harvest levels, the high harvest levels of P. pinaster lead to none of these stands being assigned to the unexploitable category. This result is probably because, the well-developed marketing opportunities for P. pinaster incentivize harvest, even when the terrain is sub-optimal.
The two thinning-related parameters, i.e., threshold DI and thinning ratio, have similar effects (Additional file 10: Table S2). The threshold DI for defining the overstocked stands can be raised to decrease the number of overstocked stands, hence increasing the total harvested biomass. Yet, the fraction of biomass removed from thinning at each time-step can also be increased to increase the harvest. Because the prior values of the thinning ratios were based on fragmented data, we believe they come with larger uncertainties than the priors for density index. The priors for thinning ratio were therefore used for optimization against the reported harvest data. This procedure led to Q. petraea, P. abies and P. pinaster thinning ratios being increased from 5, 5 and 10% to 8, 11 and 13%, respectively.
The harvest diameter, above which stands are harvested, was not calibrated (Additional file 10: Table S2). The prior value was thus used for all species. Indeed, this parameter was found to have a rather limited but complex impact on harvested volume because it acted on two different processes in the growth and harvest simulator. First, by setting the threshold for harvest-delayed plots, the harvest diameter controls whether a given plot will be harvested or not. Second, it also defines when all plots under management will be harvested. Hence, depending on the actual diameter distributions of the stands, an increase in the clearcut diameter can increase or decrease the total harvested volume for a given species.
Intensifying forest management
The objective of intensifying current forest management is to reach the national energy target for bioenergy production from woody biomass. The simulation experiments applied in this study intensified current forest management by considering three different approaches: decreasing the harvest diameter of actively managed plots, thinning currently overstocked plots and harvesting currently harvest-delayed plots. Each approach entails a family of 10 simulations, which represent the intensity of a given policy targeting: thinning 10 to 100% of the overstocked forest area (scenarios Ov1 to Ov10), harvesting 10 to 100% of the harvest-delayed forest area (scenarios D1 to D10), reducing the clearcut diameter of all actively managed forests by 1 to 10 cm (M1–M10).
When the scenario is set to mobilize only a fraction of the targeted forest, the mobilized plots are the ones with the largest diameters (for harvest-delayed forests) or with the highest density (for overstocked forests). Whatever the degree of intensification, the extra mobilization is reached progressively within 30 years by linearly increasing the number of plots targeted from zero to the target area.
Wood-use modelling
Wood use substantially contributes to the net carbon balance of French forests, because different wood uses come with different lifetimes, shares of wood reuse and potential substitution of fossil fuel emissions [87]. Therefore, wood flows from the forest to the final product including recycling and the products’ end-of-life were reconstructed. To simplify the scheme, recycling was assumed to occur only once. Data were extracted from different literature sources [16, 88,89,90,91,92,93,94,95,96] for the year 2010 and reconciled into a unified scheme (Fig. 6).
The wood product scheme starts by assigning harvest to either the lumber or the pulp and energy pathway. The proportions of soft and hardwood going down the timber, pulp and energy pathways were calibrated to match the reported national production of 18.6, 10.4 and 28.3 Mm3 including self-consumption energy-wood, respectively [90]. For hardwood, 16%, 14% and 70% of the harvest, were directed to the timber, pulp and industry, and energy pathways, respectively. For softwood, 53%, 23% and 24% of the harvest were directed to the timber, pulp and industry, and energy pathways, respectively. The assignment of harvest from different management practices to timber or pulp and energy was also adjusted to match the reported national values, assuming that harvest from high stand clearcut is preferentially used as timber, whereas the wood from thinning and coppicing is preferentially used in the production of pulp and energy. For softwood, all clearcut and 26% of thinnings were directed to timber, and 74% of thinnings were directed to pulp and paper. For hardwood, all thinnings and 43% of clearcuts were directed down the pulp and energy pathway, and 57% of the clearcuts were directed to timber.
Carbon balance
The carbon balance is calculated by summing all components of the life cycle of wood: sequestration in the forest, carbon storage in wood products, carbon emission at the end of the product’s life and the avoided emissions from substitution. Each wood use is assigned an average life expectancy and substitution potential (Table 4). The BaU scenario is used as the reference for substitution calculation hence assuming that literature-based substitution coefficients used were derived for reference and scenario conditions consistent with the wood harvest and wood-use practices implicitly described in our scenarios. The avoided emissions under each scenario were calculated by multiplying the extra wood products with respect to BaU by their respective substitution coefficients.
Carbon storage in wood products is only accounted for in the first wood use (timber end-product, pulp and paper end-product or energy) and the carbon contained in the products is emitted progressively at each time-step following an exponential decay function. However, for substitution effects, consecutive wood uses are accounted for by multiplying their respective volumes with product-specific substitution coefficients all supposed to occur at time of harvest. This assumption is considered to have a minor effect on the results since most consecutive wood uses occur in a time lapse below the 5-year time-step of the forest growth and harvest simulator. This time compression approach leads to a small overestimation of substitution potential of harvest directed to timber with the anticipation of the substitution benefits from product recycling.