Introduction

Forests are subject to a variety of conflicting interest at the individual (Blanco et al. 2015; Karppinen et al. 2020), institutional, and societal level (Eggers et al. 2019; Blattert et al. 2022). An often-dominating interest is the extraction of raw materials and direct economic benefits from timber production (Lindahl et al. 2017; European Commission 2018; Finnish Government 2022). Alternatively, the environmentally oriented forest owners and conservation oriented administrative strategies (European Commission 2020; Convention on Biological Diversity 2022) prioritize values other than timber. The clear conflict of interests requires the recognition and careful consideration of multiple objectives when planning the forest management. One possible way to search for synergies and reduce the unavoidable trade-offs (van der Plas et al. 2018; Eyvindson et al. 2019; Morán-Ordóñez et al. 2020) is the use of multi-objective forest management optimization (Li and Parrott 2016; Ezquerro et al. 2019; Díaz-Yáñez et al. 2021) that applies the idea of forests simultaneously providing multiple beneficial ecosystem services from timber extraction to carbon storage and maintenance of suitable habitats, i.e., forest multifunctionality (Jönsson and Snäll 2020).

In privately owned forest landscape, the harvesting intensity (the amount of harvested timber) will vary temporally and spatially depending on the economic demand and landowners’ willingness to sell timber (Hahtola et al. 1973). Resent timber removals in Finland (2016–2022) have been relatively high, at 91% of the average maximum maintainable annual felling potential (Ministry of Agriculture and Forestry of Finland 2022). For many administrative regions the regional harvesting intensity has been near or even above the estimated maximum felling potential (Hirvelä et al. 2023). Earlier studies have shown that a high harvesting intensity compromises non-timber ecosystem goods and services (Verkerk et al. 2014; Sing et al. 2018), and forest ecosystem resilience (Pohjanmies et al. 2021). Thus, to promote multifunctionality a notable reduction in the harvesting intensity below the current levels could be reasonable (Blattert et al. 2022).

The landscape level multifunctionality and biodiversity impacts of forestry are likely dependent on harvest intensity (Heinonen et al. 2017; Brockerhoff et al. 2017; Simard et al. 2020). Forestry operations reduce structures and habitats typical for unmanaged forests and thus negatively impacts many forest-dwelling species (Larsson and Danell 2001). In managed forests in Southern Finland the amount and diversity of dead wood, known to be important for biodiversity (Müller and Bütler 2010; Junninen and Komonen 2011), has been reduced by over 95% relative to natural forest stands (Mönkkönen et al. 2022). This has reduced the amount of suitable habitat for dead wood dependent species and accelerated their endangerment (Hyvärinen et al. 2019). In addition, the homogenic forest structure of commercial forests (Nilsson et al. 2001; Kuuluvainen and Gauthier 2018) is likely to further diminish the species diversity.

Forests additionally provide a wide range of other ecosystem services (ESS), like carbon storage, recreation value and collectable goods. Sustainable forest management provide a consistent flow of these services. These non-timber ESS are likely sensitive to variation in harvesting intensity in specific ways. For example, simulation studies have shown that increasing harvesting intensity decreases forest carbon storage (Sing et al. 2018; Seppälä et al. 2019; Simard et al. 2020) and the yields of some collectable goods while increasing the economic benefits of forestry (Heinonen et al. 2017; Eyvindson et al. 2018).

Restricted conservation resources require the use of cost-effective conservation actions. Previous studies indicate that alternative forest management regimes (determining the style, intensity and timing of harvest actions) differ in their impacts to ESS and biodiversity. For example, management regimes that promote continuous forest cover could offer more multifunctional forest landscape (Gustafsson et al. 2012; Pukkala 2016; Peura et al. 2018), whereas rotation forestry may offer more timber or economic benefits (Tahvonen and Rämö 2016; Peura et al. 2018). This variation in the effects of different management regimes offers an opportunity to cost-efficiently diminish the conflict between management objectives through multi-objective landscape level planning, as previously suggested by Mönkkönen et al. (2014) and Eggers et al. (2022).

In this study, we use the combination of forest growth simulation and linear programming optimization approach to evaluate the cost-efficiency of forest management to aid in ESS provisioning and biodiversity conservation by investigating how the selected management objective (maximal economic benefit vs. maximal multifunctionality) in combination with varying harvesting intensity affects the ability of a landscape level planning to combine conflicting objectives. Specifically, our aim was to answer following study questions:

  1. (1)

    What are the effects of harvesting intensity and management objective on forest multifunctionality?

  2. (2)

    How does the optimal combination of management regimes differ among economically- versus multifunctionality-targeted forestry planning and with varying harvesting intensity?

  3. (3)

    Does harvest intensity affect the cost-efficiency of multifunctional management?

Material and methods

Study area and description of the data

Our study area, Central-Finland, covers 16 700 km2 of land (Maanmittauslaitos 2023). Central-Finland is located at a transitional zone between south- and mid-boreal regions. In this region, wood-based industry has a high regional importance, and thus timber harvest intensity has been above the national average (95% of the maximum maintainable felling potential (Natural Resource Institute Finland 2023)). This together with the locally recognized need to comprehensively consider sustainability of forestry (Finnish Forest Center 2020), makes it suitable region for this kind of study. In Central-Finland approximately 82% of the land area is forest land (Natural Resource Institute Finland 2021), from which production forests cover more than 90% (Finnish Forest Center 2020). Protected areas, where no harvesting is allowed, cover approximately 3.7% of the land area (Natural Resource Institute Finland 2022). Most of the forests are privately owned.

We used open-source geographic information data on forest resources. These data are produced by the Finnish Forest Center and can be freely downloaded from (Metsäkeskus.fi). The downloaded data set covers 785 000 hectares (54%) of all the forestland in Central-Finland. The general quality of the data set can be considered relatively good. For example, in intermediate or mature forest stands the accuracy of estimates of basal area, diameter, height and total volume is within ± 20% of the correct value in 8 out of 10 stands (Finnish Forest Center 2016). However, the accuracy of remote sensing decreases for younger stands and areas with varying forest structure.

Clustering of the sampling frame

First, we divided the study area into two different geographical regions, south and mid-boreal. In both regions, we used stratified sampling where stratification was based on stand level characteristics (For details, see Appendix A). In each stratum, we targeted randomly selecting 10% of the stands but for practical reasons set a maximum of 100 stands per strata. This approach resulted 243 strata and 11 743 stands.

Description of simulation process and the extrapolation of simulated stands

All simulations in this study were made with an open-source forest management simulation program SIMO (Rasinmäki et al. 2009). Based on the simulated forest structure, we were able to determine the maximum maintainable harvest yield and the state of multifunctionality indicators.

For each forest stand, several alternative silvicultural management regimes were simulated. These alternatives included Business as usual (BAU) regimes and Continuous cover forestry (CCF) regimes. BAU regimes, also known as periodic rotation, are currently by far the dominant management regime in Central-Finland (Finnish Forest Center 2023). These regimes included artificial regeneration, commercial thinnings and the final felling, where all the trees (except retention trees) were removed from the felling site. In continuous cover forestry the final felling was not implemented but was replaced with more frequent felling of large trees and natural regeneration. Both BAU and CCF regimes included variation in the timing and intensity of management actions. This variation considered extended and shortened rotation length relative to BAU management, green tree retention and modified thinning intensities relative to BAU management. Also set aside (SA), abstention from all management actions, was always an alternative during the simulation process. In total, there were 22 alternative management regimes for each stand. For more detailed description of the management regimes, see the supplementary material in (Eyvindson et al. 2018).

For management regimes targeting periodic harvesting cycle, the predictions were based on Hynynen et al. (2002) growth model. As the Hynynen et al. (2002) model is developed based on periodic harvest, also an alternative model, suitable to predict forest growth under continuous cover forestry, was used (Pukkala et al. 2013). Both growth models are frequently used and compared (Heinonen et al. 2017; Triviño et al. 2023).

After the simulation process, simulated stands were extrapolated to cover whole data set. Extrapolation was based on a systematic approach, where from each sampled stand, the one having the lowest variation to a certain unsampled stand, based on mathematic model (Appendix A), was selected to represent specific stand. To achieve this, unsampled stands were compared to all sampled stands inside the stratum (maximum number of comparisons being 100). This systematic selection of sample stands improved the accuracy of the extrapolation.

Indicators of the forest multifunctionality

We used several different indicators to represent multiple forest-related values. A recent survey regarding landowner’s preferences revealed that economic orientation is the dominant value among landowners (Haltia and Rämö 2017). Because of this, the economic value was used as an independent preference when creating optimal solution. In this study, the forestry income was evaluated by estimating NPV with 4% discount rate. The calculation for NPV was based on the average tree species specific roadside prices from 2009 to 2018 for the simulated proportions of pulp wood and log wood (Mela: http://mela2.metla.fi/mela/tupa/index.php).

We used multiple indicators to reveal other than timber-based economic values of forests. For non-timber ecosystem service indicators, we used bilberry (Vaccinium myrtillus) yield and carbon storage. Bilberry yield (Miina et al. 2009) represented collectable goods, and hence hold an economic and recreational value. For carbon storage, an important atmospheric CO2 regulating service, the predictions considered both, the carbon stored in woody biomass (50% of dry biomass) of the growing stock and dead wood, as well as carbon in soil. Increased carbon storage is a sign of a forest absorbing and storing more carbon than releasing it, i.e., positive carbon balance. For estimating the soil carbon flux, separate models were used for mineral soils (Liski et al. 2005; Tuomi et al. 2009, 2011) and for peat lands (Ojanen et al. 2014).

For indicators of biodiversity, we used the amount and quality of dead wood as well as the amount of suitable habitat for six different vertebrate indicator species. Dead wood is known to be an important feature for many threatened species in boreal forests (Hyvärinen et al. 2019; Sandström et al. 2019), and due to forest management it has become a scarce resource in boreal production forests. In addition to dead wood amount, also the quality of dead wood is an important feature for many highly specialized species (Juutilainen et al. 2011; Sandström et al. 2019). Hence, the dead wood availability was measured as a function of total dead wood volume multiplied by the diversity of dead wood (Eyvindson et al. 2018). In this approach, the diversity was measured by the proportion of dead wood under different classes (species, diameter, decay class) as an inverse of Simpson´s diversity index (Triviño et al. 2017). Diversity weighted dead wood availability is high when a stand has a high amount of dead wood distributed evenly across all classes. Moreover, dead wood dependent species do not respond linearly to increasing dead wood availability, and particularly many threatened species only occur in forests that have dead wood more than 20 m3/ha (Junninen and Komonen 2011). Hence, we developed a specific function to describe a stand suitability for dead wood dependent species:

$$Q_{DW}~ = \left\{ {\begin{array}{*{20}c} {0,} & {if\;DW_{i}~ \le 5~m3/ha} \\ {0.067*DW_{i} - 0.33,} & {if~\;5\; < DW_{i}~ \le 20} \\ 1 & {if~\;DW_{i} > 20~m3/ha} \\ \end{array} } \right.$$
(1)

According to this function, stand suitability is zero on a stand with diversity-weighted dead wood volume (DWi) under 5 m3/ha, then increases linearly with the dead wood volume between 5 to 20 m3/ha and achieves the maximum suitability on a stand with dead wood volume over 20 m3/ha.

To achieving a more comprehensive estimate for the forest biodiversity, we included estimates for the habitat availability for six different vertebrate indicator species. This was achieved by determining species-specific habitat suitability index (HSI) for six species: the Capercaillie (Tetrao uralensis), Hazel grouse (Bonasa bonasia), Tree toed woodpecker (Picoides tridactylus), Lesser-spotted woodpecker (Dendrocopus minor), Long-tailed tit (Aegithalos caudatus) and Siberian flying squirrel (Pteromys Volans). The index relates to the probability of a stand being occupied by a selected species and varies from 0–1, where 0 represents unsuitable habitat with a smallest probability for a species to occupy a stand, and vice versa. Based on the HSI-index, areas with a relatively high probability of species occupancy (HSI > 0.5) were determined.

These selected indicator species cover a wide range of habitat requirements while also serving as umbrella species and hence give, together with dead wood, a reasonable estimate for the overall biodiversity. The ecological significance and habitat suitability modelling of the indicator species is described in more details in the appendix of Mönkkönen et al. (2014).

Description of the management scenarios, optimization and maximum maintainable yield

We analysed the effect of priority in forest management and the possibility to improve forest multifunctionality by creating optimal set of management regimes, utilizing two different objective functions:

1. Maximum NPV scenario represents economic objective of forest management when diverse set of management regimes can be used (Eq. 2).

2. Maximum MF describes non-timber multifunctionality objective for forest management when diverse set of forest management regimes can be used (Eq. 3).

$$\text{max NPV}_{0} = {\sum }_{i=1}^{m}\frac{{\sum }_{j=1 }^{h}\sum_{t=1}^{T}(\sum_{q=1}^{k}({P}_{q}-{C}_{q}){W}_{ijt}^{q}{X}_{ij}){\left(1+r\right)}^{5T-\left(10t-s\right)}+{PV}_{i}}{{(1+r)}^{5T}}$$
(2)
$$\text{max}MF= {\sum }_{e=1}^{E}\left(\frac{{M}_{e}-{M}_{e}*}{{M}_{e}^{*}-{M}_{e*}}\right)$$
(3)

Both objective functions are subject to a harvest constraint that requires harvest levels to meet a specific percentage of the maximum sustained harvest yield (Eq. 4).

$${\sum }_{i=1}^{m}{\sum }_{j=1}^{h}{W}_{ijt}^{q}{X}_{ij}-{F}^{q} *Z\pm \text{0,001}\ge 0, \forall t=1,\dots ,T, q\in k$$
(4)

where \({P}_{q}\) and \({C}_{q}\) are respectively the price and cost to harvest timber assortment q, \({W}_{ijt}^{q}\) is the quantity of the timber assortment q from set k, obtained from management unit i, according to management regime j at time period t, \({X}_{ij}\) is the decision variable for stand i to conduct management regime j, s is the time within the period when harvest actions are assumed to occur, T is the number of periods in the planning horizon, \({PV}_{i}\) is the productive value of the forest stand i at the end of the planning horizon, calculated using the models of Pukkala (2005) r is the discount rate applied, \({M}_{e}\) is the ecosystem service measure for criteria e of the set of E criteria, with \({M}_{e}^{*}\) and \({M}_{e*}\) are respectively the maximum and minimum possible ecosystem service measures possible for ecosystem service measure e, \({F}^{q}\) is the maximum maintainable yield of timber assortment q, and \(Z\) is a parameter to adjust the quantity of maximum maintainable timber yield.

To respond to the research questions, we evaluated the impact varying harvest intensities have on forest management according to both scenarios. We evaluated the impact of requiring specific harvesting intensities (Eq. 3, where Z is 60%, 80%, and 100%) of the maximal harvesting yield that can be sustained for long period of time. For research question 1, we assessed how harvest intensity impacts forest multifunctionality, both when prioritizing NPV (scenario 1) and MF (scenario 2). For research question 2, we evaluated how the change in management is impacted by the shift in priority and harvest intensity. While for research question 3 we evaluated the cost-efficiency of prioritizing MF over NPV maximization.

For both scenarios, the selected simulation time was from year 2016 to 2116, resulting in a 100-year time scale with 5-year steps. Relatively long-time scaling with shorter steps, promoted us to avoid conclusions that would seem optimal only in short time scale.

We used a linear programming optimization approach using the model I approach (Johnson and Scheurman 1977) to select the proportion of the stand to be managed according to the 22 simulated management alternatives. The optimization process used two different objective functions: 1) to maximize net present value or 2) to maximize the multifunctionality by producing the best possible compromise between different non-timber multifunctionality indicators (dead wood, carbon, bilberry and six biodiversity indicator species). Both objective functions were constrained to provide a specific proportion (60%, 80%, 100%) of the maximum maintainable harvesting yield.

The specific optimization models can be found in Appendix B.1. The optimization process creates a Pareto optimal solution where no objective could be improved without diminishing the outcome of at least one other objective.

The Maximum maintainable harvest was determined by optimization problem that sustained the maximum even flow of timber (See Appendix B.1).

Our analysis considers local forest landscape comprehensively including commercial forests as well as protected areas. Protected area cover in forest land in Central-Finland is 3.7% (Natural Resource Institute Finland 2022).

The maximum multifunctionality (max MF) is the optimal combined result for seven different biodiversity measures and two ecosystem service criteria (described below). The max MF is determined through an objective function that maximises the normalized sum of these nine different multifunctionality measures under the same constraints (Equations B.2-B.7, Appendix B.1) with max NPV (For details see Appendix B.2). To ensure equal consideration of all ESS in the analysis, no ESS specific weights are used.

Comparisons

Following simulation and optimizations multiple comparison was conducted to reveal the effect of harvest intensity and set management objective. First, to reveal the effect of set management objective, indicator specific relative gains were calculated considering the change from NPV to MF management with equal harvest intensity. Second, to reveal the effect of harvest intensity, indicator specific gains from different harvesting intensity scenarios with consistent management objectives, were calculated. Third, to reveal the combined effect of change in management objective and altered harvest intensity, indicator specific gains from changing NPV management objective with higher harvest intensity to MF management objective with reduced harvest intensity were calculated. Fourth, to enable cost-efficiency comparisons (total relative MF gain/opportunity cost), the opportunity cost for reduced harvesting intensity and changed management objective was calculated based on the difference in NPV. In all comparisons the total multifunctionality effect is the sum of biodiversity (BD) effect and ecosystem service (ESS) effect.

We calculated cost-efficiency based on two equations. First, the indicator specific relative gain as percentage (GMF) was calculated according to equation:

$${G}_{MF}=\left(\frac{{MF}_{end }- {MF}_{Start}}{{MF}_{Start}}\right)\cdot 100$$
(5)

where MFend is the indicator specific value at the end and MFstart is the equal value at the start. Second, as an estimate of cost-efficiency (CE), a measure of relative MF gain per euro was calculated (Eq. 6).

$$CE=\frac{{G}_{MF}}{{NPV}_{end}-{NPV}_{Start}}$$
(6)

Results

The effect of harvesting intensity and selected forest management objective

The maximum net present value of forest harvesting was 6000 €/ha over the 100-year planning horizon irrespective of the management objective (Fig. 1). The NPV was linked with the shift in harvesting intensity. As intensity was reduced, NPV reductions occurred with both scenarios, however with the NPV objective, the reductions were smaller, i.e., 40% or 20% reduction in harvesting intensity reduced NPV by 32% or 14%, respectively. The cost for multifunctional management is relatively low (90 €/ha or − 1.5%) at the 100% harvest intensity level. At 80% and 60% harvest level, NPV was 447 €/ha or 8.7% and 560 €/ha or 13.7% lower with multifunctional objective than with NPV-targeted objective.

Fig. 1
figure 1

The net present value for varying harvesting intensities and management objectives. The result is an average value over the whole simulation period (100 years). Stars represent 100%, triangles 80% and dots 60% harvesting intensity. MF (multifunctionality) and NPV (net present value, dashed line) indicate management objectives

The variation in harvesting intensity affected multiple ecosystem service indicators (Figs. 1, 2, 3 and 4). The exact magnitude and direction of the impact varies among indicators and depends on the forest management objective (Table 1). For most indicator species, decreasing harvesting intensity will result in an increase in habitat availability (Figs. 1, 2, 3 and 4). This is especially true when considering the multifunctional management scenario, that showed largest positive change with decreasing harvesting intensity from 100 to 60% in the habitat availability for the Siberian flying squirrel (habitat increment 328%) and Capercaillie (habitat increment 185%). When forestry management targeted primarily economic benefits (NPV) reducing harvesting intensity typically resulted in much smaller increase in habitat availability, and in some case no biodiversity benefits at all (Fig. 2). This pattern is evident also in dead wood density: with multifunctional management decreasing harvesting intensity from 100 to 60% resulted in 167% increment in the availability of dead wood rich forests, while with NPV targeted planning there was an opposite effect (reduction of dead wood rich forests) (Fig. 3 and Table 1).

Fig. 2
figure 2

Average suitable habitat area where HSI > 0.5 for six different biodiversity indicator species as a proportion of forest land. Figures show the whole simulation period (100 years). In all panels X-axis represents different timber harvesting intensities. Stars represent 100%, triangles 80% and dots 60% harvesting intensity from maximum maintainable yield. MF (multifunctionality) and NPV (net present value, dashed line) indicate management objectives

Fig. 3
figure 3

Amount of high-quality dead wood area. For the dead wood, the results are based on the average volume of diversity weighted dead wood in a hectare from which the dead wood quality index Q is determined. Q = 1 relates to the total volume of dead wood being equal to the area (ha) with dead wood volume ≥ 20 m3. Q > 0.5 represents areas with notably larger dead wood volume than on average in Central-Finland. Picture show the result for the whole simulation period (100 years). X-axis represents different timber harvesting intensities. Stars represent 100%, triangles 80% and dots 60% harvesting intensity from the maximum maintainable yield. MF (multifunctionality) and NPV (net present value, dashed line) indicate management objectives

Fig. 4
figure 4

The effect of harvest intensity and forest management objective for two different ecosystem service indicators. For bilberry (a panel), the result is an average berry yield as kg/ha (Y-axis). Panel b shows the effect of the harvest intensity to the average value of carbon stored in soil and in wood biomass as tons per ha (Y-axis). In both panels X-axis represents different harvesting intensities. Stars represent 100%, triangles 80% and dots 60% harvesting intensity from maximum maintainable yield. MF (multifunctionality) and NPV (net present value, dashed line) indicate management objectives

Table 1 Cost efficiency of biodiversity enhancement. Table illustrates the cost efficiency of multifunctional (MF) management (A), reduced harvest intensity when net present value (NPV) is targeted (B), reduced harvesting intensity while using multifunctional management (C) and the combination of reduced harvesting intensity and multifunctional planning (D) for improving biodiversity and the production of non-timber ESS in the landscape with varying harvest rate. Columns E, F, G, H, I, J and K follow the same logic with intermediate and lowest harvest intensity

Non-timber ESS respond to variation in harvesting intensity in idiosyncratic ways (Fig. 4). Carbon storage decreased linearly in both forest management scenarios with increasing harvesting intensity. Bilberry yields showed no response to variation in harvesting intensity with multifunctional management but increased non-linearly with increasing harvesting intensity with NPV targeted management.

The optimal set of forest management regimes

For all scenarios, the optimal outcome is a combination of multiple management regimes, and no single management regime can be considered optimal for all the objectives (Fig. 5). All the optimal solutions are composed of both BAU and CCF regimes. The relative share of CCF regimes is higher, and that of alternative rotation forestry regimes lower, when the forestry planning is targeting multifunctionality than with NPV targeting management. The relative share of CCF regimes decreases with decreasing harvesting intensity, particularly with NPV objective.

Fig. 5
figure 5

Relative share of different forest management regimes in varying harvest intensity and management objectives. Y-axis represents proportions (%) of management regimes for different harvesting intensities and optimization objectives (X-axis). NPV stands for Net Present Value and MF for Multifunctionality. SA refers to Set Aside (no management), BAU = Business as Usual (clear cut based periodic rotation), Ext.Rot = Extended rotation length, GTR = Green Tree Retention, Stand.Rot = Standard rotation length with thinning modifications, Short.Rot = Shorter than BAU rotation length, CCF = Continuous Cover Forestry (no clear cutting). Management regimes are grouped with details included in method section

At the level of maximum maintainable harvesting intensity, the relative share of set-aside (SA) areas is very low (< 0.5% in both NPV and MF models) and increases with decreasing harvesting intensity. This increase is stronger with multifunctional planning.

The cost-efficiency of multifunctional forest planning and reduced harvest intensity for improving forests biodiversity and production of non-timber ecosystem services.

The analysis showed harvest intensity specific cost-efficiency differences associated with varying management objectives and reduced harvest intensity (Table 1). Decreasing harvest intensity (columns B, C, F, G, J and K in Table 1) is always beneficial to multifunctionality, and multifunctional management always provides higher multifunctionality than NPV targeted management (A, E, I). In general, larger multifunctionality gains can be achieved by shifting planning focus from NPV to multifunctionality than by simply reducing harvesting intensity (A vs. B and E vs. F). For most indicators (excluding dead wood), the largest multifunctionality gains were achieved by combination of MF-planning and reduced harvest intensity (D, H and L).

The most cost-efficient (evaluated as the relative gain for MF divided by the change in NPV) enhancement of multifunctionality was achieved through the change of management objective (A, E and I). The relative cost-efficiency of shifting the planning objective is harvesting rate specific, being lowest in intermediate harvest intensity and increasing with lowering or increasing harvest rate. The combination of reduced harvest rate and changed management objective (D, H and L) reached intermediate cost-efficiency while decreasing harvest rate with fixed management objective (NPV or MF) produced lowest cost efficiency (B, C, F, G, J and K).

Discussion

The effect of varying harvesting intensity and the relative multifunctionality loss of NPV management objective

Our results indicate that both the harvesting intensity and management objective affect forest biodiversity, the production of non-timber ecosystems services and the economic benefits of forestry, as previously shown by (Verkerk et al. 2014; Biber et al. 2015; Eyvindson et al. 2018). For timber related economic benefits, with multifunctional management, the relative change in NPV was proportional to the level of reduction in harvested timber. However, with NPV targeted management the reduction was less than proportional (-40% harv. rate, -32% NPV), indicating how prioritizing a single optimization target can soften the impact of changing constraints.

Decrease in harvesting intensity will result in an increase in biodiversity (Fig. 2, 3 and Table 1). The effect is not uniform for all indicator species, but the increment is notable, and the general trend holds with both management objectives. This result is consistent with earlier studies (Esseen et al. 1997; Biber et al. 2015; Heinonen et al. 2017; Naumov et al. 2018; Eggers et al. 2019) that have identified a positive biodiversity effect related to the decreased forest harvesting intensity, and thus concluded that increasing harvest intensity will lead to a loss of biodiversity.

Commercial forest management has been recognized as the main reason for the decline of many forest species in Finland (Hyvärinen et al. 2019). Hence, our results confirm the trade-off between high intensity of commercial forest use and forest biodiversity. However, not all species respond similarly. For instance, Long-tailed tit and Lesser-spotted woodpecker showed increased habitat availability with high harvest intensity when management was targeted towards NPV (Fig. 2). This is likely caused by the increased share of CCF management in NPV 100%—scenario (Fig. 5). Since both species are known to benefit from CCF management (Peura et al. 2018) the increased harvest intensity increases the amount of most suitable management for both species. Nevertheless, both species have the highest habitat availability with multifunctional management and reduced harvesting intensity.

Multifunctional management offers a highly efficient way to enhance sustainability of forestry (Table 1). The positive multifunctionality effect of shifting the management focus can even exceed the negative multifunctionality effect of increased harvest intensity (Appendix C), as also shown by (Eggers et al. 2022). From the practical forest management perspective, this highlights the importance of planned, science-based actions (Sutherland et al. 2004) for biodiversity and the production of non-timber ESS. This potentially beneficial role of multifunctionality targeted management is also shown in earlier studies (Pukkala 2016; Tahvonen et al. 2019; Moor et al. 2022). Consistent with Eyvindson et al. (2018), our results showed the highest potential for multifunctional management when the harvest intensity was not maximized.

Despite the evident benefits, multifunctional management of commercial forests needs to be supported by effective conservation network, as suggested by (Nolet et al. 2018). This is because multifunctional management results uneven habitat gains benefitting only a subset of species and habitats (Siberian flying squirrel, Lesser spotted woodpecker, Capercaillie and Dead wood), while other species remain not supported. In addition, these large proportional gains, most notable in dead wood, are easiest to achieve when the reference point is low, as is the case with the habitat availability of our indicators (Fig. 2, Fig. 3 and Supplementary dataset). Despite high relative multifunctionality gains, the maximum coverage of high dead wood quality stands achievable through multifunctional management at the maximum harvesting intensity, is about 1% of the total forest area. Although, due to the lack of spatial variation in the amount of dead wood at the start point of the simulation period, our result might underestimate the amount of suitable habitat for dead wood dependent species. Nevertheless, since dead wood is scarce resource in commercial forests (Sandström et al. 2019; Mönkkönen et al. 2022) and the conservation network coverage is restricted, the general conclusion of predicament condition for dead wood dependent species holds without further actions. The same is true for other species as well, since with the MF 100%—scenario, 5 out of 6 indicator species reached the amount of suitable habitat less than desirable (< 20–30%, Fig. 2) for the long-term persistence (Andrén 1994; Hanski 2015).

Our analysis showed a relatively linear negative relationship between carbon storage and harvesting intensity (Fig. 4) and revealed that there is practically no potential (0.2%) to increase carbon stock through multifunctional management, with the current regional harvesting intensity of 95% (Natural Resource Institute Finland 2023). In contrast, reducing harvest rate resulted larger increases in carbon stock. This is in line with e.g. Biber et al. (2015), Triviño et al. (2017), and Sing et al. (2018) who have shown that the least intense harvesting scenario results in the largest carbon storage. Similarly, Seppälä et al. (2019) demonstrated a trade-off between timber harvesting and forests carbon regulating service, when considering current displacement factors of wood-based products in Finland and reasonable time scale for climate change mitigation. Although our result is consistent with previous studies, it should be noted that larger increment in carbon storage could possibly be available through independent optimization of carbon storage (Jörgensen et al. 2021; Walker et al. 2022; Mäkelä et al. 2023), as biodiversity and carbon storage can have partially conflicting relationships (Triviño et al. 2017; Sabatini et al. 2019; Mönkkönen et al. 2022). Nevertheless, in summary, the abstention from harvesting maximum maintainable annual felling potential seems to be the more determinant factor increasing carbon stocks in boreal forest in comparison to changes in management strategies, as also shown by (Mäkelä et al. 2023).

Optimal set of management regimes

Our results showed, that at the regional scale, the optimal management plan always consists of a diverse set of different management options (Fig. 5). This is well in line with other similar studies examining the optimal way of producing forest related ecosystem services. For example, Nolet et al. (2018), Díaz-Yáñez et al. (2019), Eyvindson et al. (2021), and Eggers et al. (2022), have all reported about the benefits of using a diverse set of forest management. Even at the 100% harvest intensity and with NPV-targeted planning, continuous cover forestry regimes constitute almost 50% of the total harvesting in the optimal solution (Fig. 5). The optimal combination differs drastically from the current situation where the CCF covers about 2% of all the fellings in Central-Finland (Finnish Forest Center 2023). This difference may be due to CCF only recently being legally allowed in Finnish forests and could be due to biases in how CCF operations are modelled. These issues are possible to assess and quantify as more people enact CCF in the forest and as more data becomes available under CCF forest management. The relative utility of CCF stems from more stable timber flow, larger proportion of log wood and lower regeneration costs than in rotation forestry (Peura et al. 2018). Also, Tahvonen and Rämö (2016) have showed that CCF may be economically more profitable than rotation forestry particularly at high interest rates such as in our case (4%).

Our results showed an important role of set-asides (SA) when targeting multifunctionality, but this is only possible if we give up on striving for maximum maintainable annual felling potential. For example, it is optimal to set aside one third of the forests at 60% level of harvest intensity when maximizing multifunctionality. A high proportion of SA areas when targeting multifunctionality is conceivable in the light of previous studies, which have shown the high relative importance of SA management option for providing ecosystem services and suitable habitats for many species (Nolet et al. 2018; Tahvonen et al. 2019; Lamothe et al. 2019; Pohjanmies et al. 2021).

The cost and the relative cost efficiency of MF-management and reduced harvest intensity

Our analysis illustrated harvest rate specific cost-efficiency differences associated with changing management objective and / or reducing harvest intensity (Table 1). Since the economic cost of enhanced multifunctionality is tightly linked to extracted timber (Fig. 5), the most cost-efficient improvement is always achieved through the change of management objective (NPV to MF). The potential of multifunctionally targeted management to cost-efficiently promote multifunctionality has also been recognized earlier (Díaz-Yáñez et al. 2021; Eggers et al. 2022; Mazziotta et al. 2023). The cost of multifunctional management was lowest, with maximal harvesting intensity, but notably increased with decreasing harvest intensity (Table 1). The modest monetary cost with maximal harvest intensity might make MF-management suitable option even for relatively economic oriented landowners.

Similar to increased costs, the gains in multifunctionality from shifting the management focus (NPV to MF) were much higher at intermediate and low harvest rate than at 100% harvest rate (Table 1, cf columns A, E and I). The relative cost-efficiency of multifunctional management was however dependent on harvest intensity. The multifunctionality gain per invested € first decreased with decreasing harvest rate from 100 to 80% but peaked at 60% harvest rate. This peaking, however, is caused by dead wood, that reaches high relative gains with multifunctional management at 60% harvest rate (Table 1, column I), although the absolute area of high-quality dead wood areas is lower in multifunctionality targeted scenario with 60% harvest rate in comparison to 80% harvesting rate (Fig. 3 and Supplementary dataset). Nonetheless, careful multifunctionally targeted multi-objective management is a particularly attractive societal option if there is a need to maintain a high harvest intensity. However, the largest gains in multifunctionality were generally achieved through the combination of multifunctional management and the reduction in harvest rate. Our result suggests that if the society accepts reduction in harvest rate, which is necessary to safeguard biodiversity and non-timber ecosystem services, it is recommendable to do in combination with the changed management objective.

Multifunctional management and/or reduced harvest intensity cause economic costs. Since the gains and costs of multifunctional management are higher with intermediate and low harvest rate, a monetary compensation from society might be needed to create an incentive for the large-scale introduction of multifunctional management. For the society, partial or full monetary compensation could cost-efficiently promote forest biodiversity and the flow of non-timber ecosystem services. However, even with monetary compensations, the large-scale introduction of landscape level multifunctional management can be difficult to obtain, as forest owners make decisions based on their personal preferences, instead of landscape level optimality. Therefore, multifunctional management should most likely be targeted to stakeholders with relatively large forest properties and willingness to provide multiple goods and benefits (Holm 2015), such as some common forests and municipalities. Multifunctional management could be especially suitable for buffer zones, ecological corridors of existing conservation areas or other areas where harvesting should be restricted for environmental reasons. Particularly, other effective area-based conservation measures (OECM, see: Convention on Biological Diversity (CBD) 2018, (COP 14)) that are areas out site current protected network, but are managed in a way offering positive, long-term outcomes for biodiversity and ecosystem services, might be suitable areas for multifunctional forestry.

Conclusion

Our results indicate a strong trade-off between timber related and non-timber ecosystems services (carbon storage) and biodiversity. Maximizing the even flow of timber to the maximum maintainable level will increase the per hectare net present value, as expected (Biber et al. 2015). Simultaneously, maximizing harvest rates will result in decreases in the amount of non-timber related ecosystem services and in forest biodiversity. For some biodiversity indicators, the landscape level multi-objective planning offers a way to diminish or even overcompensate this decline, with relatively low or even positive opportunity cost (Appendix C). However, without an expansion of conservation area network, mere shift in planning focus on MF may not sufficiently increase the ecological quality of forests to ensure the long-time persistence of all the species. Thus, reducing the harvesting intensity is required to also allow for more set asides.