Effect of forest management choices on carbon sequestration and biodiversity at national scale

Forest management methods and harvest intensities influence wood production, carbon sequestration and biodiversity. We devised different management scenarios by means of stakeholder analysis and incorporated them in the forest growth simulator PREBAS. To analyse impacts of harvest intensity, we used constraints on total harvest: business as usual, low harvest, intensive harvest and no harvest. We carried out simulations on a wall-to-wall grid in Finland until 2050. Our objectives were to (1) test how the management scenarios differed in their projections, (2) analyse the potential wood production, carbon sequestration and biodiversity under the different harvest levels, and (3) compare different options of allocating the scenarios and protected areas. Harvest level was key to carbon stocks and fluxes regardless of management actions and moderate changes in proportion of strictly protected forest. In contrast, biodiversity was more dependent on other management variables than harvesting levels, and relatively independent of carbon stocks and fluxes. Supplementary Information The online version contains supplementary material available at 10.1007/s13280-023-01899-0.


S1. Mortality model
The stand-level mortality models developed by Siilipehto et al. (2020) where implemented in PREBAS.The models were developed using permanent national forest inventory data from Norway, Sweden and Finland.The models are logistic functions that predict mortality in two steps: the first model predicts the probability of survival ( ), the second model predicts the proportion of basal area that remains after a mortality event ( ).
The structure of the model is the following: where P is the probability of no mortality (in  ) or the proportion of survived trees (in  );   is the linear combination of predictor variables times the parameters of estimated in Siipilehto et al. (2020).
The parameters and predictor variables of  are reported in Table 5 of Siipilehto et al. (2020).Similarly, the details of the linear part   for the proportion of basal area in surviving trees is reported in Table 6 of Siipilehto et al. (2020).
Mortality is then estimated using in combination the two models: where  is the basal area of the dead trees and  is the stand basal area of the growing stock.
The mortality model described above was applied only to the protected areas, jointly with the Reineke mortality model (Reineke, 1933), at annual time step.In the managed forests, only the Reineke function was used (Minunno et al., 2019).

Overview
The purpose of the ground vegetation (GV) module is to estimate the contribution of ground vegetation to the modelled ecosystem carbon fluxes, including (1) photosynthesis, (2) respiration, and (3) litter fall (Figure S1).We first estimate the percentage cover of different functional groups of ground vegetation, then convert this to above-ground and below-ground biomass using conversion equations.Turnover can be related to biomass, and further, with the simplifying assumption that ground vegetation is approximately at steady state, growth and respiration can also be approximated.
There is established empirical evidence that the percentage covers of different GV types depend on site fertility (Tonteri et al. 2005) and light reaching the bottom layer (Strengbom et al. 2004), which vary with the development stage of the stand (Tonteri et al. 2005).The PREBAS GV module therefore uses site fertility class and the proportion of light absorbed by the tree canopy as explanatory variables of the percentage cover of ground vegetation types.

Percentage coverage
The core of the GV module is the model of percentage covers of different functional groups of GV.In order to keep the model simple, we only use three groups, some of which are combinations of groups in a more detailed classification based on the light and nutrient demand of different plant types: (1) grasses and herbs ( ), ( 2) shrubs ( ), and (3) mosses and lichens ( ).The models are based on a qualitative understanding of the light dependence of the groups, with level parameters determined separately for different site fertility classes (Tonteri et al. 2005).Thus, we assume that grass and herb coverage increases linearly with the light entering the understorey, whereas the percentage coverage of mosses and lichens decreases with increasing light.Shrubs are assumed to reach a maximum percentage coverage at intermediate light levels and minima at full light and no light.These dependencies are described by the following equations: where the dependence on site fertility is thought to be embedded in the parameters ( ,  ,  ,  ,  ).2011) used a subsample of data in Tonteri et al. (2005) for a more detailed analysis.Her data set covers three fertility classes, including 180 sub-xeric, 423 herb-rich and 251 grove-like sites on upland soils in southern Finland.Tähkälä (2011) analysed the data with respect to stand age, stand basal area and time from the latest harvest (either clear cut or thinning).A general feature of the data was that the variability was very high, the standard deviations being approximately equal to the means in all analyses.We used the results from Tähkälä (2011) as our primary source but extended the results from three to five fertility classes using Tonteri et al. (2005).

Biomass, turnover and leaf area
Based on empirical findings (Muukkonen et al. 2006, Lehtonen et al. 2016), we assumed an exponential relationship between percentage cover of species group  and mass per unit area.In order to determine leaf area and turnover rates, the mass was divided into above-ground ( ) and below-ground ( ) components: where  ,  ,  and  are parameters for the above-ground and below-ground components, respectively (Table S2.2).The parameter values were obtained by simplifying the mixed model equations estimated by Lehtonen et al. (2016) separately for grasses, herbs, shrubs, mosses and lichens.Here we combined the original functional groups of grasses + herbs and mosses + lichens.This was done by first calculating the mass estimates for the original classes, then fitting the above simplified equations to the combined classes.
Compartment specific turnover rates were obtained from Lehtonen et al. (2016), again by combining original classes when needed.This assumed constant annual turnover rates for all above-ground and below-ground biomass components, such that turnover  ( ∈ (, , );  ∈ (, ) ) is with  denoting the component-specific turnover rate (yr -1 ).
Table S2.2.Parameters for different plant types. is expressed in kg DW ha -1 and turnover rate is yr -1 .
Leaf area  , was assumed to be obtained from the above-ground biomass through multiplication with specific lear area,  (m 2 kg -1 ): The specific leaf area of shrubs and grasses + herbs was based on Palmroth et al. 2014.They estimated total biomass, root to shoot ratio and foliage to biomass ratio of bilberry and cowberry in Rosinedal pine forest, northern Sweden.We used this information to estimate foliage area per above-ground biomass, which is shown in Table S2.2 as   .We estimated the value for grasses and herbs as somewhat higher than that for shrubs.The leaf area for mosses and lichens is a guess that produces approximately the share of photosynthesis of mosses and lichens given by Laine et al. (2016).

Contribution to ecosystem GPP and respiration
For PREBAS, we add the LAI of GV to forest canopy LAI and use this as forest ecosystem LAI.This is used in estimating ecosystem and canopy GPP in PREBAS (Peltoniemi et al. 2015).With the above parameter values, the model gives an almost linear relationship with tree canopy  and the combined proportion of absorbed radiation,  , (Figure S2.4).
The contribution of GV to canopy respiration is estimated using the assumption that GV is approximately at steady state, such that annual growth equals annual litter fall.This allows us to estimate GV respiration as where GPP and litter fall are evaluated as explained above.
Because the GV model is not dynamic, but GV biomass is determined by the current value of  , the above equation may lead to negative respiration at times of thinning or clear cuts.This is avoided by restricting respiration to non-negative values.

Model evaluation
Given that the core of the model is based on phenomenologically justified equations (Eqn S2.1) with nonsystematic parameterization, we wanted to evaluate its performance against independent studies.Two suitable studies were found: Kulmala et al. (2011), Tupek et al. (2016).Kulmala et al. (2011) presented data from an age series of five pine stands representing herb-rich or subxeric sites in southern Finland.They also reported the canopy all-sided leaf area, which can be used here to test our results.We use  = 0.18 for the canopy (Duursma and Mäkelä 2007) to obtain the respective  values.We calculate the contribution to total  of the ground vegetation as follows: where  is the proportion of light absorbed by the tree canopy and  is total photosynthesis by ground vegetation, both as estimated by Kulmala et al. (2011), and we use  = 1200 g m -2 yr -1 (Minunno et al. 2016).
A comparison of the above-ground biomass of ground vegetation reported by Kulmala et al. (2011) and the proposed model with sub-xeric sites shows that in this case, the model underestimates shrub biomass but slightly overestimates the other biomasses (Figure S2.4a).This could be indicative of the wide variability of ground vegetation between individual stands.Nevertheless, the predicted contribution of GV to  , was very much in line with that calculated from Kulmala's observations (Eqn S2.5) (Figure S2.4b).This relationship appears independent of site fertility class in the model.2016) presented a statistical model of above-and belowground biomass of ground vegetation types, divided into grasses, herbs, shrubs, mosses and lichens.The explanatory variables were stand variables available from NFI measurements, including latitude.The key variables in all models were the main tree species (pine or spruce) and latitude.However, there was a lot of uncertainty in the prediction, the goodness of fit of the models ranging from 0.1 to 0.3.
We used open-access data from the Swedish national forest inventory (Skogsdata 2020) as input to the process-based growth model PREBAS (Minunno et al. 2019) to assess site fertility and estimate  .The input data required by our model was total basal area and its proportional share among species, and mean height and mean diameter per species.We used these to estimate the biomass of the three groundvegetation species groups included in the present model.We then compared these corresponding predictions with the model by Tupek et al. (2016), which could also be evaluated for the sites.
Because of a larger number of explanatory factors, the variability in values predicted by the Tupek model was much larger than that in PREBAS combined with the current GV model.However, the order of magnitude and the overall range of variability were similar in both models (Figure S2.5).

Concluding remarks
The GV module of PREBAS is a first attempt to include GV-related carbon fluxes in the PREBAS model for regional-scale applications.Based on the large-scale summary type data used (Tonteri et al. 2005, Tähkälä 2011, Lehtonen et al. 2016), the model is largely qualitative and can be expected to provide results of appropriate order of magnitude only, rather than accurate site-specific estimates.However, the comparisons with previous empirical studies (Kulmala et al. 2011, Tupek et al. 2016) suggested that the overall order of magnitude and relations between different functional groups were consistent.Concurrently, both test data sets showed somewhat higher shrub biomass than could be predicted from our model building data, but the modelled patterns with respect to light availability agreed with data.Our model showed more variability in grasses and herbs than Tupek's model both with site type and light availability, probably reflecting the fact that these variables were not included in Tupek's explanatory factors.The understorey contribution to ecosystem GPP through light capture was remarkably similar to the results by Kulmala et al. (2011) and rather insensitive to site type, suggesting that after clear cut ground vegetation can capture 30% of the total light capture of closed stands, the ratio declining steadily as the tree canopy becomes denser.However, the results need further testing against comprehensive empirical data, such as measurements available from eddy covariance sites.

S3. Habitat Suitability Index (HSI) of biodiversity indicators
We used six previously published habitat suitability indices (HSI) for five bird species (capercaillie (Tetrao urogallus), hazel grouse (Bonasa bonasia), three-toed woodpecker (Picoides tridactylus), lesser-spotted woodpecker (Dryobates minor), long-tailed tit (Aegithalos caudatus)), and one mammal species (Siberian flying squirrel (Pteromys volans)) (Mönkkönen et al. 2014).The HSI attained values between 0 and 1 and were calculated by inserting the stand structural variables produced by the PREBAS simulations, such as tree species absolute and relative volume, basal area, density, age, and mortality or presence of dead trees (Table S3.1), to the species-specific HSI formulas described below.The value is close to one when the total volume living trees is >200 m 3 /ha and basal area (m 2 /ha) of fresh dead-wood is >2.5m 2 /ha.The latter threshold level translates into about 20m 3 /ha of deadwood in a mature stand (vtotal > 200). Lesser

S5. Stakeholder meeting
A stakeholder meeting was held by the IBC-CARBON project on 7 November 2018, where the invited 13 experts represented universities, research institutes, ministries, forest owners' associations, forest consultancy companies, companies managing nationally owned forests, and different NGOs.In addition, 9 project members and one external facilitator took part in the meeting.
The participants were sent some pre-reading material with background about the research methods and objectives of the meeting, which was expanded in introductory presentations.This included a general presentation of the research project and its objectives, and an introduction to the PREBAS model and its intended application in the project.
The main setup of the study was designed prior to the meeting, including the need to specify three alternative management strategies based on the three objectives (climate change adaptation, climate change mitigation by means of forest carbon, biodiversity conservation), and the idea of varying the realised cutting levels.The main question to the stakeholders was, what management actions they think would promote these specified objectives.
The discussion was conducted in 3 groups that were chaired by project members and included both external experts and project researchers.The groups interacted with each other during the day.When summarising the results, the researchers communicated with the experts about the possible limitations of the applied model in representing potential management actions.
The results of the meeting were summarised in a document outlining the ideas and the resulting management actions that could be used as input to the model.The document was circulated among the participants for comments.
The meeting was one in a series of stakeholder meetings organised in connection with four Finnish Strategic Research Council projects that were part of a study on knowledge and knowledge co-production (Korhonen-Kurki et al. 2022).

Figure S2. 1 .
Figure S2.1.Flowchart of the calculation model.Inputs are site fertility class and proportion of par radiation not absorbed by the tree canopy  , and outputs are ground vegetation photosynthesis, turnover and respiration.
Figure S2.4.a) Comparison of ground vegetation biomass in the observations by Kulmala et al. (2011) and the current model at sub-xeric sites.b) Comparison of the relationship between tree canopy  and combined canopy + GV  , in model and observations.The site types are OMT = grove-like, MT = herb-rich, VT = sub-xeric, CT = xeric.

Figure S2. 5 .
Figure S2.5.Comparison of predictions by PREBAS and the empirical model by Tupek et al. (2016) in Swedish NFI data.Field layer includes grasses, herbs and shrubs (kg C ha -1 ), while bottom layer includes mosses and lichens.

Figure S4. 1 .
Figure S4.1.Age class distributions evaluated from multi-source NFI (left) and plot-based NFI (right).The figure shows the total area in 20-year age classes in Finland, plus the total area of treeless sites that appear between clear cut and regeneration.

Figure
Figure S6.1a.Temporal mean values over 2021-2050 of C balance indicators under different management strategies in management-driven scenarios.Current = with current protection areas, Extended = with at least 10% set aside by all regions.

Figure
Figure S6.1.b. Left: Time development of cutting levels in management-driven scenarios from a total area of 23.4 Mha.Right: Total area-based harvest in management-driven scenarios with current and extended protected areas.

Figure S6. 2
Figure S6.2Comparison of C balance indicator variables in demand-driven scenarios under base management strategy with original and extended protected areas.Total harvests (including collected harvest residues) are shown per managed forest area (21.4 and 20.2 Mha with original and extended protected areas, respectively), other variables are shown as national averages with respect to total forest area (23.4 Mha).

Figure
Figure S6.3 a. Biodiversity indicators with current and extended (10%) protection areas in Low harvest scenarios.Miti = mitigation scenario, Prot = protection scenario.

Figure
Figure S6.3 b.Biodiversity indicators with current and extended (10%) protection areas in High harvest scenarios.
Assumptions about the mean relationship between stand age and  used in heuristic parameter estimation for percentage ground cover (Eqn S1).The left-hand side shows the mean temporal development of  , following the basic equation  =  max (1 −  .) interrupted by thinnings at times 40 and 70 yrs.The table on the right gives the maximum  ,  max , for the different site fertility classes used.
zones (hemiboreal, south, central and north boreal, northern Lapland), five site fertilities, three dominant species, and three forest development stages (regeneration, young stand, mature stand).Tonteri (

Table S3 .
1. Stand variable used to calculate HSI-values for each of the indicators.The caperceillie lekking site HSI equals to one when pine volume is high (>80m 3 /ha), spruce volume is at an intermediate level (10-20 m 3 /ha), and stem density is intermediate (600-800 stems/ha).Equals to one when the age of forest is between 20 and 60 years, proportion of deciduous trees of the total tree volume is intermediate (20-40 %) and the proportion of spruce is high (>25%).
-spotted woodpecker (Dryobates minor) HSI Equals to one when the age of forest is >60 years, proportion of deciduous trees of the total tree volume is >60%, and total basal area of trees is >15 m 2 /ha.Basal area of 15 m 2 /ha translates into total volumes well above 100 m 3 /ha in forests that are older than 60 years.Equals to one when the volume of spruce is > 175 m 3 /ha, the proportion of spruce of the total timber volume is >60%, and the volume of deciduous trees is > 15 m 3 /ha.