From Paris to Switzerland: Four Pathways to a Forest Reference Level

Introduction: Among terrestrial ecosystems, forests represent large carbon stocks threatened by changing climatic conditions, deforestation, overexploitation, and forest degradation. Close to nature forestry may help forests to continue to acting as carbon sinks by promoting their resilience against disturbances. The EU decided to carry out carbon accounting of emissions and removals from managed forests under the Paris Agreement (PA) by using a projected Forest Reference Level (FRL) based on the continuation of recent management practices. Methods: We developed four conceptual scenarios that could build the Swiss Forest Reference Level and performed simulations over 50 years using Swiss National Forest Inventory (NFI) data and the empirical forest model MASSIMO. To improve MASSIMO, we further developed a new tree species-specific model for small scale mortality that accounts for the Swiss NFI design. Then, using projected biomass and mortality from MASSIMO, carbon budgets of mineral soil, litter, and dead wood were estimated using the Yasso07 model. Results: The U-shaped mortality model performed well (AUC 0.7). Small as well as large trees had the highest mortality probabilities, reflecting both young trees dying due to self-thinning and old trees from age, pests or abiotic influences. All scenarios matched their given harvesting and growing stock targets, whereby the share of broadleaves increased in all regions of Switzerland. This resulted in decreasing biomass growth, possibly due to a species shift from typically fast growing and more shade tolerant conifers to broadleaves. The CO2-balance of the conceptual scenarios ranged from 1.06 to −3.3 Mt CO2 a–1 under Increased Harvesting and Recent Management Practices (RMP), respectively. Rotation periods are shortened under Increased Harvesting, which is an important climate adaptive management strategy, but forests were predicted to become a net carbon source. In contrast, RMP resulted in similar harvesting amounts and forests as carbon sinks, as reported in the past. Further, the RMP scenario does not involve political assumptions and reflects the idea of the CMP approach used by the EU member states, which makes it comparable to other countries. Therefore, we propose the scenario RMP as a suitable and ideal candidate for the Swiss FRL.


INTRODUCTION
Among terrestrial ecosystems, peatlands and forests constitute the main carbon stocks. However, due to changing climatic conditions, peatlands are predicted to shift from a sink to a carbon source (Loisel et al., 2021). In addition, forests may become a source due to changing climate, deforestation, overexploitation and forest degradation (Friedlingstein et al., 2019). On the other hand, close to nature forestry and climate smart forestry may help forests to continue acting as carbon sinks by promoting their resilience against disturbances (Brang et al., 2014;Gusti et al., 2020). Further, forest management can make a significant contribution to reducing CO 2 -emissions. Specifically, fossil fuels can be substituted by wood (energetic substitution) and the use of long-lived harvested wood products (HWP) can be increased at the expense of more energy-intensive materials and construction methods (material substitution; Werner et al., 2010;Braun et al., 2016). However, substitution effects are currently not considered in the reporting and accounting of the Land Use, Land-Use Change and Forestry (LULUCF) sector under the United Nations Framework Convention on Climate Change (UNFCCC) and associated agreements such as the Kyoto Protocol.
The complexity of LULUCF sector has been reason for discussions on accounting rules (e.g., Grassi et al., 2018). In the first commitment period of the Kyoto Protocol 2008-2012, the carbon accounting procedure for forests under Forest management was based on a voluntary gross-net approach taking country-specific caps into consideration (UNFCCC, 1998). A gross-net or absolute approach in this context means that countries do not include gains and losses of carbon from the LULUCF sector in their base year emissions (gross) but account for net emissions and removals from LULUCF during the commitment period (net). With this gross-net accounting approach, most managed European forests acted as carbon sinks due to legacy effects of past management and environmental changes (Nabuurs et al., 2003;Bellassen et al., 2010). In the second commitment period (CP2) of the Kyoto Protocol 2013-2020, accounting for Forest management -including harvested wood products (HWP) -became mandatory and the gross-net accounting changed to a relative or net-net accounting: only the difference between the recent forest balance and a (projected) Forest Management Reference Level (FMRL) is accounted for. The purpose of this "relative accounting with a reference level" in CP2 is to ensure that emissions or removals that occur purely from age-legacy-effects or because of effects of historical management or past natural disturbances are not accounted for Iversen et al. (2014).
With the Paris Agreement (PA), accounting rules became more flexible and in order to achieve "a balance between anthropogenic emissions by sources and removals by sinks of greenhouse gases in the second half of this century" individual countries could design their specific accounting rules based on a set of existing IPCC reporting documents (UNFCCC, 2015;Nabuurs et al., 2018). However, the PA urges its Parties to preserve and enhance existing carbon sinks including forests, integrating the "LULUCF" sector for the first time to international climate mitigation targets (Forsell et al., 2019). Moreover, given the concerns about the need to factoring out the effects of historic forest development as a result of past forest management, the UNFCCC laid out principles on situations where countries must "account" for the impact of mitigation actions (Art. 4.13 of the Paris Agreement). In response to these principles, the EU compiled rules and conditions to account for carbon credits under the LULUCF-regulation 2018/841 1 . Therefore, EU member states have to account their emissions and removals of carbon dioxide in managed forest land against a projected forest reference level (FRL), similar to the accounting of forest management under CP2 (Grassi et al., 2018). This FRL, however, must be defined more strictly and "transparently and in consideration of the business-as-usual procedures" based on the period 2000-2009 (Forsell et al., 2018). This should avoid the inclusion of increased carbon sink effects (considered as credits) due to age development of trees and forest that can be expected without any additional efforts from the forest sector (Nabuurs et al., 2018). The main difference between the FMRL in CP2 of the Kyoto Protocol and the FRL is that the FMRL projections of forest development included assumptions on future effects of policies agreed on until 2009 such as the development of harvesting levels, whereas in the FRL projections assumptions of future impacts of policies are not considered (Nabuurs et al., 2018). Since the responsibility for defining FRLs lies with the EU member states (Valade et al., 2017), member states have developed their own systems to simulate recent management practices, based on national forest inventories (NFI) and forest development models (e.g., Barreiro et al., 2016;Pilli et al., 2017;Korosuo et al., 2021).
Similar to EU member states, Switzerland decided in 2015 to account for managed forest on a net-net basis by using a reference level that includes living and dead biomass and HWP, while it excludes large-scale disturbances in forest land (see Switzerland's "intended nationally determined contribution" INDC) 2 . To do so, Switzerland uses data from its NFI in combination with the empirical forest development model MASSIMO Stadelmann, 2020) to project living biomass. Then, using MASSIMO biomass predictions, the C-cycling model Yasso07 (Tuomi et al., 2009(Tuomi et al., , 2011 is applied for estimating the carbon budget of mineral soil, litter and dead wood (Didion et al., 2014). The Swiss FRL is finally calculated as the mean annual projected CO 2 -balance of managed forest land over the period 2021-2030 by using a specific scenario.
The aim of this study is to assess how different forest management scenarios, built on the existing guidelines under UNFCCC and IPCC methodologies and consistent with the Swiss forest policy, may affect the projected mean CO 2 -balances 2021-2030. To do so, we developed four conceptual scenarios and a new model for small scale (i.e., single-tree) mortality. In this article, we addressed the following questions: (i) How can we formulate small scale mortality in an unbiased manner? (ii) How suitable are the four conceptual forest management scenarios developed based on recent forest data, forest management and forest policy objectives to represent the Swiss FRL? and (iii) How do the scenarios differ in their effect on the combined estimate of the accounting of the Swiss CO 2 balance under the Paris Agreement?

Swiss National Forest Inventory Data
Forests cover roughly one third of Switzerland's area, with higher shares in the Jura and on the Southern Alps but lower shares on the Plateau and in the Alps (Abegg et al., 2020; Figure 1) 3 . Since 1983, the Swiss NFI has repeatedly measured state and change of these forests using a systematic 1.4 × 1.4 km grid of more than 6,000 permanent plots (Fischer and Traub, 2019). Permanent plots consist of two concentric circles of 200 and 500 m 2 , in which trees are sampled if their diameter at breast height (1.3 m above the ground-DBH) exceeds the caliper threshold of 12 cm and 36 cm, respectively. In this study, we used data of the inventory campaigns NFI2 (1993/95), NFI3 (2004/06) and NFI4 (2009/17) for the fitting of the mortality model. However, for modeling the CO 2 -budgets, we only used data up to NFI3 (i.e., until 2006) in order to be in line with the KP-supplement ["not use data later than 2009"; (IPCC, 2014)] and to align with the specific requirements of the EU LULUCF-regulation 2018/841 [reference period 2000-2009(Forsell et al., 2018). However, within this 10 years reference period, only one full inventory campaign, namely NFI3 (2004NFI3 ( -2006 was conducted in Switzerland. Thus, to derive changes in forest resources we needed to use additional data of the second inventory campaign NFI2 (1993)(1994)(1995)]. For fitting of the empirical process models and for the RMP scenario, we further included data of the first inventory in order to achieve more robust models.

MASSIMO
All simulations of the living biomass presented in this study were performed with MASSIMO Stadelmann, 2020), implemented in Java SE 8. MASSIMO is individual-tree growth simulator that runs on the grid of the Swiss NFI allowing the projection of growth, ingrowth (i.e., 3 NFI result  tree recruitment), mortality and management of trees in 10year time steps, representatively. The basal area increment (BAI) model is specific to tree species groups and depends on stand basal area, basal area of trees larger than the subject tree, site index, elevation, stand age, and dominant DBH. That growth model further accounts for a growth release (i.e., increased growth) following management interventions (Thürig et al., 2005a). Models on ingrowth and non-growth above the caliper thresholds of 12 cm and 36 cm, respectively, are climate sensitive . For the simulation of ingrowth and in Yasso07, we applied a relatively warm and dry A1B scenario (CH2011, 2011Remund et al., 2016). MASSIMO differentiates between large-scale wind-throw induced mortality and individual-tree mortality. Since the last MASSIMO release Stadelmann, 2020), the individual tree mortality module has been updated to better account for stand conditions and tree species (cf. section "Mortality Module"). Forest management treatments (thinning and shelterwood cutting) can be adjusted to approximate either user-defined growing stock or harvesting targets at the level of 14 economic regions (for details cf. Stadelmann et al., 2019). While thinning is applied in all forests (except forest reserves), shelterwood cutting is only applied in even-aged high elevation forests outside of the protection forest areas (i.e., forests protecting against rockfall and avalanches). In protection forests, thinning retains stand basal area in a range defined by given recommendations (Dorren et al., 2015).
All processes were parameterized using data from the Swiss NFI (Abegg et al., 2020) ensuring accurate representation of forest dynamics in Switzerland. Model simulations are initialized with country-specific data for comparability between model results and observations and between simulated scenarios. MASSIMO and its processes have been presented beforehand (Thürig et al., 2005a(Thürig et al., ,b, 2019Stadelmann et al., 2019;Stadelmann, 2020). Since tree recruitment strongly influences forest development, a new climate sensitive ingrowth model was implemented to make MASSIMO able to cover large gradients of stand and climate conditions . However, the sub-models for growth and mortality are not climate-sensitive yet, as they are still in development. MASSIMO has been applied at regional (Temperli et al., 2017) and national scale (Stadelmann et al., 2016;Temperli et al., 2020).

Mortality Module
Two inventory intervals were used to model tree-mortality: NFI2 to NFI3 and NFI3 to NFI4. For each interval, tree-and stand-level information were extracted either at the first inventory (NFI2 for the NFI2-NFI3 interval and NFI3 for the NFI3-NFI4 interval), the second inventory (NFI3 for the NFI2-NFI3 interval; NFI4 for the NFI3-NFI4 interval), or both. For each inventory interval, only trees that were alive at the first inventory were included in the analyses. The status of trees at the second inventory was defined as dead or alive. The inclusion probability of each tree was extracted at the first inventory. It is defined based on the DBH of the tree as well as on whether the plot could be entirely measured. Tree species were classified into 12 groups based on the most represented tree species or genera in Switzerland: spruce (Picea abies), fir (Abies alba), larch (Larix sp.), pine (Pinus sylvestris, Pinus nigra, Pinus strobus, Pinus mugo arborea), stone pine (Pinus cembra), beech (Fagus sylvatica), oak (Quercus sp.), ash (Fraxinus sp.), maple (Acer sp.), and chestnut (Castanea sativa); others were classified as other conifers or other broadleaves. The diameter at breast height (DBH, 1.3 m above ground) was extracted from the first inventory. Finally, stand basal area (in m 2 ha −1 ) was used as a plot-level variable of the first inventory. The normalized number of vegetation years between each inventory (number of vegetation years minus 10) and the production region of each plot were also noted.
The tree mortality analyses were performed using R version 3.5.1 (R Core Team, 2018). Models were built using the survey R package (Lumley, 2020) that is able to account for complex survey samples. The svydesign function was used to specify the data design including plot number as the cluster id and trees' inclusion probabilities as sampling weights. Binomial regression models were then fitted with the svyglm function using tree status at the second inventory (dead/alive) as the response variable and DBH (second-degree polynomial term), species group, production region, stand basal area, and normalized vegetation years as explanatory variables. An interaction between DBH and species group was also added to account for size differences between species. The goodness of fit of the model could be assessed using AUC.

Yasso07
Yasso07 (Tuomi et al., 2009(Tuomi et al., , 2011 is a model of C cycling in mineral soil, litter, and dead wood. The model was developed for general application with low parameter and input requirements. For estimating stocks of organic C in mineral soil up to a depth of ca. 100 cm and the temporal dynamics of the C stocks, Yasso07 requires information on C inputs from dead organic matter components (i.e., litter and dead wood) and climate (temperature, temperature amplitude, and precipitation). Decomposition of the different components of C inputs is modeled based on their chemical composition, size of woody parts and climate (Tuomi et al., 2009(Tuomi et al., , 2011. Decomposition rates of C that is either insoluble (N), soluble in ethanol (E), in water (W), or in acid (A), and flow rates of C between these four compounds, to a more stable humus compartment (H) and out of the soil were derived from a global data set (Tuomi et al., 2009(Tuomi et al., , 2011. Parameters defining C decomposition rates and C cycling between carbon compartments are obtained probabilistically using Markov Chain Monte Carlo sampling and are fitted based on data from litterbag studies (Tuomi et al., 2009) and woody litter decomposition experiments (Tuomi et al., 2011). The approach allows calculating the maximum a posteriori point estimate and confidence sets for the model parameter vector, which consists of a unique combination of 24 correlated parameters. At the time of this study, three independent parameter sets, i.e., vectors of 24 parameters have been developed and published (Tuomi et al., 2009(Tuomi et al., , 2011Rantakari et al., 2012). Didion et al. (2014) examined the validity of the Yasso07 model in Swiss forests based on data on observed decomposition of (i) foliage and fine root litter from sites along a climatic and elevation gradient and (ii) of dead trees from plots of the Swiss National Forest Inventory. The study analyzed, among other, the accuracy of Yasso07 for reproducing observed C decomposition in litter and dead wood in Swiss forests. The best agreement of Yasso07 results with observed C decomposition was obtained with the parameter presented in Rantakari et al. (2012), which was obtained by restricting the global dataset on dead wood, litter, and soil decomposition set (Tuomi et al., 2009(Tuomi et al., , 2011 to European sites. The Yasso07 model is used to obtain estimates of C stocks in dead wood and litter and of C stock changes in dead wood, litter, and mineral soil for reporting in Switzerland's annual National Inventory Reports to the UNFCCC. Austria, Finland and Norway also apply the Yasso07 model for their reports to the UNFCCC. Yasso07 simulations for this study are based on the Yasso07 Fortran source code of the Yasso07 release 1.0.1 4 .

Testing Four Conceptual Scenarios for the Swiss Reference Level for Forest Land
Together with representatives of the Swiss Federal Office for the Environment (FOEN) we developed four conceptual scenarios to be evaluated for establishing a reference level for Forest Land and ran them from NFI3 (2004/06) for 5 time-steps of 10 years (Table 1). Thus, the simulations started with the NFI3 conditions with a growing stock of 360 m 3 ha −1 which corresponds to a living biomass stock of 121 TC ha −1 (data per NFI production region given in Supplementary Table 2). The "Forest Policy (FPOL)" scenario, however, ran for 30 years only as scenario assumptions could not be met in longer simulations. Thereby, we ran each scenario with 30 replicates to account for stochastic differences between each run. Except for the "Recent Management Practices (RMP)" scenario, we simulated all scenarios with standard management algorithms and a scenario approximation routine to fulfill given growing stock or harvesting, respectively, per economic region (for details on the scenario adaption routine, cf. Stadelmann et al., 2019). The RMP scenario did not use this routine nor the standard management algorithms, as it has been developed to be selfapproximating via its single-tree selection model (Stadelmann et al., in prep).
As default settings in MASSIMO, we defined heavy storms to occur every 15 years as it has been observed in the past (Pfister, 1999;Thürig et al., 2005b) and 80% of storm damaged timber was salvage logged to avoid bark beetle outbreaks (cf. Stadelmann et al., 2013). However, in the RMP scenario we shortened the storm interval to 10 years as observed between NFI1 and NFI3. This is in line with Usbeck et al. (2010) that found increasing storm damages over time. Further, we activated protection forest management as presented in Stadelmann (2020) and respected existing forest reserves.
The amount of mortality simulated with the new mortality module agrees well with the empirically measured mortality (Abegg et al., 2020) 5 . However, to meet given targets of the background mortality defined for Swiss forest under management under the Kyoto Protocol (for details on the estimation of background mortality cf. FOEN, 2016), the empirical mortality had to be increased slightly for FPOL, "Increased Harvesting (IH), " and "Retain Growing Stock (RGS)" by additionally simulating harvesting losses (i.e., we attributed part of the harvested tree compartments to remain in the forest as harvesting losses). The background mortality, estimated as 16.39%, is defined as a share on all losses (harvesting of wood + mortality) between the years 1990 and 2009 (cf. FOEN, 2016) that include two heavy storm events -Vivian/Wiebke in 1990 and Lothar in 1999 -both followed by large scale bark beetle outbreaks. While the storm damages by Lothar were treated as outlier, the mean of all other damages build the background mortality of 16.39%, which must be reached as input in Yasso07. Besides simulated mortality, stumps and roots of harvested trees entered this pool. But additional deadwood was needed to reach the given target for all scenarios but RMP. By assumptions on harvesting losses, we could ensure that exactly the predefined amount of background mortality is used as input in the model Yasso07. Harvesting losses differed for conifers and broadleaves as well as for stem, thick (>7 cm in diameter) and small (<7 cm in diameter) branches ( Table 1).
The initial state of the dead wood, litter, and mineral soil pools in 2006 (Supplementary Table 3) was taken from Yasso07 simulations for reporting in the Swiss Greenhouse Gas Inventory report 2020 (NIR; FOEN, 2020) with annual litter and dead wood data derived from measured NFI data. The plausibility of simulated C stocks and stock changes in these pools is thus continuously evaluated (see section 3.1 and 3.2 in Didion, 2019 for the data submitted to NIR 2020). In this study, the simulated development of those pools is driven by annual input of dead wood (incl. stemwood and large branchwood that are not removed from the plot and all coarse-root and stump biomass) and litter (incl. fine woody litter from small branchwood and non-woody litter like foliage and fine roots) derived from simulated MASSIMO results for each NFI plot following the methodology presented in Didion and Zell (2019). The Yasso07 simulations until 2006 were driven by temperature and precipitation data extracted from the gridded data prepared by MeteoSwiss (2016a,b). From 2007, we used simulations from MASSIMO to derive data on growth and mortality, as well as litter and dead wood remaining in the forest as inputs for the Yasso07 simulation. To account for model uncertainty, the Yasso07 simulations were based on 10 randomly selected simulated initial states for 2006 (out of 100 replicates prepared for the reporting to the NIR, cf. Didion and Thürig, 2019) and 10 randomly selected MASSIMO replicates to obtain the time series of dead wood and litter production starting 2007.

Scenario
Specifications and scenario settings within MASSIMO

Recent management practices (RMP)
Instead of targeting given harvesting or growing stocks, RMP aims at empirically adapting observed harvesting from NFI1 (1983/85) to NFI3 (2004/06). To this end, we used single-tree selection probabilities that depend on DBH, tree species group, production region, stem count, stand basal area, elevation, slope, site index, increment, owner categories (public vs. private), harvesting infrastructure and occurrence of storm damages. Predictions of RMP have been validated against independent NFI4 observations (cf. Supplementary Figure 4; Stadelmann et al., in prep) This scenario did not include special forest management treatments such as protection forest management or forest reserves, as this would lead to an underestimation of harvesting. Further, salvage logging was increased by 3% from the default ending with 83% of storm-damaged timber that will be harvested in this scenario Except for small branchwood (100% remains in forest) beside the simulated mortality no additional harvesting losses were needed to reach a background mortality of 16.39% RMP reflects the concept of "continuation of forest management practices"(as documented in the period 2000-2009) (CMP; as it is defined by the LULUCF-regulation 2018/841), but it is not identical since not all requirements are fulfilled Forest policy (FPOL) Increased merchantable timber harvesting from 6.69 Mio. m 3 a −1 observed in NFI3 (2004/06) to 8.14 Mio. m 3 a −1 in 2036 as envisaged by the Swiss forest policy (FOEN, 2013). Thereby, harvesting favors conifers against broadleaves and take place mainly in regions where harvesting costs are low (i.e., good harvesting infrastructure). Due to increased shares of broadleaves on the Plateau and a decrease in growing stock in that region, harvesting of 8.14 Mio. m 3 a −1 is not possible after 2036, when still favoring conifers and regions with lower harvesting costs.

Projected CO 2 -Balance
In the context of climate reporting, Switzerland announced in its "Nationally determined contribution" (NDC) 6 to account for its forest sector on a net-net basis, i.e., relative accounting against a FRL. The FRL is a projected mean CO 2 -balance of all forest pools for the compliance period 2021-2030 based on a conceptual scenario. The conceptual scenarios for the Swiss FRL were evaluated by calculating the total CO 2 -balances per scenario ( Table 2). To this end, we averaged simulation results of the decades 2017-2026 and 2027-2036 by weighted means for the compliance period 2021-2030 for all modeled carbon pools of the Swiss forest sector. The CO 2 -balance is the sum of the carbon stock changes in living biomass, soil carbon, litter, dead wood and "harvested wood products (HWP)" and emissions from controlled fires. HWP calculations followed the methodology described in FOEN (2021a) and based on the methodology internationally agreed on and applied during the second commitment period of the Kyoto Protocol (IPCC, 2014).

Mortality Model
The mortality model fitted a positive relationship between tree mortality probability and stand basal area, used as a proxy for competition. The relationship to DBH had a U-shaped form for most species groups, where small and large trees had a larger mortality probability than average-sized trees (Figure 2)

Scenario Simulations
All scenarios (but RMP that has no defined target) matched their defined harvesting and growing stock targets ( Supplementary  Figure 1), whereby the shares of broadleaves increased in all production regions and for all scenarios. Consequently, we found the same pattern for living biomass (Figure 3). RMP resulted in increasing biomass for broadleaves, but also for conifers in the Alps and the Southern Alps, where harvesting is lower than growth. In these regions, FPOL followed the same pattern as RMP, while the other regions -that have better harvesting infrastructure -did not. On the Plateau, where harvesting is similar to IH, biomass and growing stock strongly decreased. With IH, biomass and growing stock also decreased in all other regions except for the southern Alps that currently has relatively low growing stocks (Supplementary Figure 1), due to legacy effects for past coppice forest management (i.e., chestnut coppice forests). IH is the only scenario that showed similar trends for conifers and broadleaves in all regions. While RGS resulted in constant growing stocks, biomass slightly increased with mounting shares of broadleaves but the biomass of conifers slightly decreased. Gross growth of biomass decreased for all scenarios and production regions for total timber, as well as conifers and partly also for broadleaves ( Figure 4A and Supplementary Figure 2), 6 Avl. Online at: https://tinyurl.com/SwissNDC with decreases being most pronounced in IH that resulted in lowest biomass stocks. Similarly, total losses of biomass were highest for IH ( Figure 4B). These trends were related to the given scenario targets. For RMP, however, we found a steep increase in losses of broadleaves biomass, mainly in Jura and on the Plateau (cf. Supplementary Figure 3).

Yasso07 Simulations
To assess the relative impact of the different scenarios on the changes of carbon stocks in DOM (dead organic matter, i.e., the sum of litter, dead wood) and mineral soil pools, simulated C stocks over time are presented relative to the RMP scenario ( Figure 5). Because the initial state of the DOM and soil C pools were based on litter and dead wood input data derived from NFI measurements (cf. section "Testing Four Conceptual Scenarios for the Swiss Reference Level for Forest Land"), the change in inputs for the Yasso07 model from observed NFI data to predicted MASSIMO data was primarily responsible for the dynamics of DOM and soil carbon stocks in the first decades of the modeling. Thereby, the application of a relative level of mortality (16.39% background mortality) led to more pronounced absolute mortality rates, especially for scenarios with increased harvesting rates compared to RMP. But mortality of stemwood also increased for RMP, since we simulated shorter time intervals (i.e., 10 instead of 15 years) between heavy storm events. Carbon stocks in DOM and soil changed inversely to living biomass stocks, as 16.39% of the total losses (harvest and mortality of biomass) entered the pool of litter and dead wood, and, in the longer term, mineral soil. At the national level, living biomass continuously increased under RMP, while this increase was less pronounced for FPOL and RGS and decreased for IH until 2036 (Figure 3). This contributed to the increase in DOM and soil carbon stocks for all scenarios compared to RMP (Figure 5) for the first 15-30 years, with regional differences. DOM and soil carbon stocks then decreased compared to RMP for the remaining simulation period. The only region with a different trend from the national level was the Plateau, where the increase of carbon stocks in these pools compared to RMP was higher for all scenarios, but decreased between 2036 and 2046 for IH, ending with similar DOM and soil carbon stocks at the end of the simulation period. For the other regions, non-living carbon developed similarly to the national scale.

Mean CO 2 -Budget 2021-2030
Similarly to the observed CO 2 -budget of managed forest land during the period 2013-2019 reported by the Switzerland's Greenhouse Gas Inventory (NIR 2021: FOEN, 2021bThürig et al., 2021), the scenarios RMP, RGS, and FPOL resulted in a CO 2 -sink for the next compliance period 2021-2030, while IH resulted in a CO 2 -source. The highest sink occurred with RMP, resulting in harvesting amounts similar to those reported over the period 2013-2019 (NIR21). Harvesting rates obtained with the scenarios FPOL, RGS, and IH are remarkably higher than those of RMP. With increasing harvesting rates, the total CO 2 -balances declined or even turned into a net source (IH). Harvest (merchantable timber in Mio. m 3 a −1 ) 6.5* 6.6 7.9 7.3 8.4 The annually harvested timber volume is also presented. *Observed in the change inventory NFI3-NFI4: doi: 10.21258/1475230. The CO 2 -balance (bold) is the sum of living biomass, combined DOM and soil, HWP, controlled fires.
FIGURE 2 | Predicted mortality probabilities per tree species group along gradients of stand basal area and DBH per production region of Switzerland. Other variables were kept constant at their mean species-specific values.

DISCUSSION
The new mortality model performed reasonably well, not only when assessing the goodness-of-fit of the model directly, but also within the modeling framework MASSIMO where growth and losses were within the range of NFI observations . Under RMP, where future management (harvesting amounts) remains within the range of past and recent management, we did not need to simulate additional harvesting losses to reach a given background level of mortality as defined by FOEN (2016) in accordance with the accounting rules given in the second commitment period of the Kyoto Protocol. This background level of mortality includes mortality of single trees and large scale disturbances (see section "Testing Four Conceptual Scenarios for the Swiss Reference Level for Forest Land") that have been observed between 1990 and 2010. Since the background mortality is defined in relation to the total losses (harvest + mortality) any scenario with increased harvesting FIGURE 3 | Switzerland's living biomass stock (above and below ground) for all Switzerland and its production regions. Results of conifers (middle row) and broadleaves (lower row) sum up to the total biomass stock in the upper row. Simulation results include mean (thick line) and Standard Error of 30 replicates (shaded area). The lines between each end of a time-step (e.g., between 2026 and 2036) are linearly interpolated.
inevitably results in increased absolute mortality. Thus, to reach the 16.39% background mortality predefined by FOEN (2016), we needed to simulate harvesting losses for the scenarios IH, FPOL, and RGS that account for stumps, thick and small branches.
In RMP, where harvesting is similar to recent observations, the empiric mortality is well predicted. Thus, the new species-specific mortality model allows minimization of additional artificial mortality (i.e., harvesting losses). Moreover, it plausibly reflected the effects of site and stand conditions. Increasing basal area, for example, led to a well-known density-dependent mortality (Larson et al., 2015). The U-shaped mortality response to tree size reflects both young trees dying due to self-thinning (e.g., Westoby, 1984) as well as old trees dying for various reasons (age, pests, or abiotic influences). However, regional or national increases in harvesting intensity, as represented by the IH scenario, reduces stand density and the number of large trees, resulting in turn in decreasing mortality rates (Stadelmann et al., 2016;Temperli et al., 2017).
In all scenarios, the up to 2056 simulated share of broadleaved trees increased, while the gross growth of volume and biomass decreased. This can result from a shift in tree species from fast growing and typically more shade-tolerant conifers to broadleaves, a pattern that has been observed since the first NFI, and which is expected to increase due to climate warming. Broadleaves typically need more stand space, resulting in lower potential stem densities. Further, all scenarios except RMP resulted in an increase of harvesting that led in the shortterm to a rapid decrease in growth. The species shift from conifers to broadleaves was most pronounced in RMP, with the consequence that harvesting and therefore biomass losses (i.e., harvest and mortality) of broadleaves increased steeply. However, if harvesting costs decrease and the demand for wood fuel increases (c.f. Thees et al., 2020), the removal of thick branchesespecially from broadleaves -will need to be intensified. This would result in lower harvesting losses that would, when summed up with background mortality, not reach the predefined 16.39% threshold. Consequently, conceptual scenarios could be further developed to include a varying demand in wood fuel resulting in varying harvesting losses and varying background mortality.
Our scenarios strongly influenced the amount of harvesting losses needed to reach the predefined background mortality and volume of branchwood that entered the Yasso07-pools for carbon change simulations. Indeed, scenarios other than RMP were based on higher harvesting rates, which resulted in an increase in the combined DOM and soil carbon stock. Intensifying management led to a decrease in stem density compared to RMP and therewith, litter volumes from living trees decreased. Since climate is an important driver of processes such as forest regeneration, growth and mortality, as well as decomposition of DOM, it would be preferable to use several climate scenarios to obtain a corresponding range of carbon sequestration trajectories. However, since so far, only the regeneration module in MASSIMO is climate sensitive, we decided to apply only one climate scenario for the simulation of ingrowth and Yasso07, assuming that growth and mortality are not affected. Compared to current climate, the use of this scenario likely resulted in an increased C sink of the DOM and soil pools due to more optimal conditions (i.e., warmer temperatures with FIGURE 5 | Carbon stock relative to the recent management practices (RMP) scenario. Values above 1 indicate higher combined dead organic matter (DOM) and soil carbon stocks than in RMP and vice versa. still sufficient precipitation in the majority of Swiss forests) lead to increased tree growth but also to a faster decomposition of the DOM. Similar forest growth models applied in Europe are not climate sensitive (Barreiro et al., 2016;Vauhkonen et al., 2019) with the assumption that on relatively short time scales of 50 years management effects are more important than climate ones (Perin et al., 2021). For future improvements of MASSIMO, however, the effects of climate change on forest development should be accounted for, as they are becoming increasingly apparent (e.g., Schelhaas et al., 2015;Reyer et al., 2017;Seidl et al., 2017;Rohner et al., 2021). Consequently, including climate-sensitivity into increment and mortality models will be inevitable to allow for analyses of scenarios assuming changing climatic conditions (Matala et al., 2005;Kindermann, 2010;Rohner et al., 2018;Schelhaas et al., 2018).
Swiss forests are predicted to become a net source of carbon under scenarios with fast reductions in growing stock in the short-term (e.g., FPOL and IH), while they may lead to increasing gross growth in the long-term (Stadelmann et al., 2016). These growing stock reductions are obtained by shortening the rotation periods, which is a climate adaptive management strategy (Brang et al., 2014). Additionally, increased harvesting allows for material and energetic substitution effects that contribute in mitigating climate change (Mehr et al., 2018;Churkina et al., 2020;Jonsson et al., 2021). However, according to reporting requirements under the international climate agreements such as the Kyoto Protocol and the Paris Agreement, substitution effects of the forestry sector cannot be directly accounted for, although they indirectly are through the resulting lower emissions in other sectors such as energy. Nevertheless, these substitution effects have a considerable impact on the CO 2 -balance of forest management including harvested wood products (Werner et al., 2010). Substitution effects might even overcome harvestinginduced carbon emissions in some cases (Braun et al., 2016).
The four presented scenarios were evaluated for their plausibility to represent the Swiss FRL for the accounting of emissions and removals on managed forest land under the Paris Agreement. While RGS is a good base scenario to assess changes in growth and losses, it provides no path-breaking information, nor a direct reference to a realistic future development of Swiss forests. Thus, it does not suit as Swiss reference level. While RMP is fully empiric (i.e., the single-tree harvesting model depends on observed management only), IH, RGS and FPOL are based on transparent and reproducible management assumptions such as thinning frequency and intensity and rotation periods. FPOL includes harvesting targets of the Swiss Forest Policy, which envisages a distinct increase in harvesting rates in order to obtain stable, resilient, multifunctional and climate-adapted forests in the future. However, in this scenario a strong priority setting of harvesting costs was implemented, i.e., forests with very high harvesting costs were not harvested. Therefore, we simulated harvesting mainly on the Plateau in the Jura and at the pre-Alps where harvesting infrastructure is better than in the Alps and the Southern Alps, resulting in a decrease of available timber that prevented simulating the scenario for more than 30 years. In fact, current harvesting rates are significantly lower than simulated in FPOL and without a substantial increase in timber prices, this scenario is not realistic in the near future. This limitation also emerges for the IH scenario, where harvesting is even higher than in FPOL but more evenly distributed across Switzerland. Consequently, IH and FPOL both include policy assumptions that lead to an overestimation of short-term harvesting. These policies could seem arbitrary and unrealistic, and overestimate carbon emissions in the reference level. Consequently, choosing these scenarios to be used for a Swiss reference level would clearly result in unwarranted credits and thus jeopardize the environmental integrity of these scenarios.
Recent Management Practices, eventually, is the most impartial scenario as it is empirically based on country-specific NFI observations. We consider RMP as the most environmentally sound scenario, since only CO 2 -effects directly related to changes in management practices are accounted for. Additionally, the EU has chosen to base its accounting of managed forest land on a FRL relying on a scenario called "continuation of forest management practices (as documented in the period 2000-2009) (CMP), " on which RMP is conceptually based. Choosing RMP therefore makes the Swiss results comparable to those of neighboring countries. In December 2020, Switzerland submitted its carbonemission reduction pledge under the PA (cf. Switzerland's full NDC) 7 and announced to account for managed forest land on a net-net basis, i.e., relative accounting against an FRL with an underlying scenario excluding assumptions on the impact of future forest policies. Since the Swiss RMP is empirically driven from NFI data covering the period 1983-2006, this approach is very similar (but not identical, since years outside [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009] are used) to the approach applied by the EU. To conclude, the RMP scenario for Switzerland is reliably based on observed data, comparable but not identical to the EU approach and results in plausible estimates, making it the ideal candidate for the Swiss FRL.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
GS developed the conceptual scenarios to be tested for the Swiss Reference Level in collaboration with NR and ET, ran all simulations with MASSIMO and Yasso07, and led the writing of the manuscript with contributions by all other authors. JP developed the new mortality module together with GS and ET. NR calculated the CO 2 -balances. MD supported all work related to Yasso07. All authors contributed to the article and approved the submitted version.

FUNDING
This study was performed with the support of Swiss Federal Office for the Environment.