Shifting from even-aged management to less intensive forestry in varying proportions of forest land in Finland: impacts on carbon storage, harvest removals, and harvesting costs

Many studies have reported increased multi-functionality and financial profits due to a shift from even- to uneven-aged forest management. However, little is known (from long-term experiences or predictions) how alternative management systems could affect national-scale wood production and carbon storage, if adopted over very large areas. We analysed these effects using an area-based framework, in which multiple Markov chain models were used to simulate the development of forests according to different management systems. Classification of forests to wood availability categories was used to determine the system to be applied. We enhanced the framework to allow shifts between management systems that correspond to enforced or voluntary changes in forest use. Simulations of extensive shifts from conventional even-aged management to alternative silvicultural systems revealed interesting developmental patterns that cannot be directly deduced from studies that upscale from smaller areas. Our results show that the amount of carbon stored by Finnish forests can be increased by applying less intensive management systems, although this has trade-offs in terms of harvests and associated financial costs. The level of trade-offs differed depending on the type of forest that shifted between management systems and whether areas were also assumed to be completely set aside from forestry. These differences were further pronounced if the desired harvest levels and their allocation changed along with the management system. If the studied attributes were considered at the same relative scale and with equal weighting, the extensive shifts to alternative management systems exhibited the strongest impact on harvesting costs.


Introduction
Currently, several international and national strategies strive for a transition from a fossil-to a bio-based economy, which calls for an increase in the use of (forest) biomass for products, such as biofuels and energy, chemicals, polymers, and wood-based structures. However, concerns in regard to fragmentation, degradation, and loss of forest habitats have been increasingly expressed and international and national agreements have been signed to reverse this trend. In particular, the restoration of at least 15% of degraded ecosystems by 2020 (CBD 2010), also known as Aichi Target 15, is widely accepted as a national and global conservation target. However, maintaining both high economic forest yields and the viability of forest species involves trade-offs. In Finland, for example, Kotiaho et al. (2016) have estimated that meeting Aichi Target 15 would cost between 12 and 23 billion euro, or 368-658 million euro per annum, when only forests and peatlands are considered and if the conservation actions were prolonged until 2050. For comparative purposes, annual stumpage earnings amount to approximately 1.5 billion euro (MAF 2015).
Communicated by Miren del Rio.

3
Rather than segregating the forest and conservation sectors, a good compromise could be achieved by combining cost-effective conservation actions with the sustainable use of forests and peatlands (see also Kotiaho et al. 2016). For example, the management of a proportion of forest land as "multi-use conservation landscapes" (MUCLs; Hanski 2011) is a practical and cost-effective means to aid forest conservation strategies, in addition to the strict protection of forest networks (see also Stevens and Montgomery 2002). More precisely, Hanski (2011) has proposed that management of a third of the land as MUCLs, and the strict protection of a third of this area, would be sufficient for the conservation needs of specialist species, particularly if the MUCLs were composed of aggregated habitat patch clusters and could be connected to existing protected area networks (see also Rybicki and Hanski 2013). However, the work cited above does not specify feasible silvicultural practices or management intensities for the proportion of MUCLs (twothirds) that are not strictly protected. Yet, we can assume that such areas could be managed according to the principles of multiple-use (Fürstenau et al. 2007) or integrated forest management (Diaci et al. 2011), which support the production of ecosystem services other than just species conservation based on a more diverse set of silvicultural practices, compared to conventional rotation forestry or even-aged management systems.
Forestry practices in Finland have been based on evenaged management since World War II (Kuuluvainen et al. 2012), but options for forest management practices will clearly increase in the future. By options, we refer to the various forms of uneven-aged management, such as continuous cover forestry as defined by Pukkala (2016a). Continuous cover forestry essentially differs from even-aged management in that it avoids clear felling and planting by utilizing thinnings from above and by promoting natural regeneration. These choices may convert stands towards uneven-aged forest structures, although converging to a steady-state structure of any kind is not required (Pukkala 2016a). Continuous cover forestry is expected to become more common, because of its potential to supply multiple ecosystem services (Pukkala 2016b; Peura et al. 2018) and reduce the financial costs related to regeneration and other silvicultural operations (Pukkala 2016a) compared to even-aged management (see also Knoke 2012;Kuuluvainen et al. 2012;Puettmann et al. 2015;Nieminen et al. 2018).
Comparisons of alternative forest management systems and subsequent trade-off analyses are typically based on long-term observations (Sutherland et al. 2016;Strengbom et al. 2018), metamodelling (Lafond et al. 2017) or simulations. The latter have been carried out at the forest stand or small forest holding level (Pukkala et al. 2011;Pukkala 2016a, b;Carpentier et al. 2017), the landscape level (> 100 km 2 ; Triviño et al. 2015;Diaz-Balteiro et al. 2017;Peura et al. 2018), and the regional level (> 1000 km 2 ; Schröter et al. 2014;Pang et al. 2017). However, long-term experiences or predictions as to how alternative practices may affect national-scale wood production if adopted over very large areas are not known. With the exception of recently formulated growth (Bollandsås et al. 2008;Pukkala et al. 2013) and thinning models (Pukkala et al. 2015;Vauhkonen and Pukkala 2016), most conventional forest simulators and projection tools have been developed for even-aged forestry systems. Conventional models for forest development would, therefore, extrapolate outside the original population if applied in uneven-aged forests. Moreover, detailed forest-specific thinning instructions might not reconcile with large-area projections based on aggregated spatial scales (cf., Verkerk et al. 2014;Creutzburg et al. 2017;Mouchet et al. 2017). As such, there is a need for flexible tools that can combine detailed instructions with projection capabilities for large areas.
From the perspective of regional or national-level wood production, the areas subject to conservation or integrated management reduce the amount of Forests Available for Wood Supply (FAWS; Alberdi et al. 2016). From the point of view of provisioning of non-wood forest products or other ecosystem services, it is useful to also simulate the development of area and growing stock in the remaining areas-i.e. in Forests Not Available for Wood Supply (FNAWS) and Forests with Restrictions on Availability for Wood Supply (FRAWS; see also Vauhkonen and Packalen 2017). Extensive forest inventories, such as the National Forest Inventory (NFI), also provide data for the simulation of transitions due to growth or management for forests prioritized for uses other than solely wood production. Several approaches have been presented for the projection of the future development of forest resources based on transition probability matrices of stand-specific diameter classes (e.g. Bollandsås et al. 2008;Schou and Meilby 2013;Roessiger et al. 2016). However, corresponding Markov chain models based on transition matrices of forest size and structure classes (e.g. Vauhkonen and Packalen 2017) could be more suitable for area-based projections of forest dynamics based on the NFI data. Vauhkonen and Packalen (2017) simulated the development of forest size and structure classes derived from NFI data by combining multiple Markov chain models for different silvicultural systems. In their study, classification of forests to wood availability categories determined which system was applied. Wood availability depended on administrative forest use restrictions as recorded in the NFI data. Areas where forestry operations were forbidden (FNAWS) or restricted to selective harvests for the enhancement of ecosystem services other than wood supply (FRAWS) were distinguished from FAWS. The future development of FAWS was simulated with even-aged management; FRAWS with 1 3 continuous cover forestry, in which the final felling was replaced by thinning from above; and only the natural processes were simulated for FNAWS. The assumption that restrictions on forest use determine the silvicultural system is a simplification that might be unrealistic, especially, if large proportions of FAWS diverge from even-aged management in the future, as reasoned above. However, the FAWS in our inventory data are largely composed of even-aged stands; forest management instructions were based on even-aged silviculture until the final measurements of the 11th Finnish NFI in 2013. Thus, increasing the proportion of alternative management systems in FAWS fundamentally means that more forests will be managed in a similar fashion to FRAWS or FNAWS are obligatory managed because of administrative forest use restrictions. We propose that simulation of shifts between these categories should allow further examination of the national-scale effects of moving from evenaged, rotation forest management to alternative silvicultural regimes, which may result from enforced political decisions to promote less intensive forestry or increased nature conservation or voluntary changes in the use of forests.
The objectives of this study are (1) to enhance the Markov chain modelling framework by the simulation of shifts between forest management systems and (2) to use the developed framework to assess the trade-offs between carbon storage, harvest removal, and harvesting costs due to these shifts. We first derived benchmark projections of the future development of forests according to silvicultural systems associated with current wood availability categories as described in the previous paragraph. We then repeated the simulations with varying proportions of even-aged forests managed according to continuous cover forestry or set aside to assess the effects of diverging from business-asusual management in an increasing proportion of forest land.

Simulation framework and data
We simulated the future development of forests in Finland using an area-based Markov chain model (cf. Vauhkonen and Packalen 2017). Our analyses considered almost the entire forest land of Finland, excluding northernmost areas and forests located in the southern archipelago, which are considered to have a low importance in terms of wood supply of Finland. The total area of forests on productive and poorly productive forest land was 21.28 million ha, with an initial growing stock of 2234 million m 3 (the 11th Finnish NFI-NFI11;measurements in 2009, which equates to approximately 95% of the entire growing stock in Finland. As in Vauhkonen and Packalen (2017), the forests were assigned to wood availability categories according to forest use restrictions in the NFI data: initially, 10.1% and 10.6% of the total forest area were classified as FNAWS and FRAWS, respectively, with the remainder classified as FAWS. In order to study the national-scale effects of diverging from the current management system, varying proportions of forest were removed from FAWS and assigned to either FRAWS or FNAWS. These were subsequently simulated for future development along with forests that are currently assigned to these categories due to administrative forest use restrictions. Specifically, forests that received even-aged management in the benchmark simulations were either assigned to continuous cover forestry (in the proportion transited to the FRAWS category) or set aside (FNAWS) at the beginning of the simulations.
The selection of land transited from FAWS to the other categories was simulated according to four different strategies that mimic the different drivers and land availability in the transitions. The aim was to mimic the establishment of either (a) MUCLs as a variant of the concept proposed by Hanski (2011), which consist of both continuous cover forestry and set-aside areas; or (b) multi-use landscapes (MULs) managed by continuous cover forestry without setaside areas. The four strategies were obtained by emphasizing the selection of forests with high (strategies MUCLhigh, MULhigh) or low (MUCLlow, MULlow) proxy conservation values, determined as a function of the maturity of the trees, tree species composition, and site fertility (Appendix 2; Lehtomäki et al. 2015). Because of its formulation, the proxy might act as a surrogate not only for the biodiversity features of conservation interest, but also the more general multiple-use potential of a forest. As further discussed in Sect. 4.1, the four strategies generally cover the range of production possibilities, by reducing the management intensity in extensive proportions of the different forest types.
Each land transition strategy was composed of eight proportions p, where p was either 5, 10, …, or 40% of the FAWS land area, which was then added to either FRAWS or FNAWS. The proportions were selected by first computing the conservation value (consval) of each NFI plot, as described in Appendix 2, and the cumulative distribution of consval separately for each of the 15 forestry regions (or Forest Centres as distinguished by the Finnish NFI). The number of plots required to represent p% of the FAWS land area was selected at equal intervals from the cumulative distribution of each forestry region as follows: • MUCLhigh p% of FAWS was first selected according to consval. A third of this area was further selected according to the cumulative distribution of consval computed for plots in the selected p% and assumed to transit to FNAWS (no management). The remaining proportion of p% was assumed to transit to FRAWS and to be managed as continuous cover forestry. The emphasis in the new 1 3 MUCLs established according to this strategy was on FAWS with a high conservation value. One-third of this selection (with an emphasis on forests with the highest conservation values) was set aside. • MUCLlow As with MUCLhigh, but the land area to transit was selected according to the inverse of consval, excluding plots with consval = 0 (i.e. bare land). The emphasis in the new MUCLs was on FAWS with a low conservation value. One-third of the forests with the highest conservation values in the first selection were completely set aside. • MULhigh p% of FAWS was selected according to consval and assumed to transit to FRAWS to be managed as continuous cover forestry. The emphasis in the new MULs was on FAWS with high conservation value. No additional forest areas were set aside. • MULlow As with MULhigh, but the land area to transit was selected according to the inverse of consval, excluding plots with consval = 0 (i.e. bare land). The emphasis in the new MULs was on FAWS with low conservation value. No additional forest areas were set aside.

Simulations of forest development
As the Markov chain model, we used v. 2.0. of the European Forestry Dynamics Model (EFDM), which is implemented in the R statistical modelling environment (R Core Team 2016) and can be downloaded from https ://githu b.com/ecjrc/efdm as open source under the European Union Public License (EUPL). The EFDM approach is based on arranging the forest area into matrix cells according to ecological and socioeconomic factors and simulating the development of the resulting matrices. In the simulations, the area represented by each matrix cell is managed according to a set of pre-defined activities, which may transit the area to other cell(s) depending on the transition probabilities associated with the activities. In practice, the forest area distribution after one simulation step (i + 1) is obtained as a multiplication of the area in state i by the probability that the area receives one of j activities and the activity-conditional transition probabilities (for details, see Sirkiä 2012or Packalen et al. 2014).
The simulations were carried out in five-year time-steps, which correspond to the measurement interval in the data used to derive the transition probabilities. In total, ten steps were simulated (i.e. the last year in the simulation period is approximately 2060, depending on the initial measurement year). The initial forest state, (business-as-usual) transition and activity probabilities, and output coefficients for the model were derived from the NFI11 data according to the workflow presented by Vauhkonen and Packalen (2017) based on using permanent NFI plots as pairwise data for natural processes. Appendix 2 also illustrates the process of deriving the initial state and the simulation of the future development of one example forest. In the following sections, we briefly describe the parameterization, but note that parameters and their effects are explained in more detail in Vauhkonen and Packalen (2017).

Initial state
The EFDM is parameterized by matrices with dynamic (e.g. age, volume) and static (e.g. geographical region, site fertility) dimensions (henceforth "factors"). Similar to Vauhkonen and Packalen (2017), we defined the dynamic factors separately for forests to be managed using different silvicultural attributes, such as age and volume for evenaged management systems or stem number and volume for continuous cover forestry. The natural processes of FNAWS were simulated using age and volume matrices; no differences would have been observed if stem number and volume matrices were applied (Vauhkonen and Packalen 2017). The continuous measurements from the NFI were classified using age classes of 0, 5, 10, …, 120, 120+ years, whereas the class limits for both the volume and stem number were determined as the values of the 10th, 20th, …, 90th and 95th quantiles of the pairwise observations made from the permanent NFI plots. The class limits (presented in Appendix 1) were defined by Vauhkonen and Packalen (2017) with a motivation to obtain an approximately equal amount of pairwise observations per class and to reduce the value range of the last class by halving the number of observations included in it. These dynamic factor matrices were derived separately by applying the following static factors: (1) known land-use restrictions: FAWS, FRAWS, FNAWS; (2) forest ownership: private, public + other; (3) site fertility: a total of five categories that correspond to the four taxation classes that productive forests are assigned in the Finnish NFI + a fifth class that includes all poorly productive forest land; (4) dominant species: pine, spruce, deciduous trees.

Management activities and their probabilities
Possible management activities were "no management" (i.e. simulation of natural processes only), "thinning", and "regeneration harvest". A "thinning" always referred to a management thinning and was implemented as a thinning from below for both FAWS and FRAWS. A "regeneration harvest" was implemented either as a final felling in the even-aged management system (simulated for FAWS) or a thinning from above in continuous cover forestry (for FRAWS), as described in Sect. 2.2.3.
The probabilities of the activities were determined in two steps. First, the NFI data were used to compute the mutual proportions of the activities, resulting in two alternative allocations of the activities: (1) a business-as-usual allocation 1 3 (A BAU ) based on the proportions of permanent NFI plots not managed, thinned, or regenerated during the most recent 5-year period; or (2) a schoolbook-allocation (A SB ) based on the proportions of the NFI plots that should be managed within the next 5 years strictly according to forest management instructions. The proportions were based on the NFI records that were made according to the instructions that prevailed at the time of the measurements [for more details, see Yrjölä (2002)]. The two alternative allocations are class specific (cf. Figure 1 of Vauhkonen and Packalen 2018); following the initial proportions of A SB would, in general, involve proposing much more management activities than where all forests are considered as FAWS and simulated according to the even-aged management system 1 3 A BAU . Second, the activity probabilities were iterated to the level that yielded a specified total roundwood harvest level. As in Vauhkonen and Packalen (2017), the iteration was carried out by repeatedly multiplying the activity probabilities ≠ 0 by 1.01 or 0.99, depending on the sign of the difference between the goal and the harvest level given by the current activity probabilities, until the harvesting goal was met or the activity probabilities could not be changed. The specified total harvest goal was calculated from two harvest levels provided in the National Forest Strategy of Finland (MAF 2015): the business-as-usual level (65 million m 3 /a) and the desired future level (80 million m 3 /a). Because our analyses considered about 95% of the total growing stock in Finland, both goals were multiplied by this proportion to yield the final harvest targets of approximately 62 million m 3 /a or 76 million m 3 /a, respectively. Activity probabilities were iterated to meet these harvesting goals only at the beginning of the simulations. Thereafter, the same proportion of the land area was managed in every simulation step, i.e. the volume harvested after the first simulation step depended on how the forest class distribution evolved during the simulations.

Transition probabilities
The transition probabilities of the natural processes (growth) were derived using pairwise observations from the permanent NFI plots (total: 11,987 or about 23% of the plots), which were measured at approximately five-year intervals between NFI11 and the earlier inventory (NFI10). Positive differences in total volumes on plots with no treatments, based on data that could be matched with certainty between the two subsequent inventories, were recorded as the pairwise data. The estimated transitions, therefore, included growth and mortality, but not potential reductions due to calamities or natural disturbances, for example.
The transitions due to management activities were based on simulations of their expected development. The forests affected by final fellings were forced to transit to the beginning of the even-aged rotation. Their early development was simulated conservatively, such that 25% of the final-felled area moved to volume class #2 (Appendix 1) in the first simulation step after final felling, and the remainder thereafter. A thinning simulator was implemented to derive pairwise observations due to the treatments. The simulator determined the trees to be removed following two types of instructions: (1) a thinning from below corresponding to conventional instructions for forest management and (2) a thinning from above with an intensity corresponding to an interest rate of 3%, as predicted by Eq. 2 in Pukkala et al. (2015). As detailed in Vauhkonen and Packalen (2017), these simulations were applied to plots with an initial basal area > 10 m 2 /ha and a mean height > 10 m, which corresponds to commercial thinnings. In addition, pre-commercial thinnings were simulated for the less mature plots with a thinning need recorded by the NFI. These thinnings were implemented as thinning from below, but instead of applying a thinning curve for the harvest removal, the aim was always to retain a residual stand with approximately 1000 trees/ha. In mimicking an operational implementation, the trees to be cut were distributed to different parts of the diameter distribution, as described in detail by Vauhkonen and Packalen (2017). The thinnings took place at the beginning of each simulation step. The growth of the forests thinned from below was simulated by applying the transition probabilities of forests not managed in the simulation step where the thinning took place.
The transition probabilities were simulated using the same pairwise data for both age-volume and stem numbervolume classes and can be expected to develop similar to established forests (cf. Vauhkonen and Packalen 2017). However, differences may occur between the management systems due to assumptions on the early development after the regeneration harvest. To assess the sensitivity of these assumptions on the management systems, all figures of Sect. 3 were alternatively reproduced using the same methodology, but with different assumptions for the early development of even-aged and continuous cover forestry. Specifically, the parameters above were modified such that 75% of the final-felled area moved to the next volume class that was already in the first simulation step after final felling; and the forests thinned from above were simulated for growth in the same simulation step that the harvest occurred, similar to forests thinned from below, i.e. they were assumed to recover rapidly despite heavy thinning. A comparison of the full results indicated increasing differences, i.e. uncertainties towards the end of the simulation period. However, the management systems were not essentially different from each other, based on either of the results, so our analyses focused on the results from the conservative early development simulations, while the results of the more rapid development are provided as Electronic Supplementary Material for the reader who wishes to evaluate the degree of sensitivity in the simulations.

Output coefficients
As the simulations only considered the development of the forest area distribution, separate transformation coefficients were determined to derive further information for the variables of interest (carbon storage, harvest removal, and harvesting costs). The coefficients were determined as the mean values of the NFI plots within the dynamic classes (Appendix 1) as follows: 1 3 • Carbon storage The total above-and below-ground biomass of each NFI plot was first computed by estimating stem, branch, foliage, stump, and root biomass with the models described in Repola (2008Repola ( , 2009). The biomass components were multiplied by species-specific expansion factors (approximately 0.5; see Table 1 in Pukkala 2014) and summed to obtain the carbon content in each plot. The transformation coefficients were computed as mean values for the volume classes (Appendix 1) using all NFI plots. • Harvest removal The proportions of log-and pulp wood of the total volume obtained from the different harvests were first computed from the estimates of assortment volumes obtained by theoretically bucking the trees in the NFI plots, following the cut-to-length harvesting method. The total fellings were computed by summing these proportions. Coefficients for the thinnings were determined as mean volumes of logs for sawn wood and pulp wood from the NFI plots for which the thinnings were simulated. The coefficients for the final fellings were based on similar mean values in all NFI plots and determined using both age and volume classes. • Harvesting costs The time expenditure of cut-to-length logging and roadside-transportation of the trees in the NFI plots was estimated using the models described in Rummukainen et al. (1995). The models were applied with the assumption that one hectare of forest represented by an NFI plot was cut to the timber assortments as described above. The models assumed higher time expenditure for thinning-types of harvests compared to final fellings, but typical terrain conditions and with no entry time for the logging equipment, i.e. the time expenditure values only depended on the amount of timber assortments obtained as a result of the different types of harvests. As in Vauhkonen and Pukkala (2016), the harvesting costs were computed with the assumption that operating a harvester and a forwarder cost 80 €/h and 57 €/h, respectively. The transformation coefficients for harvesting costs (expressed as €/m 3 ) were computed as mean values for the volume classes based either on all NFI plots (final fellings) or plots for which the thinnings were simulated (thinnings).

Total effects under business-as-usual (BAU) management
Assuming BAU management and the current wood availability categories, the total carbon storage in the above-and below-ground living biomass was 883.2 million tonnes carbon at the end of the simulation period. In total, 4057.1 million m 3 of roundwood was harvested during the simulation (average: 73.8 million m 3 /a, with a variation from 61.9 to 79.5 million m 3 /a between the simulation steps). The unit costs of harvesting varied from 16.9 to 21.5 €/m 3 (average: 18.85 €/m 3 ) between the simulation steps. Figures 1 and 2 show how these values developed according to the varying proportions of wood availability. All land transition alternatives resulted in increased carbon storage in the above-and below-ground living biomass at the end of the simulation period, compared to current wood availability (Fig. 1, top panel). The MUCLlow alternative resulted in the highest levels of carbon storage (903.8-1090.1 million tonnes carbon, depending on p%). The higher the p% value assigned to MUCLlow, the greater the difference to the other land transition alternatives. The latter behaved similarly when compared to each other in terms of the increase in carbon storage and the magnitude of this increment as a function of p% (Fig. 1, top panel).
The reduction in potential harvest removal was greatest in the MUCLhigh alternative, followed by the MUCLlow alternative ( Fig. 1, middle panel). Total harvest removals (computed for 55 years by temporally allocating the harvests to the beginning of the simulation steps in EFDM) varied from 3984.5 to 3371.2 million m 3 in MUCLhigh, and from 3980.0 to 3508.6 million m 3 in MUCLlow. Removal depended on the p%: the higher the value, the greater the difference to the harvest removals of BAU management and wood availability. Harvest removals in the other two landuse transition alternatives varied from 4051 to 3887 million m 3 according to p%, i.e. the difference to BAU management was much less compared to the MUCLhigh and MUCLlow alternatives. The difference in harvest removals produced by the other two strategies remained constant despite the increase in p%. The MULhigh alternative resulted in the least reduction in harvest removals among the considered land transition alternatives.
Harvesting costs increased in conjunction with the increasing proportion of land transited from FAWS, compared to BAU management and wood availability. This result applied to all land transition alternatives. However, both the MUCLlow and MULlow alternatives increased costs slightly, with unit costs varying between 19.0 and 21.15 €/ m 3 , depending on the p% value. In addition, these costs were not affected to any extent by the increase in p%. On the contrary, there was a strong positive correlation between the p% value and costs in the MULhigh and MUCLhigh strategies. The latter strategy also resulted in the highest unit costs (19.9-30.2 €/m 3 , depending on the p% value).

Temporal effects under BAU management
Carbon storage in the above-and below-ground living biomass, harvest removals, and harvesting costs developed 1 3 over time until the end of the simulation period (Fig. 2). As mentioned in Sect. 2.2.2 and further discussed in Sect. 4, the harvests of the first simulation step were iterated to the same level in all land transition alternatives, so the differences begin to show immediately after the first simulation step. Carbon storage started to increase in approximately two steps (10 years) after the land transitions, whereas the increase in the harvesting costs took place immediately. With respect to carbon storage, the land transition alternatives differed in magnitude, but the trends were similar over time. The land-use transition alternatives differed more with respect to harvest removals and costs.
While Fig. 1 suggests that the increase in p% almost equally affected the total harvest levels based on the MULhigh and MULlow land transitions, these alternatives clearly differed in terms of temporal development pattern over the entire simulation (Fig. 2). In particular, the harvests of MULlow were closer to BAU over the first few simulation steps, whereas the harvests of MULhigh approached those of BAU towards the end of the simulation, eventually exceeding that level. MULhigh also differed from the other strategies in that the highest p% value allowed the greatest increases in harvest towards the end of the simulation. The temporal development of harvesting costs also differed depending on land transition strategy. In alternatives with emphases on high consval, the costs increased considerably due to the land transition, but their temporal development resembled BAU (i.e. the costs decreased slightly over time). In alternatives with emphases on low consval, the costs did not initially differ from BAU, but instead increased over time.

Joint effects due to shifts between wood availability categories and changes in total harvesting levels and allocation
Altering harvest levels and their allocations affected the levels of carbon storage in the above-and below-ground living biomass, harvest removals, and harvesting costs under current wood availability, as illustrated in Fig. 3 in a scale normalized to the level of the initial values of these attributes assuming BAU management and current wood availability (grey bars in the leftmost column of Fig. 3). Relative to that situation, more harvest removals were obtained at lower costs when the harvest allocation was changed from A BAU to A SB . This change also increased carbon storage towards the end of the simulation. Carbon storage decreased and harvesting costs increased when the A BAU harvest allocation was maintained, but the proportion of harvested areas was increased at the beginning of the simulation to correspond with the harvest goals outlined in the National Forest Strategy. However, the allocation of increased harvest levels (according to A SB ) somewhat compensated for these effects Fig. 3 The temporal development of carbon stock (top row), harvest removals (middle), and harvesting costs (bottom), assuming different degrees of wood availability and harvesting. The y-axes are presented in a scale normalized to the level of the initial values of these attributes assuming business-as-usual activity and transition probabilities, the development of which is illustrated in sub-figure "A BAU ; 62 mill" according to the applied allocation of harvests and the total amount of harvests in m 3 /a for the first simulation step. The grey bars depict the development assuming current wood availability, while the coloured lines indicate p = 25% of the land area of Forests Available for Wood Supply (FAWS) shifted to Forests with Restrictions on Availability for Wood Supply (FRAWS) or Forests Not Available for Wood Supply (FNAWS) according to the different strategies (refer to Fig. 1 caption for the interpretation of the symbols). The black lines indicate the hypothetical situation where all forests are considered as FAWS and simulated according to the even-aged management system 1 3 (the rightmost column of Fig. 3). Notably, the distribution of harvests differed over time (cf., Fig. 3), in that reduced harvesting activities at the beginning of the simulations resulted in faster ingrowth of forests in the older and more densely stocked classes. Therefore, levels of harvested roundwood tended to increase over time with lower harvesting targets and, conversely, decrease with higher targets.
A comparison of land transition alternatives with respect to the aforementioned effects and against each other with p = 25% [selected as this level yields the total area of restricted wood availability of 8.7 million ha, i.e. approximately double the initial level (4.4 million ha)] is shown in Fig. 3. In addition to the temporal development, the total effects of increasing wood availability restrictions can be presented as trade-off curves (Figs. 4, 5), using a similar relative scale as above. Figures 3, 4, and 5 indicate that the level of trade-offs and the relative effects due to adding restrictions depended on which attribute of interest, land shift strategy, and harvesting scenario was considered. In general, increased restrictions on wood availability produced the most severe effects on harvesting costs, followed by harvest removals and carbon storage in the above-and below-ground living biomass. These effects were especially pronounced with the MUCLhigh alternative ( Fig. 5; up to 30-84% increase in costs with p = 40%, depending on the harvesting scenario), but were considerably less with the other alternatives (up to 24-46%). Compared to MUL, the MUCL alternatives produced the greatest relative effects in all the studied attributes. However, more relevant here was whether the land that was shifted came either from forests Fig. 4 Trade-offs between carbon stock and harvest removals. Refer to Fig. 3 caption for the scaling and differences between sub-figures 1 3 with high or low consval. When the studied attributes with equal weightings (at the same relative scale as in Figs. 3, 4, and 5) were considered, the MULlow and MUCLlow strategies resulted in greater carbon storage and fewer effects on harvests and their associated costs than those with a stronger emphasis on high consval.

Using NFI data to simulate transitions between wood availability categories
Our work contributes to the development of the EFDM approach, in that we extended the area-based Markov chain modelling framework by considering how NFI data could be used for the simulation of shifts between forest management systems. Throughout the study, we use the term "shifts between wood availability categories" to refer to the principle of associating FAWS, FRAWS, and FNAWS (cf., Alberdi et al. 2016;Vauhkonen and Packalen 2017) with different silvicultural systems and shifted areas, represented by NFI plots within these categories. Such large-area changes could result from the adoption of less intensive silviculture practices due to voluntary or enforced changes in forest use. Similar principles could also be employed for simulating shifts between afforestation/reforestation and deforestation categories (e.g. Grassi et al. 2012), which are often included in the modelling of carbon dynamics due to joint forestry and other land-use changes. The areas shifted between the wood availability categories were selected based on conservation value predicted for each NFI plot as a function of site characteristics and species-specific volume and mean diameter of the growing stock . Similar proxy values have been used for conservation (Lehtomäki et al. 2009;Arponen et al. 2012;Sirkiä et al. 2012) or management prioritization (Vauhkonen and Ruotsalainen 2017), when plant mapping data have not been available for more detailed analyses. The application of the functions yielded the highest conservation values for forests with large trees, more than one species, and high site fertility. Aside from ecology and conservation, the maturity, average tree size, and species composition of a forest are more indicative of the general multiple-use potential of a forest, because similar attributes are also used to predict recreational and scenic values and yields of non-wood forest products (Pukkala et al. 1988;Miina et al. 2016).
Our MUL and MUCL alternatives were variants of the "third-of-third rule-of-thumb" designed for conservation planning (Hanski 2011). By assuming proportions other than the constant third and different management practices for the non-protected area, we can draw conclusions on the likely impacts due to the shift between alternative forest management systems beyond what is solely related to conservation. Following Hanski (2011), the land area that shifted between the wood availability categories was "located as evenly as possible across regions and countries to guarantee representativeness". However, a spatially explicit approach was not used in allocating MUCLs close to existing protected area networks, as recommended by Rybicki and Hanski (2013). It could also be possible to assess national-scale impacts due to such decisions, based on spatial conservation prioritization frameworks (Lehtomäki et al. 2009. Both the MUCL and MUL alternatives were simulated with the assumption that emphases on the shifted land be placed on forests with either high or low conservation proxy values. For example, forest conservation would obviously benefit if the forests with the highest conservation values were set aside, which might not be feasible for at least two reasons. First, as mentioned above, detailed plant mapping data are not available for the prioritization of forests for conservation (when examined at broad geographical scales), so analyses must instead be based on indirect conservation value proxies derived from general forest inventory data (cf., Lehtomäki et al. 2009Lehtomäki et al. , 2015. Second, based on such information, forests that are mature and highly stocked are considered of high conservation value, but are obviously valuable for alternative uses as well. Due to the opportunity costs and the restricted availability of forests with the highest conservation values (as measured by the consval index), any practical implementable extension to an existing conservation network must also include areas with lower values (see also Schröter et al. 2014;Lundström et al. 2016). Forests with lower values could be considered valuable for the production of specific ecosystem services but also for conservational purposes, if assessed after a few growing seasons (cf., Lundström et al. 2016).
From the discussion above and the results obtained in this study, it is clear that future forest projections will vary depending on assumptions of land availability for integrated management. Yet, opportunity costs or land availability are rarely considered, even in spatial prioritization studies (but see Moilanen et al. 2011;Schröter et al. 2014). Pang et al. (2017) assumed that the proportion of the landscape known to have high nature conservation, cultural, or recreational values would be subject to less intensive management. Peura et al. (2018) assumed that alternative management was adopted over an entire landscape, which is not realistic in practice due to fragmented land ownership and, consequently, non-uniform management objectives across the landscape. In their simulations based on the Finnish NFI10 data, Alrahahleh et al. (2016) increased the conservation area by 10% or 20%; when the sample plots used for conservation were selected in random (but with a probability related to basal area), sites not feasible or available for conservation may be emphasized, as reasoned above. In Solberg et al. (2017), the allocation was solved by means of optimization. It is unclear to what degree the aforementioned studies may exaggerate future development, if the shifts between management systems do not occur at the assumed magnitude. In contrast to many earlier large-area analyses (see also Verkerk et al. 2014;Creutzburg et al. 2017;Mouchet et al. 2017), our simulations were run with varying proportions of land shifts determined by different strategies. Using this approach, we largely circumvent the problem related to placing assumptions on the type of management applied in a specific area. The curves based on the simulations with the emphases on high and low consval could be interpreted as a range of production possibilities that are feasible depending on land availability.

Limitations and strengths of the Markov chain modelling approach in forest development simulations
Our projection of the forest development that follows the simulated shifts between wood availability categories was based on a similar parameterization of the Markov chain model, based on permanent NFI plot data, described by Vauhkonen and Packalen (2017). Here, we implicitly assume that even if the shifts between the wood availability categories were to take place, the future development of the forests in these categories corresponds to NFI observations 1 3 from current FAWS, FRAWS, and FNAWS. The realism of these assumptions should be considered, especially from the perspective of the increased use of thinnings from above due to the shifts to the FRAWS category.
As discussed in more detail by Vauhkonen and Packalen (2017), the data used for computing the transition probabilities should include an adequate number of observations of plots thinned during earlier management periods to correctly reproduce post-thinning recovery and development in the subsequent simulation steps. The validity of this assumption should be examined from two aspects. First, in this study we did not have the means to explicitly determine how many plots in the NFI data were at a developmental stage comparable to forests thinned from above. Instead, we relied on the extent of the NFI sample, but acknowledge that as most of the managed forests in Finland are traditionally thinned from below, the predicted post-thinning stand dynamics might not perfectly coincide with uneven-aged forests. Second, a fundamental principle related to Markov chains is the assumption that the future forest state depends only on the present state and not on the preceding events. Therefore, in principle, there is nothing to prevent the same areas from being thinned more often than would be feasible, with respect to a proper recovery period, in the subsequent simulation steps. However, similar considerations could also be applied to other studies that have compared even-aged management and continuous cover forestry. For example, while Peura et al. (2018) used growth models formulated for uneven-aged forests (Pukkala et al. 2013) to simulate post-thinning dynamics, it is unclear whether the forests in the model fitting data of Pukkala et al. (2013) were selectively harvested in sequences, e.g. approximately every 15 years that was assumed to be continued for 100 years (Peura et al. 2018). In other words, even if a growth model suggests that selective harvests can be continued in such sequences, it is unclear if this is feasible in practice. The Markov chain model only projects the development of the forest area distribution, whereas all other outputs must be derived indirectly via coefficients. In addition to the averaging of the discrete classes used in the matrices, as further elaborated by Vauhkonen and Packalen (2017), there may be inaccuracies in the models used to derive the output coefficients. First, the degree of variation explained by the models may vary, especially in the case of biomass components and conversion factors (cf., Neumann et al. 2016). Second, the coefficient values for the harvesting costs cannot be directly obtained as averages of NFI data, but are based on a modelling chain from time expenditure to costs. We acknowledge that the unit costs reported in Sect. 3 are higher than those that would result from realworld harvesting operations or reported in other studies (e.g. Kärkkäinen et al. 2018). Assuming business-as-usual management, our modelling chain yielded unit costs of 21.5 and 18.85 €/m 3 for the first simulation step and the entire simulation, respectively, whereas the realized costs varied from 10 to 12 €/m 3 in 2007-2016 (Strandström 2017). This difference may be related to the assumption that the harvests were focused on individual NFI plots, whereas the allocation of real-world harvests would obviously be based on aggregating forests that are to be harvested or otherwise more detailed economic reasoning. If considered as relative values between the management systems, however, the time expenditure should correctly indicate the impacts of acquiring an increased amount of wood from thinnings rather than from final fellings, and coefficient values expressed as costs are likely to be more informative than time expenditures. The time expenditure and costs will vary between thinning from below and thinning from above approaches (Appendix 2), but the models formulated by Rummukainen et al. (1995) might not be adequately representative of modern forest operations.
An overall strength of the EFDM approach with respect to all the studies discussed above is the flexibility of the implementation and, therefore, the potential to run multiple scenarios with slightly changed assumptions. Due to the use of NFI data and the area-based matrix model for simulating the development of the forest area distribution, our results are valid for large-area level projections (here for the entire Finland). Yet, the growth and harvesting simulations are based on the individual NFI plots, thereby facilitating a very detailed modelling of the transitions according to different forest types, for example. Therefore, the results of the simulations are less affected by ecosystem model assumptions (cf., Alrahahleh et al. 2016) or optimized allocation of silvicultural treatments (cf., Hynynen et al. 2015;Heinonen et al. 2017Heinonen et al. , 2018Solberg et al. 2017). As demonstrated above, expert opinion on changes in activities or allocations can be added to the modelling framework, corresponding to future land-use-climate policies. The activity/transition probability matrices and the output coefficients (i.e. the entire parameterization of EFDM) can be fundamentally modified according to the needs of the analysis, which might not be true for all corresponding simulators without accessing their source code.

Further aspects for comparisons between evenand uneven-aged management systems
Our simulations do not include the effects of biotic or abiotic damages. However, this should not affect the comparisons of the simulated alternatives, unless the disturbances are more frequent in some of them. In this sense, different forms of uneven-aged management, such as continuous cover forestry, may be more problematic than even-aged management, especially if applied intensively (Nevalainen 2017). If the 1 3 occurrence and degree of damages were modelled in terms of the axes of the matrices used in this study (e.g. volume and age) and included as additional activity and transition probabilities, respectively, we could assess the impacts of disturbances based on the Markov chain model framework (see also related discussion in Vauhkonen and Packalen 2018).
The use of the 3% rate of return when simulating thinnings results in relatively large harvest removals, even when maintaining the legislative minimum basal area. An ecological perception of feasible management practices in the proportion of the MUCLs that are not strictly protected could differ from those simulated. A more diversified set of management or natural interventions could have been considered here to mimic forms of uneven-aged management other than continuous cover forestry (Puettmann et al. 2015;Pukkala 2016a). Our simulations could also be extended to cover the responses to different harvests (Mehtätalo et al. 2014;Montoro Girona et al. 2017;Bose et al. 2018; see also Roessiger et al. 2016) or improved silvicultural performance due to operations, such as ditching, fertilization and the use of improved genetic material for regeneration (cf., Hynynen et al. 2015;Heinonen et al. 2018).
We did not include any ecological measure similar to those related to carbon and harvests. The additional set-aside area could be used as such a measure, because ecological values often benefit from no management (e.g. Sutherland et al. 2016). However, this conclusion is not exclusive for provisioning services or even all habitat services (cf., Peura et al. 2018). Our analyses could be extended using speciesspecific habitat suitability indices as in Pukkala (2016b) and Peura et al. (2018). Overall, a specific management system might not benefit all services and, moreover, the benefits could be site dependent (Biber et al. 2015;van der Plas et al. 2018), although alternative management may be focused on forests based on totally different criteria than suitability for a specific site. For example, continuous cover forestry could be a logical choice for poorly productive forests because of the low economic profitability of even-aged management systems. Similar scenario analyses, as carried out in this study, could be run to compare different management systems on specific sites, such as poorly productive, drained peatlands in Finland (cf. Nieminen et al. 2018), and this could be carried out with respect to a wider selection of ecosystem services and choices related to management regimes than considered in our study.

Concluding remarks on the implications of large-area forest management
The simulated large-area shifts from conventional evenaged management to alternative silvicultural systems revealed interesting development patterns that cannot be directly deduced from studies that are based on individual forest stands or holdings and upscaled to larger areas. At the national scale, the simulated development of carbon storage in the above-and below-ground living biomass, harvest removal, and harvesting costs differed depending on the forest type and the size of the area assumed to diverge from the business-as-usual management approach. Our results would suggest that it could be beneficial, with respect to the national wood supply, if large-scale adoption of the alternative management practices initially commenced in the less productive forests.
The forests selected to transit from even-aged to alternative management according to the higher values of the consval index (i.e. MUCLhigh, MULhigh strategies) would provide large amounts of timber as FAWS. When assigned to alternative management systems, harvesting becomes less efficient in terms of both removals and costs, because a portion of the trees remain in the forest and the trees that are removed must be acquired through selective harvests. In addition, the carbon storage in the above-and below-ground living biomass increased most when the shifts between management systems were realized according to the MUCLlow strategy. In this strategy, the forests that shifted to continuous cover forestry were generally less densely stocked, less fertile and had a low number of species, i.e. the selection emphasized the low values of the conservation value proxy. However, one-third of every p% selected by this strategy was always set aside and this proportion was emphasized in forests that were densely stocked, fertile and had a higher number of species, i.e. highest values of the conservation value index within the selected p%. With such an allocation, the forests that were set aside effectively increased carbon storage. The effects of harvest removals and costs were minor, because the majority of the most densely stocked FAWS remain for wood production in evenaged management rotations. On the other hand, the amount of thinnings from above increased towards the end of the simulation, i.e. when the growth of the initially less densely stocked forests permitted an increasing number of harvests; exhibited as a slight increase in harvesting costs at the end of the simulation. A period of altogether 50 years was simulated; if the simulations were to be extended beyond 2060, it would lead to an increase in other uncertainties that should be considered in separate analyses.
Analyses of our results from the point of view of different decision makers can provide further insights, illustrated here by three examples: • A MUCL alternative with p = 33% would correspond to the "third-of-third rule-of-thumb" (Hanski 2011), according to which one-third of the land area should be managed as a "multi-use conservation landscape" and one-third of this proportion be strictly protected. Even if this alternative could meet the conservational aims as reasoned 1 3 by Hanski (2011), the solution is inefficient with respect to carbon storage or harvesting. Solutions that are more feasible in aspects other than conservation could be found by examining the data in Figs. 1, 2, 3, 4, and 5 by means of multi-criteria decision analyses or corresponding tools involving the preferences of the different stakeholder groups. We did not optimize forest management for any of the objectives mentioned in this paragraph. • In addition to the simulated shifts of FAWS to FRAWS or FNAWS, Figs. 1, 3, 4, and 5 include a reference projection, where all forests are managed as FAWS. This demonstrates that about 20% of the forest area studied here is already subject to restrictions that prevent the use of heavy harvesting operations, such as final felling. If additional enforced restrictions were considered, it could be beneficial to compare their favourable effects with those already achieved under the current restrictions (cf., Sect. 3.3.). • Studies based on individual forest stands or holdings often present uneven-aged forest management alternatives as highly attractive to forest owners due to the potential to obtain higher profits by avoiding regeneration costs and by reducing capital costs through the removal of the largest trees in selective harvests. If alternative management systems are adopted over large areas and similar harvesting totals (as realized in the past; cf. MAF 2015) are still expected, the harvesting becomes inevitably more costly with selective harvests than final fellings. It is not clear as to who will pay the increased costs of wood procurement. The proportion of forests managed under different silvicultural systems is likely to depend on the resulting market mechanisms that probably differ depending on whether the shifts between management systems are enforced or voluntary.

Appendix 2
The parameterization of the EFDM requires that NFI data are classified into factors that represent the initial forest state. Transition probabilities, management activities, and output coefficients are specified for each factor combination (Vauhkonen and Packalen 2017). The present study differs in that more than one initial state was defined for the simulations by transferring different proportions of FAWS to the other wood availability categories according to a conservation value proxy (cf. Sect. 2.1). This appendix aims to (numerically) exemplify the computations related to the conservation value (1), classification of the NFI data to obtain the initial state (2)

Conservation value proxy
The proportion of land transferred from FAWS to the other categories was selected according to a conservation value proxy (adapted from Lehtomäki et al. 2015): where diameter and volume are the species-specific mean diameter and growing stock volume of a plot, f() is a speciesspecific transformation function used to convert the diameter to a conservation value index based on expert knowledge, w is a weighting that corresponds to site fertility, and sp is a species index. To compute consval, the trees measured from each NFI plot were assigned to groups of pine, spruce, birch, and other deciduous species. The median and maximum values of the mean diameters of these species were computed for each of the 15 forestry regions in Finland. Species-specific asymptote and scaling parameters (specifically, mod_asym, and parameters with suffices lavdia and ravdia from the file parameters-esmk.csv; Lehtomäki 2015) were related to the median and maximum values (as carried out in file gis.calculate.index.R; Lehtomäki 2015). The parameters produced sigmoidal transformation functions (Fig. 6) that were used to transform the species-specific mean diameters to conservation indices between 0 and 1. Figure 6 shows the species-specific mean diameters of the example plot, the resulting conservation index values and sigmoidal transformation functions for every forestry region and species. Smaller diameters yield greater transformed values in deciduous species than in coniferous species. For birch, this transformation was strongly forestry region specific, in that occurrences of large birch trees yielded a higher conservation value in regions where the median diameter of birch is low. The same species-specific weightings for site fertility, as used by Lehtomäki et al. (2015) in their "Coarse with classes" workflow for conservation prioritization based on multi-source NFI data, were extracted. The example plot above was assigned site weightings of 3.0, 4.0, and 7.0 for spruce, birch and aspen, respectively, as it was located on a Units in m 3 /ha for volume, 1/ha for stem number, and years for age  6 Sigmoidal transformations from species-specific mean diameter to conservation value. Black, red, blue, and green lines show the transformation functions for pine, spruce, birch, and other deciduous trees, respectively, and the coloured dots show how the trees measured from the example plot related to the transformation functions. Separate lines are drawn to depict the different forestry regions. Note that for some species, the transformation is determined piecewise, with different asymptote and scaling parameters applied for diameters below and above the region-specific median value. (Color figure online) 1 3 highly fertile soil. To obtain the final consval for the plot, the transformed diameters were multiplied by species-specific total stem volume and weighting for site fertility and summed over the species in a plot. The total consval value (907.12) for the example plot was thus obtained as a sum of Eq. 1 applied per species as follows: According to the example above, the conservation value would be lower for plots with the same volume but with less species. Site weightings affect the result considerably; the weightings are reduced to 1.0-1.5 for conifer species on less fertile soils. The weightings are also reduced for deciduous species, although less severely, and the occurrence of deciduous species in the least fertile soils are also weighted slightly higher. Overall, the consval index was found to reflect well the growing stock, species and site variation, relative to variation within forestry regions (Fig. 6).

EFDM input data
The example plot represents an area of 350 ha when computed according to NFI methodology for forest area estimation (cf. Vauhkonen and Packalen 2017). It is classified in the initial forest area distribution matrix as a spruce-dominated, private FAWS in the highest taxation class according to the static factors (Sect. 2.2.1). It is further classified to volume class #12, stem number class #6, and age class #17 according to its total volume (396 m 3 /ha), stem number (730 stems/ha), age (78 years), and class limits (Appendix 1). A similar classification is applied for each plot in both the full and pairwise data, and these data sources are used to derive the initial state, activity/transition probabilities, and output coefficients as described in Sect. 2.2 and in more detail by Vauhkonen and Packalen (2017).

EFDM simulations
The forest area distribution after one simulation step is obtained by two matrix multiplications (see Sect. 3 in Packalen et al . 2014): (1) the current area distribution is multiplied by the activity probabilities, which yields separate matrices that represent the area of each cell affected by each activity; and (2) the aforementioned intermediate matrices are multiplied by the transition probabilities, which yields the area distribution in the next step. Figure 7 depicts the Fig. 7 The proportion of different management activities applied to the example plot on the first step of the simulations. The activity probabilities vary depending on the choice between businessas-usual allocation (A BAU ) or schoolbook-allocation (A SB ) and whether the forest is considered as Forests Available for Wood Supply (FAWS) or Forests with Restrictions on Availability for Wood Supply (FRAWS) and, therefore, simulated using agevolume or stem number-volume matrices, respectively 1 3 activity probabilities that are applied to the example plot, depending on whether the plot is considered either as FAWS or FRAWS and managed according to either A BAU or A SB . Finally, the proportions are adjusted by two alternative total harvesting amounts, but these effects are omitted here for clarity purposes.
The area managed by thinnings or final fellings produces harvested removals. According to the output coefficients for the plot in the first simulation step, final felling of the plot described above would yield in total 200.5 m 3 /ha and 162.8 m 3 /ha of logs for sawn wood and pulp wood, respectively, at an estimated (total) cost of 8.76 €/m 3 . The corresponding values are 57.8 m 3 /ha and 89.1 m 3 /ha for thinning from below (total cost of 15.6 €/m 3 ) and 129.9 m 3 /ha and 95.8 m 3 /ha (total cost 12.3 €/m 3 ) for thinning from above.
The transition probabilities shift the initial area of the plot to multiple matrix cells. Using age-volume classes, the forest that is not managed gains age, but remains in the (highest) volume class. Thinnings from below transit the area to volume classes from #7 upwards, the most common class being #10. Using stem number, about half of the volume remains in the same stem number and volume class. Harvests shift the area from stem number/volume classes #3/#6 upwards, the most common target class being #4/#7. Notably, natural processes with stem number may shift the area both up and down in terms of classes, and shift the area downwards from both thinning from below and above, although based on different thinning rules. In the next simulation step, the forest area of the receiving classes is further updated using the activity and transition probabilities of the specific classes. Figure 8 shows how the initial area of 350 ha in one class evolved during the 10 simulation steps based on either agevolume and stem number-volume classes. Because of the multiple transitions that take place due to activity and transition probabilities, it might not be feasible to track the development of forests in a single class, but the overall development of the class structure is deemed realistic with respect to the properties of the forest.  Fig. 8 Distribution of age-volume (above) and stem numbervolume (below) classes after simulating 10 time-steps and applying class-specific activity and transition probabilities to the initial 350 ha of forest. Refer to Appendix 1 for the interpretation of class numbers on the x-axes. The grey tones within the bars depict volume classes, the darkest grey referring to the lowest number of volume class and vice versa