Management Strategies for Wood Fuel Harvesting—Trade-Offs with Biodiversity and Forest Ecosystem Services

Bioenergy is expected to contribute to mitigating climate change. One major source for bioenergy is woody biomass from forests, including logging residues, stumps, and whole trees from young dense stands. However, at increased extraction rates of woody biomass, the forest ecosystem, its biodiversity, and its ability to contribute to fundamental ecosystem services will be affected. We used simulation and optimization techniques to assess the impact of different management strategies on the supply of bioenergy and the trade-offs between wood fuel harvesting, biodiversity, and three other ecosystem services—reindeer husbandry, carbon storage, and recreation. The projections covered 100 years and a forest area of 3 million ha in northern Sweden. We found that the development of novel and cost-effective management systems for biomass outtake from young dense stands may provide options for a significant supply of bioenergy to the emerging bioeconomy, while at the same time securing biodiversity and important ecosystem values in future stand developments. In addition, there is potential to increase the extraction of harvest residues and stumps while simultaneously improving conditions for biodiversity and the amount of carbon stored in forest ecosystems compared to current levels. However, the projected continuing trend of increased forest density (in terms of basal area) has a negative impact on the potential for reindeer husbandry and recreation, which calls for researching new management strategies on landscape levels.


Introduction
There are high expectations that bioenergy will contribute to mitigating climate change [1,2]. Today, woody biofuels account for about half of the EU's total renewable energy consumption, and the demand for bioenergy will most likely increase [3] and references therein. Sweden is one of the main wood producers within the EU and is considered to be an important source for providing forest residues to the European bioenergy market [3]. In addition, the Swedish government adopted a climate goal of net zero emissions of greenhouse gases by 2045, and negative emissions thereafter [4]. The forest is expected to contribute to reaching the goal, inter alia by supplying more bioenergy [5]. It is therefore anticipated that Sweden's woody biofuel demand will increase considerably, with an additional demand for forest-based feedstock in the emerging chemical-, petrochemical-, and energy industries of about 40-55 TWh in 2030 and 60-70 TWh in 2050 [6]. For comparison, the total energy supply in Sweden amounted to 565 TWh in 2017, of which 143 TWh consisted of biofuels [7]. Sweden and Finland are the two countries within Europe that utilize the largest volumes of woody biofuels, and as the demand for bioenergy increases it is expected that many other countries will follow suit [8].
(2) Identify a compromise solution for attaining the largest possible harvest levels of wood fuels combined with the least negative impact on other ecosystem services and biodiversity.
The analysis is done using the county of Västerbotten in northern Sweden, with more than 3 million ha productive forest1 (>10% of the total Swedish forest area; productive forests are forests with a mean annual increment > 1 m 3 /ha) as the case study area. Västerbotten county represents the forest conditions in northern Sweden, where only a small proportion of the available primary woody biofuels is utilized currently. We use an advanced forest decision support system and assess the trade-offs between biofuel extraction and other ecosystem services and biodiversity under two different biofuel extraction settings: extracting only logging residues from final fellings, and more intensive extraction practices additionally, including logging residues from thinnings, stumps from final fellings, and whole-tree extraction from dense young stands. All analysis is done with the assumption that the harvest of timber volume and pulpwood volume should be in line with the potential sustainable harvest levels based on projections done in the latest national scenario analysis [26].

Case Study Area
The current forest conditions in the case study area were represented by a total of 2738 national forest inventory (NFI) plots, measured between the years 2008 and 2012. The NFI is an annual survey of all forest land in Sweden, initiated in 1923, where circular plots with a radius of 7 or 10 m placed along borders of squared tracts distributed in a systematic cluster design are inventoried [27]. Each plot represents a forest area of up to 1800 ha. Of the total productive forest area in Västerbotten, 3.0% is located in nature reserves, 7.5% in voluntary set-asides, and 6.7% in retention patches (based on the latest national forest impact analysis [26]). The mean age of the forest is 68 years, with almost 40% of the forest being younger than 40 years. Forests in nature reserves are primarily older than 80 years ( Figure 1). The mean standing volume is 100 m 3 /ha, and the main species are Scots pine (Pinus sylvestris, 43%), Norway spruce (Picea abies, 38%), and birch (Betula spp., 16%). The remaining volume is divided almost equally between other broadleaves (1.4%) and lodgepole pine (Pinus contorta, 1.5%).
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 21 Figure 1. Location of the study area (a), volume distribution by age class (b), and age class distribution (area) of the productive forest in the study area (c).

Ecosystem Services
We followed the ecosystem service classification of the Millenium Ecosystem Assessment [28], and in addition to woody biofuel extraction we included one provisioning service (reindeer husbandry), one cultural service (recreation), and one regulating service (carbon sequestration). We also included biodiversity, which can be regarded as a cultural service or supporting service but also as a foundation for all other ecosystem services, i.e., as a separate quality [29] (Table 1).  Figure 1. Location of the study area (a), volume distribution by age class (b), and age class distribution (area) of the productive forest in the study area (c).

Ecosystem Services
We followed the ecosystem service classification of the Millenium Ecosystem Assessment [28], and in addition to woody biofuel extraction we included one provisioning service (reindeer husbandry), one cultural service (recreation), and one regulating service (carbon sequestration). We also included biodiversity, which can be regarded as a cultural service or supporting service but also as a foundation for all other ecosystem services, i.e., as a separate quality [29] (Table 1). Reindeer husbandry is an important industry for the Sami people in northern Sweden, which justifies a consequence analysis of increased levels of woody biofuel harvesting. As an indicator for reindeer husbandry, we used the forest area that has the potential to support the growth of ground lichen based on [30]. Ground lichen, which is abundant especially in old, open pine forest, is an important food source for reindeer during autumn and winter, and has declined substantially due to modern forest practices [31].
Recreational activities in forests, including walking, jogging, berry and mushroom picking, as well as hunting, are very common both in Sweden and other boreal countries, and the importance and use of forests for recreational activities has been frequently recognized [32,33]. As indicator for recreation we used an index between 0 and 1 describing a plot's suitability for recreation, based on a set of forest variables such as average height, trees/ha, and tree species distribution [34].
Carbon storage in forests contributes to the mitigation of climate change. In this study, carbon storage was calculated as the total forest carbon stock, including carbon in above-and belowground tree biomass and in the forest soil. Soil carbon storage was based on the Q-model, which is implemented in the Heureka system, and predicts the amount of carbon and nitrogen present in the litter and humus layer at a given time period [35].
We had no direct measure of biodiversity and therefore used proxies that are known to correlate with high species richness and/or the occurrence of rare species, and that are used as indicators to monitor the Swedish environmental quality objective Sustainable Forests [36]. In this study, we used the amount of dead wood, old forest, mature broad leaf forest, and number of large trees as indicators for biodiversity. Numerous studies support the fact that dead wood is one of the most important prerequisites for preserved biodiversity [37]. Dead wood in this study included all stemwood resulting from natural mortality with a minimum diameter of 10 cm, and in all decay classes. Old forest is important to some species relying on continuity [38]. Large trees are amongst the most important substrates for red-listed species in boreal forests [39]. During the last century, there has been a change in the landscape structure in that conifer species (pine and spruce) have become more dominant. This has led to a situation where biodiversity components tied to the presence of deciduous trees are under pressure [40].

Analytical Framework
The trade-offs between increasing the extraction of woody biofuels and the impacts on the indicators for biodiversity and ecosystem services were investigated through a two-step procedure using long-term simulation of future forest development under different management alternatives combined with optimization based on LP. To this end, we used the Heureka application PlanWise, which is a decision support system for analysis and planning of the forest landscape [41]. PlanWise uses an optimization approach to determine the best possible combination of management strategies to meet user-defined objectives and constraints. This is done in two steps: treatment generation, and treatment selection (optimization).
During treatment generation, a large number of treatment schedules are created for each plot, for one or several management strategies. Treatment schedules vary in the timing of silvicultural activities, while management strategies can vary in management system (unmanaged, even-aged, and uneven-aged) and the type of silvicultural activity applicable (for example, natural regeneration/planting, number of thinnings, minimum rotation length, or different logging residue extraction practices, or different regeneration measures). For each plot, treatment schedule and period, PlanWise simulates the development of the tree layer by using a large set of models, e.g., empirical growth and yield models, models for stand establishment, diameter and height growth, ingrowth, and mortality. These models were developed by means of regression analysis using data from the National Forest Inventory, long-term experiments, and yield plots [42][43][44]. Other components of the forest ecosystem, like deadwood and carbon, were then calculated based on the development of the tree layer.
In treatment selection, an optimization problem was solved using linear programming (LP) or mixed integer programming (MIP) aiming at finding the optimal combination of treatment schedules for each plot, based on a user-defined objective that was to be minimized or maximized, and constraints that were to be met.
The following sections describe the details of the treatment generation and selection process applied in this study. Figure 2 gives an overview of the analytical framework. to meet user-defined objectives and constraints. This is done in two steps: treatment generation, and treatment selection (optimization). During treatment generation, a large number of treatment schedules are created for each plot, for one or several management strategies. Treatment schedules vary in the timing of silvicultural activities, while management strategies can vary in management system (unmanaged, even-aged, and uneven-aged) and the type of silvicultural activity applicable (for example, natural regeneration/planting, number of thinnings, minimum rotation length, or different logging residue extraction practices, or different regeneration measures). For each plot, treatment schedule and period, PlanWise simulates the development of the tree layer by using a large set of models, e.g., empirical growth and yield models, models for stand establishment, diameter and height growth, ingrowth, and mortality. These models were developed by means of regression analysis using data from the National Forest Inventory, long-term experiments, and yield plots [42][43][44]. Other components of the forest ecosystem, like deadwood and carbon, were then calculated based on the development of the tree layer.
In treatment selection, an optimization problem was solved using linear programming (LP) or mixed integer programming (MIP) aiming at finding the optimal combination of treatment schedules for each plot, based on a user-defined objective that was to be minimized or maximized, and constraints that were to be met.
The following sections describe the details of the treatment generation and selection process applied in this study. Figure 2 gives an overview of the analytical framework.

Treatment Generation
For each NFI plot, a set of possible treatment schedules was generated covering the entire planning horizon, which in this study was set to 100 years, divided into 20 five-year periods. The treatment schedules for a given management strategy differed in the timing of forestry activities within the framework given by the management strategy. For each schedule, the outcome in terms of indicator values for the selected ecosystem services and biodiversity was simulated.
For the production forest (83% of the forest area, excluding nature reserves, voluntary set-asides, and retention patches), we defined six different possible management strategies: wood-production, no thinnings, broadleaves, long rotations, continuous cover forestry (CCF), and unmanaged ( Table  2). The wood-production oriented strategy represents standard management practices, comprising regeneration through planting, cleaning, several thinnings, followed by final felling with the retention of three high stumps and ten live trees per hectare. The no thinning strategy includes cleanings, but no commercial thinnings. The broadleaves strategy aims at increasing the proportion of broadleaves in the landscape by increasing the share of retained broadleaves in cleaning and

Treatment Generation
For each NFI plot, a set of possible treatment schedules was generated covering the entire planning horizon, which in this study was set to 100 years, divided into 20 five-year periods. The treatment schedules for a given management strategy differed in the timing of forestry activities within the framework given by the management strategy. For each schedule, the outcome in terms of indicator values for the selected ecosystem services and biodiversity was simulated.
For the production forest (83% of the forest area, excluding nature reserves, voluntary set-asides, and retention patches), we defined six different possible management strategies: wood-production, no thinnings, broadleaves, long rotations, continuous cover forestry (CCF), and unmanaged ( Table 2). The wood-production oriented strategy represents standard management practices, comprising regeneration through planting, cleaning, several thinnings, followed by final felling with the retention of three high stumps and ten live trees per hectare. The no thinning strategy includes cleanings, but no commercial thinnings. The broadleaves strategy aims at increasing the proportion of broadleaves in the landscape by increasing the share of retained broadleaves in cleaning and thinning operations, and allowing for longer rotation periods. In the long rotations strategy, pine is regenerated naturally and rotation time is prolonged by up to 60 years after reaching the minimum legal age for final felling. In the CCF strategy, which was only applied in spruce-dominated stands, the forest is managed with a series of selection fellings. Finally, a set-aside strategy was simulated for all production forests, in which the forest was left unmanaged. Table 2. Description of the management strategies applied in the business-as-usual (BAU) and bioeconomcy (bioE) extraction settings for the productive forest. Harvesting of woody biofuel is an option in all strategies except for the set-aside strategy. To be able to compare the effect of different ambitions regarding woody biofuel harvesting, the management strategies were simulated under two different settings for biofuel extraction: business-as-usual (BAU), representing current practices, and bioeconomy (bioE), with intensified biomass extraction options. In BAU, the only primary biofuel sources are logging residues (tops and branches) from final fellings. In bioE, additional biofuel sources are small diameter trees from whole tree harvesting in early thinnings (termed "biofuel thinning" in the remainder of the text), logging residues from regular thinnings, and stumps (including large roots) from final fellings. In biofuel thinnings, all the above-ground tree biomass is used for biofuel. Due to economic consideration, logging residues and stumps are only extracted in plots with a minimum spruce proportion of 50%, in both biofuel extraction settings. We applied environmental restrictions to logging residue and stump extraction based on recommendations of the Swedish Forest Agency [26]. These included that no logging residue and stump extraction take place on wet soils or peatland, due to the risk for rutting caused by forest machine trafficking. On moist soils, stumps are not extracted at all, while logging residue extraction is allowed on sandy soils, but not on clay-or fine silt soils. Further, stumps are not extracted on lichen-rich plots due to considerations The remaining forest, 17% of the area, was set aside for nature conservation without any woody biofuel extraction, in both settings. Nature reserves were left unmanaged. In voluntary set-asides and retention patches, some nature-conservation focused management was implemented, by applying a management strategy consisting of selection felling on 10% (voluntary set-asides) or 20% (retention patches) of the area, leaving the remaining area unmanaged.

Treatment Selection
After simulating alternative treatment schedules for each NFI plot, a set of optimization problems were formulated and solved with LP for both biofuel extraction settings (BAU and bioE). All problems consisted of assigning treatment schedules to each NFI plot so that the value of the objective function was maximized while the defined restrictions were met. All optimizations included a constraint specifying the demand for roundwood harvest, based on projections done in the latest national scenario analysis on potential sustainable harvest levels [26]. Annual roundwood harvest demands range from 6.1 in the beginning, to 6.9 million m3 (under bark) in the end of the planning horizon ( Figure 3). The roundwood demand includes firewood and unmerchantable stemwood, as a proportion of regular roundwood harvests. We did not include these volumes in our assessment of potential woody biofuels. We applied the same harvest demand constraint in all optimizations, allowing at most a +/−1% deviation from the harvest demand in each period. In addition, a constraint for evenness in woody biofuel extraction was included in all problems by accepting a maximal deviation of 20% in each period from average woody biofuel extraction over the entire planning horizon, which is in line with the deviation in harvest volume. residue extraction is allowed on sandy soils, but not on clay-or fine silt soils. Further, stumps are not extracted on lichen-rich plots due to considerations for reindeer husbandry. Only spruce and pine stumps are extracted, stumps of all other tree species are left in the ground. The remaining forest, 17% of the area, was set aside for nature conservation without any woody biofuel extraction, in both settings. Nature reserves were left unmanaged. In voluntary set-asides and retention patches, some nature-conservation focused management was implemented, by applying a management strategy consisting of selection felling on 10% (voluntary set-asides) or 20% (retention patches) of the area, leaving the remaining area unmanaged.

Treatment Selection
After simulating alternative treatment schedules for each NFI plot, a set of optimization problems were formulated and solved with LP for both biofuel extraction settings (BAU and bioE). All problems consisted of assigning treatment schedules to each NFI plot so that the value of the objective function was maximized while the defined restrictions were met. All optimizations included a constraint specifying the demand for roundwood harvest, based on projections done in the latest national scenario analysis on potential sustainable harvest levels [26]. Annual roundwood harvest demands range from 6.1 in the beginning, to 6.9 million m3 (under bark) in the end of the planning horizon ( Figure 3). The roundwood demand includes firewood and unmerchantable stemwood, as a proportion of regular roundwood harvests. We did not include these volumes in our assessment of potential woody biofuels. We applied the same harvest demand constraint in all optimizations, allowing at most a +/−1% deviation from the harvest demand in each period. In addition, a constraint for evenness in woody biofuel extraction was included in all problems by accepting a maximal deviation of 20% in each period from average woody biofuel extraction over the entire planning horizon, which is in line with the deviation in harvest volume.

Maximum potential per indicator
The first set of LP problems aimed at determining the maximal potential production of the ecosystem services and biodiversity indicators (reindeer husbandry, carbon storage, recreation, old forest, mature broadleaf-rich forest, large diameter trees, and deadwood) for each scenario, subject to fulfilling a pre-defined level of roundwood harvest. The objective in each problem consisted of maximizing the total indicator value over all planning periods (separately for each of the indicators), subject to a demand for a certain amount of roundwood harvest for each five-year period.

Maximum Potential per Indicator
The first set of LP problems aimed at determining the maximal potential production of the ecosystem services and biodiversity indicators (reindeer husbandry, carbon storage, recreation, old forest, mature broadleaf-rich forest, large diameter trees, and deadwood) for each scenario, subject to fulfilling a pre-defined level of roundwood harvest. The objective in each problem consisted of maximizing the total indicator value over all planning periods (separately for each of the indicators), subject to a demand for a certain amount of roundwood harvest for each five-year period.

Trade-offs between Woody Biofuel Extraction and One Indicator at a Time
To identify trade-offs between the extraction of woody biofuel and the other indicators, a second set of LP problems was formulated and solved, identifying a set of Pareto-efficient solutions, i.e., solutions where the result for any one of the objectives cannot be improved without compromising Sustainability 2020, 12, 4089 9 of 20 the result for the other objective. The objective in each problem consisted of maximizing the total indicator value over all planning periods, subject to a gradually decreasing demand for woody biofuel, i.e., decreasing the woody biofuel demand step-wise (99%, 95%, 90%, 80%, etc.) until a further decrease did not change the results anymore.

Compromise Solution, including all Indicators
Finally, we used GP to identify a compromise solution, aiming at minimizing the maximum relative deviation from the maximal potential production (resulting from the first set of LP problems) for all of the indicators at the same time.
Each problem was formulated with linear objectives and constraints and solved with LP. Consequently, the models are examples of a standard model I formulation [45]. All optimization models were formulated within the Heureka system, using the ZIMPL optimization modeling language [46] and solved with Gurobi 7.0. See Appendix A for the mathematical formulation of the problems.

Results
The two biofuel extraction settings differ widely in the woody biofuel extraction potential: BAU produces at most 1 TWh woody biofuels per year, while the maximum potential is more than 4 TWh per year in bioE. Woody biofuel amounts also varied strongly depending on which indicator was maximized (Table 1). However, in the bioE management settings, woody biofuel amounts were always higher than the maximum potential in BAU, independent of which indicator was maximized.
There was a large variation in the potential delivery of the different ecosystem services and biodiversity indicator levels (Table 3). In both biofuel extraction settings, there are trade-offs between woody biofuel extraction and ecosystem services, as well as biodiversity (Figure 4). The trade-off is particularly strong between woody biofuel on the one hand and mature broadleaf-rich forest and old forest on the other, while carbon stock, deadwood volumes, and recreation are less affected. However, the shape of the trade-off curves-i.e., the almost linear part, where an increase in biofuel harvest results in relatively small losses in the other indicator-suggests that careful management planning allows for increasing woody biofuel extraction with only small losses in all biodiversity and ecosystem service indicators, up to a certain point (Figure 4). The bioE settings allowed for considerably higher woody biofuel extraction, while reaching the same maximum levels of ecosystem services and biodiversity indicators as in BAU, except for recreation. The compromise scenarios allowed for the extraction of 0.8 (BAU settings) and 3.3 TWh (bioE settings) of woody biofuel per year. At the same time, the compromise scenarios also led to higher average ecosystem service and biodiversity indicator levels compared to current conditions (given by the input data and reflecting the state of the forest in 2010), except for reindeer pasture and recreation index, which decreased in all optimizations compared to current conditions, for both biofuel extraction settings (BAU and bioE). This is pre-dominantly due to forests becoming continuously denser, which has a negative impact on ground lichen and recreation.  The C stock increased over time in all scenarios, and about half of the increase was due to an increase of C stock in tree biomass, and the other half in forest soils. Woody biofuel extractions affected the total C stock only marginally. The C stock increased over time in all scenarios, and about half of the increase was due to an increase of C stock in tree biomass, and the other half in forest soils. Woody biofuel extractions affected the total C stock only marginally.
The average annual woody biofuel extraction potential exceeded current extraction levels considerably in both the BAU and bioE management when maximizing biofuel, as well as in the compromise scenarios ( Figure 5). Most of the biofuel potential under the bioE settings comes from biofuel thinning, an assortment that is hardly harvested currently. The extraction potential in the bioE compromise scenario corresponded to approximately 27% (BAU compromise scenario: 6%) of the total energy consumption of about 12.2 TWh in the Västerbotten region in the year 2017 [47]. The extraction potential for logging residues was somewhat higher in BAU compared to bioE as the average rotation length for spruce-dominated stands is longer in BAU, which means that stands have more biomass when final-felled. The average annual woody biofuel extraction potential exceeded current extraction levels considerably in both the BAU and bioE management when maximizing biofuel, as well as in the compromise scenarios ( Figure 5). Most of the biofuel potential under the bioE settings comes from biofuel thinning, an assortment that is hardly harvested currently. The extraction potential in the bioE compromise scenario corresponded to approximately 27% (BAU compromise scenario: 6%) of the total energy consumption of about 12.2 TWh in the Västerbotten region in the year 2017 [47]. The extraction potential for logging residues was somewhat higher in BAU compared to bioE as the average rotation length for spruce-dominated stands is longer in BAU, which means that stands have more biomass when final-felled. The production strategy was the dominant strategy in both the BAU and bioE management settings, applied on more than half, and up to 82%, of the forest managed for wood supply in all scenarios ( Figure 6). The production strategy was most prevalent when the focus was on maximizing biofuel, which at the same time had the lowest prevalence of CCF. CCF was more commonly applied in all other scenarios, including the compromise solutions (both under BAU and bioE management settings), as it supports biodiversity and several ecosystem services better than even-aged management [48]. There were only small differences in the distribution of management strategy between the BAU and bioE management settings. For example, the no thinning strategy was less common in the compromise scenarios under the bioE settings, as strategies that allow woody biofuel thinning were chosen instead. Notably, in all scenarios, a variety of strategies was applied-pure production-oriented management is not optimal, not even for maximizing woody biofuel extraction. The production strategy was the dominant strategy in both the BAU and bioE management settings, applied on more than half, and up to 82%, of the forest managed for wood supply in all scenarios ( Figure 6). The production strategy was most prevalent when the focus was on maximizing biofuel, which at the same time had the lowest prevalence of CCF. CCF was more commonly applied in all other scenarios, including the compromise solutions (both under BAU and bioE management settings), as it supports biodiversity and several ecosystem services better than even-aged management [48]. There were only small differences in the distribution of management strategy between the BAU and bioE management settings. For example, the no thinning strategy was less common in the compromise scenarios under the bioE settings, as strategies that allow woody biofuel thinning were chosen instead. Notably, in all scenarios, a variety of strategies was applied-pure production-oriented management is not optimal, not even for maximizing woody biofuel extraction. Sustainability 2020, 12, x FOR PEER REVIEW 14 of 21 Figure 6. Proportion of management strategies in the BAU and bioE scenarios, when maximizing each indicator individually and in the compromise solution. Forest area set-aside fore nature conservation (nature reserves, voluntary set-asides, and retention patches) is not included.

Discussion
Society is gradually evolving from a dependence on fossil fuels to sustainable alternatives like wood fuels. For many countries, there are good opportunities to increase the production of wood fuel in combination with timber and pulpwood production. However, at increased production rates, the forest ecosystem and its ecosystem services and biodiversity will be affected. As a result, studies investigating the possibilities to increase the production of wood fuel in a sustainable way are important. In this study, we used long-term simulations of future forest development for different management and biofuel extraction strategies, combined with an optimization approach to analyze how much woody biomass forest ecosystems can provide in a sustainable way.
The results from this study demonstrate that there is considerable leeway for forest management in pursuing desirable objectives, while maintaining current roundwood harvest levels. With careful management planning, the extraction of woody biofuel (i.e., biomass from logging residues, stumps, and thinnings in young dense stands) could be increased considerably from current levels without a significant negative impact on most of the studied ecosystem services and biodiversity. However, reindeer grazing potential and recreation index decreased compared to current levels, in all scenarios. This is to a large extent due to the projected increasing density in future forests, which continues an on-going trend [49]. This could be alleviated with increased thinning intensity, which would reduce the potential for wood production and was not assessed in this study.
Even if these results seem promising from a production perspective, a few issues have to be discussed. First, the indicators included in our analysis did not capture all effects of woody biofuel extraction. Our assessment of reindeer grazing potential and recreation focused primarily on forest structure and did not consider the potentially positive effects of residue removal on accessibility for reindeer or humans. Another uncertainty concerns the potential negative impact of residue extraction on wood production in the long-term, due to the additional extraction of nutrients [17,50]. In the

Discussion
Society is gradually evolving from a dependence on fossil fuels to sustainable alternatives like wood fuels. For many countries, there are good opportunities to increase the production of wood fuel in combination with timber and pulpwood production. However, at increased production rates, the forest ecosystem and its ecosystem services and biodiversity will be affected. As a result, studies investigating the possibilities to increase the production of wood fuel in a sustainable way are important. In this study, we used long-term simulations of future forest development for different management and biofuel extraction strategies, combined with an optimization approach to analyze how much woody biomass forest ecosystems can provide in a sustainable way.
The results from this study demonstrate that there is considerable leeway for forest management in pursuing desirable objectives, while maintaining current roundwood harvest levels. With careful management planning, the extraction of woody biofuel (i.e., biomass from logging residues, stumps, and thinnings in young dense stands) could be increased considerably from current levels without a significant negative impact on most of the studied ecosystem services and biodiversity. However, reindeer grazing potential and recreation index decreased compared to current levels, in all scenarios. This is to a large extent due to the projected increasing density in future forests, which continues an on-going trend [49]. This could be alleviated with increased thinning intensity, which would reduce the potential for wood production and was not assessed in this study.
Even if these results seem promising from a production perspective, a few issues have to be discussed. First, the indicators included in our analysis did not capture all effects of woody biofuel extraction. Our assessment of reindeer grazing potential and recreation focused primarily on forest structure and did not consider the potentially positive effects of residue removal on accessibility for reindeer or humans. Another uncertainty concerns the potential negative impact of residue extraction on wood production in the long-term, due to the additional extraction of nutrients [17,50]. In the study, we assumed that residue extraction after thinnings reduces the basal area increment by 5% during the following 15 years, while residue extraction after final felling does not have any impact on future wood production. Wood ash recycling may counteract nutrient depletion and minimize negative effects of residue extraction on growth, at least on some soil types [51,52].
While both logging residue and stump extraction could be increased without large adverse effects on other ecosystem services and biodiversity, the largest potential and currently under-used resource is biomass from biofuel thinning. Thus, the early stage of rotations seems to give a high flexibility in management without jeopardizing future stand values. As a recent literature review indicated, biomass harvest in dense early thinnings enables early outtake of bioenergy, while simultaneously maintaining the stand structure's vertical heterogeneity and thereby supporting biodiversity [53]. The potential may also increase if cost-effective management systems are developed and implemented [54].
Our assessment was based on currently available technology for residue removal. The future development of technology and methodology for woody biofuel extraction operations may increase the available potential without compromising other forest values. For example, smarter extraction methods that are better able to separate residues into nutrient-rich parts-to be left in the forest-and nutrient-poorer parts-to be extracted-could lead to higher biomass extraction levels, better fuel quality and more nutrients left in the forest [55]. Another option is the development of precision forestry that leaves just enough of, e.g., deadwood/nutrients to keep the highest ecosystem service levels. In addition, our assessment did not account for transport distance from the forest to the industry. Energy wood and wood chips from harvest residues are predominantly transported by truck, which is not economically viable over distances longer than 100 km [15]. Long transport distances are currently one of the main factors behind the low residue extraction rate in northern Sweden. The future development of woody biofuel harvest is mainly influenced by economic preconditions, i.e., future energy prices and transport distances. Therefore, extraction levels may change quickly as a result of changes in energy prices or a shift from truck to train transport [12,15].
Our estimates of woody biomass potential from harvest residues are lower than those assessed in the latest Forest Impact Analysis-at most circa 0.5 TWh per 10 6 ha managed forest in our study compared to 1.8 TWh per 10 6 ha managed forest in northern Sweden [10,26]. However, the Forest Impact Analysis did not include any economic constraints for residue extraction, while our analysis included several such constraints, in particular by setting a lower limit for the proportion of spruce (50%) for forests where residues can be extracted.
The focus of our analysis was on the long-term impact of increased woody biofuel extraction, in particular the indirect impact through changes in management strategies. For example, a higher proportion of mixed forests can lead to the reduction in woody biofuel potential, since a minimum spruce proportion is needed for biofuel extraction to be economically feasible. Our results demonstrate that the increasing demands on the forest call for careful management planning. All scenarios had a variation in management strategies, i.e., the production strategy was not optimal on all sites even when the sole focus was on woody biofuel production.
The reliability of long-term projections naturally depends to a large degree on the quality of the used growth and yield models. The decision support system we used in our analysis has been shown to give reliable growth projections for time periods of 100 years [42], in particular for even-aged forests. Growth projections in forests managed with continuous cover forestry are less certain [56]. We did not project potential impacts of climate change, as our modeling approach only allowed for including climate change effects on growth, but not expected increased disturbance risks. Both the increased growth that is expected for northern European forests as a result of a warmer climate [57], and the need for salvage cuttings after disturbances may increase future biofuel potentials compared to our assessments.

Conclusions
To conclude, the development of novel and cost-effective management systems for biofuel outtake in early rotations may provide options for a significant supply of bioenergy to the emerging bioeconomy, while at the same time securing biodiversity and important ecosystem values in future stand developments. In addition, there is potential to increase the extraction of harvest residues and stumps while at the same time improving conditions for biodiversity and the amount of carbon stored in forest ecosystems compared to current levels. However, the suitability of forests for recreation and reindeer grazing may decrease from the current level due to increasing future stand density-irrespective of the intensity of wood fuel harvest.  x ij = [0, 1] ∀i ∈ I ∀j ∈ J i (A9) Equation (A1) specifies the objective function, i.e., to maximize the total value over all planning periods for indicator e. Equations (A2)-(A5) are connected to the timber and pulpwood demand in each period. Equations (A6) and (A7) are connected to the demand of evenness in woody biofuel extraction. Finally, Equations (A8) and (A9) ensure that all plots are assigned in total one treatment schedule.

Trade-offs between Woody Biofuel Extraction and One Indicator at a Time
The second set of problems aimed at identifying the trade-offs in each scenario between the extraction of woody biofuel and the other indicators subject to fulfilling a pre-defined level of roundwood harvest and to a demand of evenness in woody biofuel extraction.
The mathematical formulation of the second set of problems is as follows: Equation (A10) specifies the objective function, i.e., to maximize the indicator value over all periods. Equation (A11) specifies the gradually decreasing demand for woody biofuel. This means that the stated problem was solved a number of times for each indicator and scenario. Each time the value of α is gradually decreased from 0.99, 0.95, 0.90, 0.80, etc., until a further decrease in the demand of woody fuel did not change the results anymore. MaxZ w is the potential maximal extraction of woody biofuel. The value of MaxZ w is the result of solving the first set of Equations (A1)-(A9) with e ijp = B ijp .

Compromise Solution, including all Indicators
The third set of problems identified a compromise solution, aiming at minimizing the maximum relative deviation from the maximal potential production (resulting from the first set of LP problems) for all of the indicators, at the same time subject to fulfilling a pre-defined level of roundwood harvest and to a demand of evenness in woody biofuel extraction. This means that the problem described below was solved two times, once for both biofuel extraction setting.
The mathematical formulation of the last set of problems is as follow: Min Z c = max Z * * e − Z * e e = 1, ..E.
Subject to: and Equations (A2)-(A9). Equations (A12) and (A13) specify the objective function, i.e., to minimize the maximal relative deviation from the maximal potential production for all of the indicators at the same time. Z * * e is the normalized value of the maximal production of indicator e resulting from solving of the first set of problems and Z * e is the normalized value under decision variable x ij . Table A1. Explanation of variables, sets and constants included in functions 1-13.

Model Part Explanation
Variables x i j decision variable that takes a value between 0 and 1, indicating the proportion of NFI plot I that is assigned to treatment schedule j Sets P Set of all periods I Set of all NFI plots J Set of all treatment schedules J i Set of treatment schedules for NFI plot i, J i ⊆ J E Set of indicators S Set of scenarios Constants e ijp Amount of indicator e in in NFI plot i, treatment schedule j and period p T ijp Harvested timber volume in NFI plot i, treatment schedule j and period p M ijp Harvested pulpwood volume in NFI plot i, treatment schedule j and period p B ijp Wood fuel biomass extracted in in NFI plot i, treatment schedule j and period p T p Timber demand in period p M p Pulpwood demand in period p a Proportion of the maximal potential production of one the other indicators