The E ﬀ ect of Harvesting on National Forest Carbon Sinks up to 2050 Simulated by the CBM-CFS3 Model: A Case Study from Slovenia

: With the advent of global warming, forests are becoming an increasingly important carbon sink that can mitigate the negative e ﬀ ects of climate change. An understanding of the carbon dynamics of forests is, therefore, crucial to implement appropriate forest management strategies and to meet the expectations of the Paris Agreement with respect to international reporting schemes. One of the most frequently used models for simulating the dynamics of carbon stocks in forests is the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3). We applied this model in our study to evaluate the e ﬀ ects of harvesting on the carbon sink dynamics in Slovenian forests. Five harvesting scenarios were deﬁned: (1) business as usual (BAU), (2) harvesting in line with current forest management plans (PLAN), (3) more frequent natural hazards (HAZ), (4) high harvest (HH) and (5) low harvest (LH). The simulated forest carbon dynamics revealed important di ﬀ erences between the harvesting scenarios. Relative to the base year of 2014, by 2050 the carbon stock in above-ground biomass is projected to increase by 28.4% (LH), 19% (BAU), 10% (PLAN), 6.5% (HAZ) and 1.2% (HH). Slovenian forests can be expected to be a carbon sink until harvesting exceeds approximately 9 million m 3 annually, which is close to the calculated total annual volume increase. Our results are also important in terms of Forest Reference Levels (FRL), which will take place in European Union (EU) member states in the period 2021–2025. For Slovenia, the FRL was set to − 3270.2 Gg CO 2 eq / year, meaning that the total timber harvested should not exceed 6 million m 3 annually.


Introduction
Global forests represent one of the most important terrestrial carbon stocks [1] and play an important role in global carbon cycles, particularly due to their carbon sequestration capacity and positive influence on water balance and temperature regulation [2]. Global carbon stocks in forests are estimated at 861 ± 66 Pg C, of which 44% represents forest soils, 42% forest biomass, 8% dead wood biomass and 5% litterfall. Globally, over a half of the carbon stocks in forests can be found in tropical forests (55%), followed by boreal (32%) and temperate forests (14%) [3].
Greenhouse gas emissions and removals by forests are reported annually by countries committed to the United Nations Framework Convention on Climate Change (UNFCCC) [4], in the framework of the land use, land-use change and forestry (LULUCF) sector. This is one of five economic sectors, with the other four being energy, industrial processes and product use, agriculture and waste. The European Union (EU) included the LULUCF sector in its energy and climate policy for the reduction of emissions forest types and landscape spatial entities. The simulation results are provided on an annual basis, separately for all five carbon pools [23]. Advanced options for displaying the results enable analyses of the carbon transitions between forest carbon pools, the atmosphere and harvested wood products.
The structure of the CBM-CFS3 model includes three main steps [24]: (1) a pre-processor program which prepares the inventory database and generates the pool for dead biomass, (2) a pre-processor program which calculates the carbon stock for individual pools and localities on an annual basis during the simulation, and (3) an archive index database which includes model parameters and connects them to the input data and simulation results. The changes in below-ground biomass are calculated using the methodology presented by Li, et al. [31]. The archive index database was customised for the EU member states and includes climate information which affects the annual decay rate of dead organic matter, the ecological parameters for individual bio-geographical regions of the EU and the set of volume-to-biomass conversion coefficients for European tree species which are important for the correct conversion of merchantable volume into biomass and foliage components.
The model input data are represented by seven independent matrices. The Classifiers and Values matrix defines the number of model units and the classifiers which define them. In our case, the classifiers were defined by individual forest types and further divided by the share of conifers/broadleaves. The Age Classes matrix defines the range of existing age classes and the degree of transition from one age class to another. The initial state of forests is described in the Inventory matrix, which includes the areas of individual types of forests, further divided into mixture and age classes. The development of woody biomass in forests by age for individual types of forests is determined in yield curves included in the Growth and Yield matrix. Growing stock is defined as the gross merchantable volume of biomass, which includes the volume inside the bark of the main stem, excluding tops and stumps, but including defective and decayed wood of trees or stands [32]. With the Transition Rules matrix, we can define the rules on how individual forest types transit from one to another. For example, a large-area disruption which affects an adult spruce stand can transition to another forest type with a higher share of broadleaves and thus more resilient to natural hazards. All natural and man-made disturbances, including harvesting, are defined in the Disturbance Events matrix. Disturbances can be expressed in absolute spatial units or as an area proportion on which the disturbance can occur. The Disturbance Types matrix precisely defines disturbances and their influence, and later connects them to the disturbance parameters set out in the archive index database.
The inventory database, yield curves, disturbance events and harvest schedule are key data inputs for the simulation in the CBM-CFS3. For each year of simulation, disturbance must be related to a specific forest type and its age class and connected to one of the available disturbance types from the archive index database. There are more than 100 disturbance types available, ranging from forest fires to bark beetle outbreaks and different types of harvesting. In addition, the user also has the possibility of determining their own type of disturbance and the intensity of its influence on forest carbon pools.
The CBM-CFS3 was first developed for use in even-aged forests where the age of trees is known. However, when yield curves realistically describe the development of individual stands of different species, the model can also be used in mixed uneven-aged forests which are typical for Slovenia [32]. Two basic data sources were used to run the model, both provided by the Slovenia Forest Service (SFS): (1) the Forest Compartment Database from 2014 and (2) the Permanent Sample Plot (PSP) Database. The Forest Compartment Database consists of approximately 53,000 forest compartments with a mean size of 22 ha, which are permanent forest planning categories and cover the information on all forests in the country. For each compartment, various forest attributes are available, such as forest area, forest type, growing stock and tree species composition. Using data on forest type, similar forest compartments were aggregated and represented the initial state of the forests in 2014. The PSPs are part of the control sampling method in Slovenia [33]. Each plot (500 m 2 each) is remeasured every 10 years, and common tree attributes are surveyed for each tree on the plot, such as location, tree species, DBH, height of selected trees and status between consecutive inventories (e.g., unchanged, harvested, died, ingrowth).

Preparation of Inventory Data
In our study, the model classifiers were defined as a combination of forest type and tree species mixture. Forest type was determined according to the typology of Slovenian forest sites [34], see Table 1. Tree species mixture was represented by three categories based on the prevailing tree species: (1) if the total growing stock of broadleaves within a compartment was greater than or equal to 75%, the forest compartment was classified as broadleaved forests; (2) if the total growing stock of conifers was greater than or equal to 75%, the forest compartment was classified as coniferous forests; and (3) in all other cases, the forest compartment was classified as mixed forests, which was the prevailing category in all forest types (Table 1). Finally, 43 actual model classifiers were defined.  [34], assigned climate unit from the archive index database, total share of areas among different forest types and the share of broadleaved (BRD), coniferous (CON) and mixed (MIX) stands within each forest type.

Short
Forest Forest types and model classifiers were then linked to the appropriate climate unit of the archive index database ( Figure 1; Table 1). Thus, we provided a suitable set of ecological parameters which are necessary and affect the carbon cycle, such as the decay rate of dead organic matter, the transfer of carbon between carbon pools and litterfall characteristics. The majority of forest types were classified in the following climate units: Slovenia-CLU55 (central and southern Slovenia), Slovenia-CLU45 (Savinja and Styria regions) and Slovenia-CLU54 (Montane and altimontane forests at higher altitudes in the Alps and Pohorje).

Development of Growth Curves and Determination of Age
To initialise and run the CBM-CFS3, at least the following input data are required for each classifier: areas by age classes and yield curves. Areas were calculated from the Forest Compartment Database and yield curves from the PSP database.
Since in Slovenia we do not classify forest stands into age classes when planning forest management, and moreover we manage a significant area of uneven-aged stands, we had to define the (approximate) age of each stand on PSPs, first. We did that following the ensuing methodology. The dominant diameter DBHdom per plot (i.e., mean DBH of the 100 thickest trees per hectare) was calculated for the two most recent measurements (DBHdom.0 and DBHdom.1); both were then classified into DBH classes of 5 cm. For each PSP, the average diameter increment of dominant trees (IDBH.dom) was calculated as the difference between DBHdom.1 and DBHdom.0. Based on the IDBH.dom, the transition periods that dominant trees needed to overgrow the observed DBH class were calculated for each forest type. By summing these transition periods from that of the lowest DBH class to that of the observed DBH class on a plot, we estimated the age of a stand on each PSP. When the age of a stand was calculated, we assumed that dominant trees need 20 years to achieve the DBH measurement threshold of 10 cm, regardless of the forest type they grow in. With all that data, we could finally classify stands on PSPs into age classes spanning 20 years each. We classified 10 age classes, each ranging 20 years, but the oldest age class (AGEID09) included all stands with an age greater than 181 years ( Table 2).

Development of Growth Curves and Determination of Age
To initialise and run the CBM-CFS3, at least the following input data are required for each classifier: areas by age classes and yield curves. Areas were calculated from the Forest Compartment Database and yield curves from the PSP database.
Since in Slovenia we do not classify forest stands into age classes when planning forest management, and moreover we manage a significant area of uneven-aged stands, we had to define the (approximate) age of each stand on PSPs, first. We did that following the ensuing methodology. The dominant diameter DBH dom per plot (i.e., mean DBH of the 100 thickest trees per hectare) was calculated for the two most recent measurements (DBH dom.0 and DBH dom.1 ); both were then classified into DBH classes of 5 cm. For each PSP, the average diameter increment of dominant trees (I DBH.dom ) was calculated as the difference between DBH dom.1 and DBH dom.0 . Based on the I DBH.dom , the transition periods that dominant trees needed to overgrow the observed DBH class were calculated for each forest type. By summing these transition periods from that of the lowest DBH class to that of the observed DBH class on a plot, we estimated the age of a stand on each PSP. When the age of a stand was calculated, we assumed that dominant trees need 20 years to achieve the DBH measurement threshold of 10 cm, regardless of the forest type they grow in. With all that data, we could finally classify stands on PSPs into age classes spanning 20 years each. We classified 10 age classes, each ranging 20 years, but the oldest age class (AGEID09) included all stands with an age greater than 181 years ( Table 2). The forest area was calculated for each forest type by summing the areas of the compartments of that forest type. The total area of each forest type was then proportionally divided into mixture types and age classes according to the share of PSPs of each combination of mixture type and age class. Since PSPs are generally not tallied in young forests, the youngest age class (0-20 years) included the areas of young forests obtained from the stand map of the SFS.
For each classifier, the yield curve was calculated using Equation (1), where GS represents growing stock, AGE is the middle age of an age class, and a, b, c represent model parameters. For all forest types, yield curves for 43 classifiers were developed ( Figure 2) according to the method described above.

Harvesting Scenarios
The time-related carbon sink dynamics in Slovenian forests was analysed according to five hypothetical forest management scenarios ( Figure 3): (1) business-as-usual (BAU), (2) harvesting in line with current forest management plans (PLAN), (3) more frequent natural hazards (HAZ), (4) high harvest (HH) and (5) low harvest (LH). All scenarios are based exclusively on the assumption of future harvest intensity and amount of harvested biomass, since the latter most significantly  Table 1.

Harvesting Scenarios
The time-related carbon sink dynamics in Slovenian forests was analysed according to five hypothetical forest management scenarios ( Figure 3): (1) business-as-usual (BAU), (2) harvesting in line with current forest management plans (PLAN), (3) more frequent natural hazards (HAZ), (4) high harvest (HH) and (5) low harvest (LH). All scenarios are based exclusively on the assumption of future harvest intensity and amount of harvested biomass, since the latter most significantly determines carbon stock dynamics. The BAU scenario is based on the realised timber harvest in the past 10 years. Importantly, the BAU scenario does not account for additional harvesting due to natural disturbances. Based on the amount of timber harvested in the period 2004-2019 [9] and omitting data from years in which natural disturbance events occurred (2014-2016, 2018), the harvest intensity was regressed as a function of time (R 2 = 0.789). Prior to 2030, we used the calculated function to estimate annual timber harvest, while it was assumed that the harvest amount would increase by 1% per year in the period 2030-2050. Thus, we assumed a progressive increase in harvesting, which is expected to exceed 8 million m 3 in 2050. When determining the PLAN scenario, we followed the recommendations from forest management plans where the annual allowable cut is defined at a rate of 75-90% of the annual stand volume increment [9]. The annual harvested amount in the PLAN scenario was determined by fitting the negative exponential function, with the dependent variable being the increase in the rate of planned harvest compared to the planned harvest in the previous year (R 2 = 0.358). To calculate the adjusted values, we used the planned harvest in 2004 of 4,162,662 m 3 [35] as a starting point. Until 2035, the degree of increase of the harvested biomass took place pursuant to the calculated function; after 2035, a constant level of increase of 0.04% per year was used. The HAZ scenario was also based on the PLAN scenario and assumed a 3% increase in harvesting levels on an annual basis, as well as four extraordinary natural disturbances, appearing in an interval of approximately 10 years and assuming an increase in harvest intensity over the following three years, which is usually a consequence of bark beetle outbreaks [12]. The HH and LH scenarios are primarily based on the assumption of the annual allowable cut of the PLAN scenario, but the intensity of harvesting in the HH scenario was increased by approximately 30% and in LH reduced by 40%. With the selected increase and reduction, we wanted to obtain a higher range of annual harvesting levels, which could be potentially beneficial for the evaluation of harvesting levels on carbon stocks. The BAU scenario is based on the realised timber harvest in the past 10 years. Importantly, the BAU scenario does not account for additional harvesting due to natural disturbances. Based on the amount of timber harvested in the period 2004-2019 [9] and omitting data from years in which natural disturbance events occurred (2014-2016, 2018), the harvest intensity was regressed as a function of time (R 2 = 0.789). Prior to 2030, we used the calculated function to estimate annual timber harvest, while it was assumed that the harvest amount would increase by 1% per year in the period 2030-2050. Thus, we assumed a progressive increase in harvesting, which is expected to exceed 8 million m 3 in 2050. When determining the PLAN scenario, we followed the recommendations from forest management plans where the annual allowable cut is defined at a rate of 75-90% of the annual stand volume increment [9]. The annual harvested amount in the PLAN scenario was determined by fitting the negative exponential function, with the dependent variable being the increase in the rate of planned harvest compared to the planned harvest in the previous year (R 2 = 0.358). To calculate the adjusted values, we used the planned harvest in 2004 of 4,162,662 m 3 [35] as a starting point. Until 2035, the degree of increase of the harvested biomass took place pursuant to the calculated function; after 2035, a constant level of increase of 0.04% per year was used. The HAZ scenario was also based on the PLAN scenario and assumed a 3% increase in harvesting levels on an annual basis, as well as four extraordinary natural disturbances, appearing in an interval of approximately 10 years and assuming an increase in harvest intensity over the following three years, which is usually a consequence of bark beetle outbreaks [12]. The HH and LH scenarios are primarily based on the assumption of the annual allowable cut of the PLAN scenario, but the intensity of harvesting in the HH scenario was increased by approximately 30% and in LH reduced by 40%. With the selected increase and reduction, we wanted to obtain a higher range of annual harvesting levels, which could be potentially beneficial for the evaluation of harvesting levels on carbon stocks.

Distribution of Felling by Forest Types and Age Classes
In a disturbance matrix, we defined the type and amount of felling by classifiers and age classes for each year of the simulation. The annual felling rate was first allocated to classifiers in line with the existing volume shares in 2014. The highest harvesting amounts were, therefore, directed towards classifiers with the largest volume share in the inventory database, such as FT4 and FT5 mixed forests (Table 1). For each classifier, the harvest amount was then distributed by age classes according to the proportions representing the actual timber harvest in the last 10 years (Table 2). These shares were estimated on the basis of information from the SFS harvesting databases. Two disturbance types, namely commercial thinning and final felling were defined, whereby for each year a portion of final felling was allocated to the oldest forests. The stand maps were used to determine the percentages of final felling, which were 16% (BAU, PLAN and LH), 18% (HAZ) and 20% (HH). Final felling affects the regeneration of stands, since the model starts a new succession stage with age equal to 0.

Evaluation of the Model Simulations
The development of carbon stock dynamics in Slovenian forests was simulated using the CBM-CFS3 for all five scenarios from 2014 to 2050. Carbon stocks were compared with the official estimates reported by Slovenia in the National Inventory Report (NIR) [36] in the framework of the UNFCCC for the period 2014-2018. The mean absolute percentage error (MAPE) [37] between the simulated and reported values of above-and below-ground biomass was calculated. In addition, we estimated and compared the net annual carbon stock change for the carbon pools. Conversion of carbon (C), expressed in tonnes, to carbon dioxide CO 2 (tonnes) was performed by use of Equation (2), where the constant 44/12 is used to convert carbon density into units of CO 2 . Negative values represent removals from the atmosphere, while positive values represent emissions to the atmosphere.

Results
The accumulation of carbon stock in living above-and below-ground biomass coincides with the predicted harvest intensity for each scenario, whereby the accumulation is the highest for the LH and BAU scenarios. Both of these scenarios envisage the lowest harvesting up to 2050 compared to the other scenarios ( Figure 4). An increase in carbon stocks in above-ground biomass is generally expected, as the anticipated harvesting levels are mainly lower than the annual volume increment in Slovenian forests of approximately 9 million m 3 [9]. After 2025, the HH scenario assumes that harvesting will exceed the total volume increment ( Figure 5). By 2050, the carbon stock in above-ground biomass increases by 28.4% (LH), 19% (BAU), 10% (PLAN), 6.5% (HAZ) and 1.2% (HH) compared to the base year of 2014. In the case of below-ground biomass, these shares were the highest in the LH (16.9%), BAU (8.4%) and PLAN (1.1%) scenarios, while stocks in below-ground biomass in the HAZ and HH scenarios decreased by 1.2% and 6.3%, respectively. The analysis of the temporal dynamics, taking into account the change in carbon stocks in living biomass (i.e., above-and below-ground), show that Slovenian forests will remain a carbon sink under all scenarios, with the exception of the HH and HAZ scenarios in some years (     Higher harvesting levels have a negative effect on carbon stock in litter and soil, but also a positive effect on the accumulation of carbon stock in deadwood (Figure 4). The effect on the carbon stock in the soil and litter becomes more pronounced after 2030 and 2035, respectively, when the differences between the scenarios begin to increase. Carbon stock in deadwood has doubled (scenario LH) or tripled (scenarios HH, PLAN and HAZ) in 2050 compared to the base year 2014. However, deadwood accumulation largely depends on conservation measures and the attitude of forest owners towards deadwood, and thus the results might not reflect the true dynamics in Slovenian forests.
The impact of harvesting on carbon stocks in stands of different age classes ( Figure 6) indicates the ageing of forests for all scenarios, which depends on the allocation of harvested biomass to specific age classes and a relatively low share of final fellings. Taking into account the mean values of the age classes, the average age of forests in 2014 is 87 years, and in 2050 it is between 110 (HH) and 119 (LH) years.   The results of the CBM-CFS3 simulations were compared with the official estimates reported by the Slovenian NIR in 2020 for the period 2014-2018 [36] in the framework of the UNFCCC ( Table 3). The results of model simulations show satisfactory match with the officially reported estimates of carbon stocks, where the MAPE is greater for below-ground biomass. MAPE was on average 4.5% for above-ground biomass and 9.16% for below-ground biomass. Better results were obtained for the BAU and LH scenarios, while for the HAZ, HH and PLAN scenarios, which assumed higher harvest intensity, the MAPE were greater. More specifically, MAPE for the comparative period 2014-2018 were lowest in the LH scenario (2.8% for above-and 7.8% for below-ground biomass) and highest for the HH scenario (5.7% for above-and 10.1% for below-ground biomass). Table 3. Comparison between simulated and reported (United Nations Framework Convention on Climate Change (UNFCCC)) estimates of carbon stocks in living biomass in megagrams for Slovenian forests in the period 2014-2018.

Pool
Year

The Impact of Harvesting on Carbon Stock Dynamics in Slovenian Forests
The CBM-CFS3 has been recognized as one of the most reliable tools for simulating carbon stock dynamics in forest ecosystems [39] and has often been used in similar simulations in Europe [29,40] and globally [27,41]. Initially, this model was developed for even-aged forests which are commonly found in far northern geographic regions. However, the model can also be applied to mixed and uneven-aged forests if the yield curves objectively describe the development of forest biomass over time [29].
The carbon stock dynamics in above-and below-ground biomass simulated by the CBM-CFS3 show a satisfactory match with the values reported by Slovenia in the framework of the UNFCCC (Table 3), whereby the results for above-ground biomass yielded lower MAPE. Greater MAPE in terms of below-ground biomass are to be expected since the entire calibration of the model is based on the inventory database which only encompasses above-ground biomass. National reports in the framework of the UNFCCC are based on the National Forest Inventory (NFI) [42][43][44], which is a different data source than that used in our study. Our simulations are based on data of the SFS which is collected in the framework of forest inventories for the purpose of forest management planning. The disadvantage of this data is that only 1/10 of forest management compartments are updated every year. Therefore, the data used are not completely up to date. Furthermore, the NFI database is based on individual tree species, which are systematically measured every 5-6 years and as such represent a more consistent data source. Further reasons for general underestimation of reported values could be related to the parameterization of the CBM-CFS3 model. Also, the selected comparison period (2014-2018) might be unsuitable for comparison due to many natural disturbances during this period. Discrepancies between the observed and simulated carbon stocks from the CBM-CFS3 model were also reported for red spruce (Picea rubens) forest in Eastern Canada [45].
The results of our simulations realistically reflect the dependence between the assumed future harvest and projected net carbon sinks. All simulated scenarios, with the exception of HH and occasionally HAZ, indicate the possibility of significantly higher annual harvesting than that seen in the previous two decades, which would still guarantee a carbon sink on a national scale. Higher harvest intensity would also make sense in terms of ensuring the sustainability of all forest functions. However, in the first accounting period from 2021 to 2025, the Forest Reference Level (FRL) will need to be considered when planning the amount of harvesting in all EU member states [46]. The proposal of a draft delegated act assumes an FRL value of −3270.2 Gg CO 2 eq per year for Slovenia in the period 2021-2025 [38], which according to our simulations, means limiting the harvest to about 6 million m 3 per year. The exact threshold at which the FRL value will be exceeded is hard to specify, since net sinks largely depend on the tree species being harvested [47]. In addition, the overall carbon balance also depends on wood use, and hence the inclusion of harvested wood products (HWP) and substitution effects in simulations of carbon pool developments can have a substantial impact on the overall carbon balance [48,49]. Finally, to meet the commitment under the LULUCF regulation [7], Slovenia should enhance the carbon sink in forests, promote the use of wood as a material and energy source, and ensure the stability of the forest land area.
The intensity and structure of harvesting are two important components when changing the age structure of forests. In the base year of 2014, the largest share of forests is represented in the age class of 81-100 years, while the largest share of forests in the final year of 2050 will be in the class of 101-120 years. The latter is due to the distribution of the harvest into age classes, where the largest share is directed towards semi-mature and old forests (approximately 90% of the harvested forests are in the age classes 80-140 years). Due to the short simulation period (37 years), all classifiers were advanced by only one age class. Nevertheless, we can conclude from our simulations that forest rejuvenation is a long-term process in which most of the harvest should be directed to the oldest forests and the regeneration of these forests, which would ensure a balanced age structure in the coming decades. Regardless of the accounting rules after 2030, the harvest in Slovenia will have to be increased towards the middle of this century to ensure that forests act as a carbon sink and remain resilient over the long term.

Sources of Uncertainty
The CBM-CFS3 has proven to be a reliable tool for simulating carbon stock dynamics at the national level. Yield curves represent key input data for estimating carbon stock dynamics with the CBM-CFS3 model. In the even-aged forests of northern latitudes, volume yield curves are often well-known [50]. In our case, these were calculated on the basis of the PSP data on growing stock and their estimated age from the calculated transition periods. The disadvantage of this procedure is that the yield curves thus calculated already contain a part of the harvest, which is reflected in lower growing stock. We estimated that the curves include an average of 2.5 million m 3 of harvest per year, which coincides with the low harvesting in the period 2000-2009 [9]. The adequacy of the correction choice was confirmed by simulations in the preliminary phase, in which the harvesting and volume increment were equalised and the carbon sinks were approximately 0. For uneven-aged heterogeneous managed forests, which are typical for Slovenia, yield curves are difficult to calculate precisely. We selected a relatively simple regression equation, while more complex approaches would include additional explanatory variables, e.g., site index, elevation and slope [43]. However, based on the satisfactory match between simulated and reported estimates of carbon stocks for the overlapping period 2014-2018, we assume the developed yield curves realistically reflect the development of stand biomass over time.
Sources of uncertainty are also related to the choice of ecological parameters for tree species mixture categories. In our study, we defined three different categories of species mixture (see Section 2.2.1), which must be connected to existing tree species in the archive index database. The broadleaf forests (BRD) were connected with the category OB_SI (other broadleaf trees of Slovenia) and coniferous forests (CON) with the category OC_SI (other coniferous trees of Slovenia). The mixed forest category (MIX) is not available for Slovenia, so we connected mixed forests with the category OB_SI. We also checked other options (e.g., OC_SI), but matching with reported estimates resulted in larger MAPE. This choice seems reasonable, as broadleaves are more common in Slovenia's mixed forests than conifers.
The results of our simulations provide a realistic picture of the impact of harvesting on net carbon sinks projections up to 2050 for Slovenian forests. However, some other studies reported a distinct turnaround of carbon sequestration after 2050 [51,52], which could also happen in our simulations, if a longer period was considered. In addition to harvesting, changes in carbon stocks in living biomass may also be affected by other processes, such as the degree of tree mortality and changes in forest areas, which were not covered in our study and should be addressed in future studies.

Conclusions
In this study, we used the well-established and tested CBM-CFS3 model to quantify the impact of various harvesting scenarios on the forest carbon dynamics of Slovenian forests. The results of the simulations showed a satisfactory match with officially reported estimates under the UNFCCC and highlighted the importance of properly selected harvesting intensities to ensure carbon sinks in future decades. Harvesting below 9 million m 3 per year will ensure that Slovenian forest land is a carbon sink, while harvesting should be around 6 million m 3 per year to comply with the FRL.