Searching for Pareto Fronts for Forest Stand Wind Stability by Incorporating Timber and Biodiversity Values

: Selecting a variant of forest regeneration cuttings that would ensure fulﬁlling multiple, frequently conﬂicting forest functions is a challenging task for forest management planning. The aim of this work is to present an e ﬃ cient and complex analysis of the impact of di ﬀ erent forest management scenarios on stand wind stability, timber production (economy), and biodiversity of a secondary mixed temperate forest in Central Europe. We evaluated four di ﬀ erent harvest-regeneration systems: clear-cutting, shelter-wood, selection cutting, and no-cutting using theSIBYLA growth simulator. We simulated forest stand development over time and applied 450 variants of 4 harvest-regeneration systems. The selected outputs from the simulator were used as indicators of the fulﬁlment of wood-production and non-wood-production functions. The calculated indicators were forest stability (height / diameter ratio), economic e ﬃ ciency (soil expectation value, SEV), and tree species diversity (Shannon index). These indicators were used as inputs for multi-criteria a posteriori decision analysis using the weighted summation method and Pareto fronts. The results revealed substantial trade-o ﬀ s among the three investigated criteria. The decision space was highly sensitive to their weighting system and included all regeneration systems. The Pareto fronts for wind stability revealed that the maximum stability could be achieved with shelter-wood based on target diameter. This variant, however, fulﬁls the other two examined functions only to a limited extent (SEV and diversity only to 9% and 27% of their absolute maxima). Other similar variants achieve high stability by sacriﬁcing the diversity and increasing SEV, simultaneously. If a high diversity level is favoured, optimal stability could be achieved by the selection system. The proposed approach enables objective testing of a large number of variants, and an objective assessment of stand management planning since it provides us with the complex multi-dimensional picture about the impact of criteria weights on the selection of optimal variants, and the relative fulﬁlment of individual criteria.


Introduction
In the past, European forests were primarily used for timber production, especially to fulfil demands for fuel and construction timber [1]. To maximise timber production, natural forests were frequently converted into pure and even-aged stands with changed tree species composition [2] and uniform structure, which nowadays cover around 6 million hectares of Europe [3]. In many countries, such secondary stands were usually managed using approaches that maximise profits from timber and minimise costs. These stands have recently become affected by natural disturbances, mainly windthrow and bark beetle attacks. Thus, management of such stands that ensures their sustainability from the point of view of production, stability, biodiversity, and climate change impacts has recently become a highly discussed issue [4]. The difficulty of selecting optimal management for such secondary even-aged forest stands is even more pronounced due to the problems that arise from their changed structure and tree species composition [5].
Currently, forest ecosystems are perceived as multifunctional systems providing a wide range of different products and services [6]. Selecting optimal harvest-regeneration methods for individual stands that can ensure multiple forest functions is a key management problem in many parts of the world [7], because forest management should meet a number of economic, social and ecological objectives, which are often contradictory [8]. For this study we adopted a definition of the "forest stand" term from the Slovak Forest Act 326/2005, in which a "forest stand" is defined as a basic spatial forest unit for forest monitoring, management planning, recording and control of applied forest management treatments with the area of at least 0.5 ha, and borders specified with regard to forest ownership, site and terrain conditions. Considering the aforementioned facts, an optimal forest stand size in the conditions of Slovakia is assumed to be 10 ha, but it can vary from 5 to 20 ha [9].
Forest potential to fulfil a specific function is frequently not achieved because management applied to a particular stand is a compromise between different demands of individual stakeholders (e.g., forest owners, governmental and non-governmental organizations of forest protection). Therefore, sustainable management of forest stands that ensures multiple forest functions requires integration of a wide range of data, information, models and methods. Multi-functional forest management planning representing an operational basis of sustainable forest management is a very complicated procedure because of the multiple dimensions it should cover [10]. A properly designed decision support system could facilitate this process, contribute to the transparency and quality of decision making and minimise conflicts in the management of natural resources [11]. At present, there are at least 100 computer decision support systems (DSS) used in many countries all over the world [12]. They are usually composed of several distinct parts, each responsible for a different task in the decision-making process including multi-criteria decision making (MCDM), simulation, optimisation, statistical evaluation, economic analysis, and management of information [13]. The overview of existing MCDM tools with regard to forest management planning is presented by [14].
One of the most widely used MCDM approaches is the weighted summation method that realizes the optimal management solutions based on the weights assigned to different criteria and sums to 1 (100%). These solutions explore the entire decision space and fail to analyze the trade-offs between criteria. Pareto fronts that belong to so-called a-posteriori multi-criteria decision analysis methods [15] define a set of non-dominated Pareto optimal solutions, i.e., solutions, for which an improvement in one objective causes a worsening of at least one other objective. By applying the Pareto front to the analyzed solutions we can obtain fast and easy to understand displays of combinations of several indicators. This can help us to examine efficient trade-offs between them, and account for a balanced combination of interests. The method has been applied in a number of studies dealing with environmental [16] and forest-management planning [15,17,18]. Several works [19,20] showed that the interactive approach can increase the efficiency of the planning process, because it may facilitate the involvement of multiple stakeholders in the decision-making process since the levels of achievements of various pre-defined objectives are easily quantified [15]. The abilities to account for multiple objectives and to involve stakeholders in decision making have been identified as highly required features of DSS [13]. The Pareto front links points that cannot be improved in one criterion without reducing the value of at least one of the remaining criteria. Thus, this method provides the information needed to assess the compromises in forest management [21].
In the conditions of Slovakia, forest management for a single forest stand is prescribed based on the generalised rules and measures derived for a so-called operational set of stands [22]. Hence, no management alternatives are given or examined a priori with regard to individual forest functions. Moreover, such a static approach does not account for the ongoing environmental changes, which should be considered in efficient forest management planning. The frequency and intensity of natural disturbances has been increasing over the last decades [23]. In Slovakia and the whole of Central Europe, wind is the main disturbance factor affecting forests of the region, and its impact is likely to increase further in future due to climate change [23,24]. Thus, selected forest management should lead to a higher static stand stability that would reduce forest susceptibility to windthrows and ensure future wood production. Decision analysis tools can provide a valuable support for multi-functional forest management planning because they can objectively quantify and compare the effect of multiple management possibilities [10]. Moreover, such tools can increase flexibility and adaptability of the forest management planning process to specific site and stand conditions.
The aims of this contribution were to (1) present a decision analysis tool based on the SIBYLA growth simulator [25], and a posteriori multi-criteria decision making techniques, (2) show its functionality by solving a stand-level planning problem (optimal harvest-regeneration system of secondary forest stands in central European), (3) analyse the impact of planning variants on three criteria: wind stability (height to diameter ratio), economic efficiency (Soil expectation value, SEV) and biodiversity (Shannon-index), and explore the effects of different combinations of weights of three criteria on the selection of an optimal harvest-regeneration variant for an individual stand, and (4) analyse the trade-offs among criteria and reveal the Pareto fronts for wind stability of the stand susceptible to wind disturbances. Secondary forests are currently most susceptible to disturbances and environmental changes [4]. Hence, optimal management of such stands is crucial to ensure their sustainability and multiple forest functions. The proposed approach of data analysis enables objective testing of a large number of harvest-regeneration variants, and an objective assessment of stand management planning since it can provide us with the complex multi-dimensional picture about the impact of criteria weights on the selection of optimal harvest-regeneration variants, and the relative fulfilment of individual criteria.

Decision Analysis Tool
We developed a decision analysis tool that could solve the stand-level forest management decision task, which of the harvest-regeneration variants can simultaneously best fulfil potentially conflicting management objectives. The tool requires specification of site and forest conditions, harvest-management variants and criteria that represent individual management objectives. In our case we selected the following three objectives: (1) the economic value of timber production, (2) stand diversity, and (3) stand stability, because they are important from the point of forest sustainability and economy, and they are all quantifiable. Based on the previous studies, e.g., [26], we can expect conflicts between timber production and non-timber objectives, while stand diversity and stand stability should be complementary. Information about the impact of different management variants on stand development is obtained from the simulations performed with the SIBYLA growth simulator. This simulator was selected because it was calibrated for central European tree species, which allows us to simulate stand dynamics using species-specific models, is able to simulate a great number of management treatments, and contains a 3D forest visualisation tool [27]. In the next step, the simulated output is processed to quantify the selected objectives over the whole simulated period (see Section 2.5 for more details). The processed output was used in multi-criteria analysis performed with the OPTIMUS software specifically developed as a part of the presented decision analysis system. OPTIMUS is a tool that enables user-specified optimisation based on multiple criteria selected by a user, different optimisation techniques (conjunctive and disjunctive method, PRIAM (PRogramme utilisant l'Intelligence Artificielle en Multicritere), order, lexicographical, score method, weighted summation, basic variant, analytic hierarchy process) and different approaches of weight calculation (proportional, manual, order, Fuller, Point, Saaty, Gradual distribution of weights) for selected criteria [28]. The output calculated by OPTIMUS was used for the construction of Pareto fronts and graphical interpretation of the results that was performed in R in order to present the results in a user-friendly and easy-to-understand design. The use of the proposed decision analysis tool was demonstrated employing the information about a typical Central European secondary forest stand dominated by Norway spruce. Such secondary forests are currently most susceptible to disturbances and environmental changes [4]. Hence optimal management of such stands is crucial to ensure their sustainability and multiple forest functions.

Case Study Area
For the purposes of this study we used the data from a secondary mixed forest stand located in the Slovenské Rudohorie Mountains of the Western Carpathians, Slovakia ( Figure 1). The forest is a representative secondary mixed forest stand in Slovakia and in Central Europe, in which Norway spruce (Picea abies (L.) Karst.) is the dominant species (Table 1). Northern aspect was prevailing within the stand. The altitude ranged from 430 m to 470 m a.s.l. The bedrock is formed by phyllites, slates, metapsammites, metavolcanic rocks, and sporadically metacarbonates. The soil is acidic to strongly acidic, and the dominant soil type is podsolic Cambisol. The annual average temperature is 8.7 • C and the average annual precipitation total is 605 mm. optimisation techniques (conjunctive and disjunctive method, PRIAM (PRogramme utilisant l'Intelligence Artificielle en Multicritere), order, lexicographical, score method, weighted summation, basic variant, analytic hierarchy process) and different approaches of weight calculation (proportional, manual, order, Fuller, Point, Saaty, Gradual distribution of weights) for selected criteria [28]. The output calculated by OPTIMUS was used for the construction of Pareto fronts and graphical interpretation of the results that was performed in R in order to present the results in a user-friendly and easy-to-understand design. The use of the proposed decision analysis tool was demonstrated employing the information about a typical Central European secondary forest stand dominated by Norway spruce. Such secondary forests are currently most susceptible to disturbances and environmental changes [4]. Hence optimal management of such stands is crucial to ensure their sustainability and multiple forest functions.

Case Study Area
For the purposes of this study we used the data from a secondary mixed forest stand located in the Slovenské Rudohorie Mountains of the Western Carpathians, Slovakia ( Figure 1). The forest is a representative secondary mixed forest stand in Slovakia and in Central Europe, in which Norway spruce (Picea abies (L.) Karst.) is the dominant species (Table 1). Northern aspect was prevailing within the stand. The altitude ranged from 430 m to 470 m a.s.l. The bedrock is formed by phyllites, slates, metapsammites, metavolcanic rocks, and sporadically metacarbonates. The soil is acidic to strongly acidic, and the dominant soil type is podsolic Cambisol. The annual average temperature is 8.7 °C and the average annual precipitation total is 605 mm.  Note: * mean of representative samples according to [29], ** height in 100 years [30].
The forest stand was managed according to the usual forestry practice of less-intensive thinning interventions from below (mostly the removal of damaged and low-quality trees). At the time of the measurement, the forest stand was 60 years old and the stocking was 0.8. The stand characteristics were assessed on six circular plots using systematic sampling. The Stem diameter at breast height  Note: * mean of representative samples according to [29], ** height in 100 years [30]. The forest stand was managed according to the usual forestry practice of less-intensive thinning interventions from below (mostly the removal of damaged and low-quality trees). At the time of the measurement, the forest stand was 60 years old and the stocking was 0.8. The stand characteristics were assessed on six circular plots using systematic sampling. The Stem diameter at breast height (DBH) was measured using a calliper with precision of 1 cm. Tree height (h) was measured using a hypsometer Vertex. The basic descriptive characteristics of the stand are shown in Table 1.
Additional data needed as input for the SIBYLA growth simulator are defined in the text below. The soil nutrient supply of the case study area was 0.8 and the soil moisture was 0.6 (a relative value in the range from 0 to 1; where 0 is infertile or dry and 1 is fertile or wet). From the climatic point of view, the stand belongs to a slightly warm, and a slightly moist hilly to upland area with the annual temperature amplitude (i.e., difference between the average monthly temperature in the warmest and coldest months) of 15.4 • C, daily mean temperature in the vegetation period (i.e., from April to September) of 14 • C and precipitation total in the vegetation period of 600 mm. The vegetation period in the selected stand lasts 160 days with the daily mean temperature above 10 • C. The concentration of nitrogen oxides (NO x ) and carbon dioxide (CO 2 ) in the air was equal to 307.8 ppb (parts per billion) and 354.8 ppm (parts per million), respectively.

Harvest-Regeneration Variants
In the presented paper, we used most frequent forms of four harvest regeneration systems [31]. Their rotation and regeneration periods were defined with regard to tree species composition of the case study area. In total, 450 different harvest regeneration variants were defined and their impact on the selected stand characteristics (Section 2.4) was simulated with SIBYLA. These variants represent the analysed decision space. The particular forms and the variants that were applied in this study are described in Table 2. According to the Slovak Forest Act 326/2005, during the period of forest regeneration a forest stand is divided into so-called regeneration elements, inside which a specific cutting is performed at a particular time (see Figure 2). For example, in the case of small-scale clearcutting, the area cut at a certain time point must be less than 3 ha, the width of the cut area must be less than the double of the mean stand height, and the distance between the two cut regeneration elements must be at least equal to the width of one regeneration element (Forest Act 326/2005). Hence, the cut area is not equal to the stand area, but it depends on stand dimensions and the type of the applied management (see Figure 2 for examples).

SIBYLA
SIBYLA is a growth simulator that belongs to the group of individual tree distance-dependent empirical models [25,32]. SIBYLA is based on the modelling principle and the algorithms of SILVA 2.2 [33], some parts of which were modified or completely revised. SIBYLA was calibrated using the data of about 9213 beech trees, 7358 spruce trees, 3444 oak trees, 1181 pine trees, and 1137 fir trees [27].
In our study, we generated the data about individual trees within a 16 ha large forest stand using the basic stand characteristics (Table 1) of the case study area described in Section 2.2., because no information on individual tree positions and dimensions was available. Hence, the generated forest stand was 60 years old, and its tree species composition, mean diameter, height and volume of individual tree species were in accordance with Table 1. The stand was divided into 16 regeneration elements, each of 1 ha (200 × 50 m). The sizes of the generated forest stand, and of the regeneration element were chosen with regard to technical and simulation possibilities of the growth simulator and to ensure that all predefined harvest regeneration variants were performed to their full extent (i.e., their regeneration periods were completed, and all planned cuttings and phases were performed during the simulation). Hence, the area of the generated forest stand was set to 16 ha to ensure that at the end of the regeneration period the whole stand is harvested with the most complex harvest-regeneration variant, which was an expanding small-scale shelter-wood with 3 cutting phases per decade ( Figure 2A). In the case of simpler variants, several regeneration elements were used as starting points to fulfil variant  Figure 2B). We simulated the development of the generated forest stand under the management of 450 predefined harvest regeneration variants (Table 2). Note: * target diameter was determined on the base of natural conditions; ** single tree cutting according to Liocourt model of selection system; *** self-development. Each row represents one decade of the regeneration period and applied cuttings within a 16 ha large forest stand divided into 16 regeneration elements. Coloured stripes in A represent cutting phases of the shelter-wood system, while vertical stripes of different color represent preparatory and seeding cuts, and horizontal stripes represent final cut of the respective regeneration element. Red and green colored elements in B represent stand parts, where first and second clearcuttings took place in 1st and 5th year of the decade, respectively.

SIBYLA
SIBYLA is a growth simulator that belongs to the group of individual tree distance-dependent empirical models [25,32]. SIBYLA is based on the modelling principle and the algorithms of SILVA 2.2 [33], some parts of which were modified or completely revised. SIBYLA was calibrated using the data of about 9213 beech trees, 7358 spruce trees, 3444 oak trees, 1181 pine trees, and 1137 fir trees [27].
In our study, we generated the data about individual trees within a 16 ha large forest stand using the basic stand characteristics (Table 1) of the case study area described in Section 2.2., because no information on individual tree positions and dimensions was available. Hence, the generated forest stand was 60 years old, and its tree species composition, mean diameter, height and volume of  Note: * target diameter was determined on the base of natural conditions; ** single tree cutting according to Liocourt model of selection system; *** self-development.

Indicators Selected for the Multi-Criteria Decision-Making (MCDM) Process
Based on the literature survey we selected three indicators for the multi-criteria decision-making process of selecting a harvest-regeneration system in the mixed forest that would best fulfil environmental, ecological and production forest functions. The indicators were derived from standard outputs of SIBYLA.

Non-Production Function
• Stand stability As a measure of stand stability we selected height to diameter (H/D) ratio, which is the most commonly applied indicator of stand stability [34,35]. The values of this measure increase with decreasing stability, i.e., low values indicate high stability [36]. Therefore, we subtracted the calculated values of H/D from 1 to obtain the increasing trend of stand stability with increasing values of the criterion.

•
Stand diversity A great number of species diversity indices can be found in the literature on biological diversity and ecological monitoring. Here we applied a common 'Shannon's' Index or 'H' index [37], which is calculated according to [38]: where p i is the proportion of basal area of ith species in the stand, and N is the number of species in the stand.

Production Function
Soil expectation value (SEV).
where NPV is the net present value, r is the discount rate (2%), and T is the rotation period. Net present value is characterised as the present value of revenues minus the present value of costs. Net present value defines an investor's willingness to pay for an asset based on estimated benefits, costs, and the desired rate of return. Thus, NPV is a powerful tool in valuing forest properties [39]. Symbolically, NPV is: where R stands for the revenues, y is year, r is the discount rate (2%) and C is the costs. In our case, we used 2% as the most commonly used discount rate in Central European forestry analyses [40].

Multi-Criteria Optimisation
We applied two a posteriori multicriteria decision making techniques: (1) the weighted summation method and (2) the Pareto front technique to examine the effect of different harvest-regeneration variants on forest stability, diversity and economic value.

Weighted Summation Method
The weighted summation method is one of the most widely used MCDM approaches [15]. The method requires standardisation of values of included criteria based on their minimum and maximum values from all tested variants. The overall performance measure of a specific variant is calculated as follows: where j is a criterion, i is the number of criteria used in the analysis, HV j is the value of a criterion j for a specific variant, MIN(K j ) and MAX(K j ) are the minimum and maximum values of a criterion j obtained from all tested variants, respectively, and VK j is the weight of a criterion j.
Since under specific weights this method converges to a single point solution, in our analysis we tested 63 different combinations of weights of the three selected criteria, starting from maximum weights (i.e., 1) of a single criterion, e.g., stability, when the weights of the other two criteria were equal to 0, up to proportional weights of all three criteria (0.333-0.333-0.334). Such an analysis can enlighten how the decision making is affected by different weights of individual criteria. In this study we applied the weighted summation method and different combinations of weights using the OPTIMUS tool [28].
Sensitivity of the selection of optimal variants was examined with a statistical test of the differences in the relative values of the indicators obtained under the optimal variant and other variants. We accounted for the theoretical error from the repetition of the simulations due to randomised processes implemented in SIBYLA. The calculation of the error of the relative value of the indicator was driven by a binomic distribution [41]. Variants were identified as alternative insignificant variants if their confidence intervals overlapped with the confidence interval of the relative values of the indicators for the optimal variant. Statistical tests were performed at 95% significance level.

Pareto Front
The Pareto front is one of multi-criteria decision making (MCDM) techniques, which helps users to analyse trade-offs between different criteria [21]. It is based on a posteriori preference analysis that facilitates the specification of the levels of achievement of various objectives in a typical multi-objective forest management planning framework [42]. The Pareto front links points that cannot be improved in one criterion without worsening the value of at least one of the remaining criteria. Thus, this method provides the information needed to assess the compromises in forest management [21].
We constructed bi-criteria and multi-criteria Edgeworth Pareto fronts as envelope curves of all simulated variants by linking all non-dominated solutions that maximise the pairs and the triplet of criteria in the analysed decision variable space. First, we created two-dimensional graphs, where each point of the global Pareto front indicates an efficient combination of levels of achievement of two criteria shown in the graph. A decision maker can examine the relationships between the pairs of the selected indicators, i.e., in our case stability, stand diversity index, and SEV, which are displayed on the vertical and horizontal axes. Furthermore, we created a three-dimensional graph showing all three criteria simultaneously. In this graph, each point of the global Pareto front presents a combination of the pair of criteria shown in the axes, and the third criterion is represented by local Pareto fronts bounded by the corresponding values.
The analyses were conducted in R [43] and figures were produced using ggplot2 [44] and ggtern [45] packages.

Results
The structure of the developed decision analysis tool is presented in Figure 3. The SIBYLA growth simulator [25] was used for simulating the growth of a forest stand managed by selected harvest regeneration variants ( Table 2). The selected outputs produced by SIBYLA were taken to calculate the values of the pre-defined indicators, which represented the input to OPTIMUS that performed the multi-criteria evaluation based on the weighted summation method.

Results
The structure of the developed decision analysis tool is presented in Figure 3. The SIBYLA growth simulator [25] was used for simulating the growth of a forest stand managed by selected harvest regeneration variants ( Table 2). The selected outputs produced by SIBYLA were taken to calculate the values of the pre-defined indicators, which represented the input to OPTIMUS that performed the multi-criteria evaluation based on the weighted summation method. The two-dimensional relationships between the selected indicators can be examined using the bi-criteria Pareto fronts (Figure 4). They showed that the stability of the forest stand and the diversity of the forest stand according to Shannon (H') were complementary objectives and competed only at high stability levels ( Figure 4B). The achievement of 100% diversity was possible at stability equal to 75%. Hence, the forest stand can be managed to achieve high levels of diversity and stability at the same time. On the other hand, diversity and SEV ( Figure 4A) as well as stability and SEV ( Figure 4C) were more competing objectives. To achieve 100% SEV, the diversity index would be reduced to 18%, and stability to 31%.
Using the WSM and weight combinations, we obtained 63 variants with the highest score sums for specific weight combinations from all 450 tested variants. Out of them, variants based on the selection cutting system were most frequent, since they occurred 30 times among the best variants (i.e., 48% of best variants), particularly in the cases with high weights on diversity and stability. The second most frequent harvest regeneration form was the shelter-wood based on target diameter, which occurred 23 times (i.e., 37%), particularly in the cases when the weight on stability was high. Large-scale clearcutting system (8 times, i.e., 13%) was found to be among the 63 best ones, when the weight on production was high ( Figure 5).
When we examined optimal variants in more detail, we found that only 7 different variants out of 450 tested harvest regeneration variants occurred among the 63 best ones ( Figure 5, Table 3). Selection cutting with a target diameter of 65 cm and 1 target tree per hectare occurred most frequently, followed by the same management system based on target diameter of 70 cm (Table 3). These variants of selection cutting maximised the diversity criterion. Stand stability was maximised The two-dimensional relationships between the selected indicators can be examined using the bi-criteria Pareto fronts (Figure 4). They showed that the stability of the forest stand and the diversity of the forest stand according to Shannon (H') were complementary objectives and competed only at high stability levels ( Figure 4B). The achievement of 100% diversity was possible at stability equal to 75%. Hence, the forest stand can be managed to achieve high levels of diversity and stability at the same time. On the other hand, diversity and SEV ( Figure 4A) as well as stability and SEV ( Figure 4C) were more competing objectives. To achieve 100% SEV, the diversity index would be reduced to 18%, and stability to 31%. with a shelter-wood target diameter variant with a rotation period of 160 years and a regeneration period of 40 years. On the contrary, the variant maximising SEV (i.e., large-scale clearcut) was found to have a shorter rotation period of 90 years and a short regeneration period of 10 years (Table 3). The analysis revealed that the variant that fulfils SEV to its maximum can fulfil selected stability and diversity criteria only to 31% and 18% of their maxima. Similarly, variants based on target diameter selection that maximise stability are able to fulfil SEV only to 9% and diversity to 27%. Diversity can be maximised with variants based on the selection cutting system, which can fulfil SEV and stability to approximately 12% and 75% of their maxima, respectively (Table 3, Figure 6).   Using the WSM and weight combinations, we obtained 63 variants with the highest score sums for specific weight combinations from all 450 tested variants. Out of them, variants based on the selection cutting system were most frequent, since they occurred 30 times among the best variants (i.e., 48% of best variants), particularly in the cases with high weights on diversity and stability. The second most frequent harvest regeneration form was the shelter-wood based on target diameter, which occurred 23 times (i.e., 37%), particularly in the cases when the weight on stability was high. Large-scale clearcutting system (8 times, i.e., 13%) was found to be among the 63 best ones, when the weight on production was high ( Figure 5).  Table 3 for an explanation of abbreviations) with regard to weights of individual criteria. Legend: CLS = large-scale clearcutting system, NoC = no cutting, STC = selection single tree cutting system, STD = shelterwood based on target diameter. Variants are described in Table 3. The sensitivity analysis showed that the maximum error of the sum of relative indicators at 95% significance level was 4.5%, i.e., maximum span of the 95% confidence intervals was ±4.5%. Within the calculated confidence intervals, a varying number of alternative variants was identified for different combinations of criteria weights. The results of alternative management variants that were insignificant to optimal variants are presented in Supplementary Materials. For seven weight  Table 3 for an explanation of abbreviations) with regard to weights of individual criteria. Legend: CLS = large-scale clearcutting system, NoC = no cutting, STC = selection single tree cutting system, STD = shelterwood based on target diameter. Variants are described in Table 3.
When we examined optimal variants in more detail, we found that only 7 different variants out of 450 tested harvest regeneration variants occurred among the 63 best ones ( Figure 5, Table 3). Selection cutting with a target diameter of 65 cm and 1 target tree per hectare occurred most frequently, followed by the same management system based on target diameter of 70 cm (Table 3). These variants of selection cutting maximised the diversity criterion. Stand stability was maximised with a shelter-wood target diameter variant with a rotation period of 160 years and a regeneration period of 40 years. On the contrary, the variant maximising SEV (i.e., large-scale clearcut) was found to have a shorter rotation period of 90 years and a short regeneration period of 10 years (Table 3). The analysis revealed that the variant that fulfils SEV to its maximum can fulfil selected stability and diversity criteria only to 31% and 18% of their maxima. Similarly, variants based on target diameter selection that maximise stability are able to fulfil SEV only to 9% and diversity to 27%. Diversity can be maximised with variants based on the selection cutting system, which can fulfil SEV and stability to approximately 12% and 75% of their maxima, respectively (Table 3, Figure 6). Table 3. List of variants with the highest sums of scores that occurred within 63 best variants identified for different combinations of criteria weights by the weighted summation method. Percentage values (% SEV, % Stability, and % H') were calculated as a ratio of the actual value of the criteria achieved under particular management variant to the maximum value from 450 variants. Forests 2020, 11, x FOR PEER REVIEW 13 of 20 Figure 6. A three-dimensional graph illustrating relationship between the relative values (%) of all three selected indicators simultaneously. Red dots refer to optimal variants (see Table 3 for an explanation of abbreviations) selected based on the weighted summation method for different combinations of criteria weights. Legend: CLS = large-scale clearcutting system, CSS = small-scale clearcutting, NoC = no cutting, SES = shelterwood expanding scale, SLS = shelterwood large scale, SSS = shelterwood small scale, STC = selection single tree cutting system, STD = shelterwood based on target diameter.

Discussion
The results of our decision aid tool that combined SIBYLA with the optimisation tool showed that the defined management objectives were competing (Figures 4-6), and therefore, it was not possible to determine a harvest regeneration variant that could ensure fulfilling all selected forest functions at 100%. Similarly to [26], stronger conflicts were found between SEV representing production function and non-production indicators (Shannon index or H/D ratio, Figures 4 and 6). Furthermore, we revealed the Pareto efficiency front for stability of the stand against wind and its trade-offs with diversity and SEV. This analysis could optimally allocate the examined harvestregeneration systems in the multi-criteria decision space and simultaneously reveal their effects on decision criteria. The Pareto curves deliver useful information for decision-makers to trade forest values in an optimal way. The proposed Pareto-optimisation approach can be adopted for other forest services based on the outcomes from SIBYLA or other simulation tools.
The analysis of our results showed that continuous-cover forestry approaches were identified as those that can best fulfil multiple criteria. Such management has been found to be more efficient in fulfilling ecosystem services than clearcutting [46,47]. In addition, continuous cover forestry was also found to increase carbon stocks in secondary forests under climate change scenarios [4]. However, our multi-criteria analysis revealed that the fulfilment of the production function was significantly reduced if variants based on single tree selection were applied to the stand of our interest, although their application could achieve high levels of both non-production functions (Table 3). Figure 6. A three-dimensional graph illustrating relationship between the relative values (%) of all three selected indicators simultaneously. Red dots refer to optimal variants (see Table 3 for an explanation of abbreviations) selected based on the weighted summation method for different combinations of criteria weights. Legend: CLS = large-scale clearcutting system, CSS = small-scale clearcutting, NoC = no cutting, SES = shelterwood expanding scale, SLS = shelterwood large scale, SSS = shelterwood small scale, STC = selection single tree cutting system, STD = shelterwood based on target diameter.
The sensitivity analysis showed that the maximum error of the sum of relative indicators at 95% significance level was 4.5%, i.e., maximum span of the 95% confidence intervals was ±4.5%. Within the calculated confidence intervals, a varying number of alternative variants was identified for different combinations of criteria weights. The results of alternative management variants that were insignificant to optimal variants are presented in Supplementary Materials. For seven weight combinations no insignificant alternative variants were found, while for the rest we found from 1 to 24 alternative variants. No alternative variants were revealed for the weight combinations, were one objective is strongly prioritized. Most commonly, we found 4 or 5 alternative variants to the optimal one. The number of alternative variants increases with the increasing balance between the weights of individual criteria. The maximum number of identified alternative variants was 24 for a combination of weights 20% diversity × 50% stability × 30% SEV. Although the analysis revealed that for 27 different weight combinations alternative variants represented different harvest-regeneration systems, in the other 36 cases the alternative variants were different variants of the same harvest-regeneration system. For example, in the case of weight on diversity of 50% or more, the prevailing number of alternative variants represented a single tree selection cutting system based on different target diameters and/or number of target trees. Similarly, for high weights on stand stability, the majority of both optimal and alternative variants were based on the target diameter shelterwood system.
For the mutual evaluation of all three indicators simultaneously, we created a three-dimensional graph, which presents all possible combinations of SEV, forest stand diversity and stability values ( Figure 6). This graph shows how individual criteria compete with each other in the whole examined space. As can be seen from Figure 6, most of the optimal variants selected by WSM for different combinations of weights occur at the outer Pareto front. The fulfilment of individual criteria is closely linked to the weights applied in WSM, e.g., the variant based on clearcutting system (CLS_90_10) maximises SEV and was found to be optimal if SEV weights were between 70% and 100%, i.e., weights of the other two criteria were low (0%-30% or 0%-20% for stability and diversity, respectively, see Figure 5). The maximum level of stability can be achieved with the single tested variant, namely shelterwood based on a target diameter with 160 year-long rotation period and 40 year-long regeneration period (Table 3), which is unlike other optimal variants not located at the outer Pareto front ( Figure 6) due to low fulfilment of other functions (9% of maximum SEV and 27% of maximum diversity). This variant was identified as optimal by WSM for 9 weight combinations and high stability weights from 60% to 100% ( Figure 5). High levels of stand stability exceeding 80% of absolute maximum can be achieved with other 17 variants of the shelterwood system based on target diameter (Figures 4 and 6). However, none of them occurred among the optimal variants selected by WSM ( Figure 5), because the application of variants fulfilling high levels of stability substantially reduces both diversity and SEV achieving less than 36% and 28% of their maxima (Figures 4 and 6). Optimal single tree-cutting variants (STC_65_1, STC_70_1, STC_75_1) primarily maximise diversity, but their application enables high levels of stability to be achieved exceeding 75% of its absolute maximum (Table 3, Figure 6) even with lower weights on stability ( Figure 5). However, these variants significantly reduce the fulfilment of the production function to 11% or 12% of its absolute maximum (Table 3).

Discussion
The results of our decision aid tool that combined SIBYLA with the optimisation tool showed that the defined management objectives were competing (Figures 4-6), and therefore, it was not possible to determine a harvest regeneration variant that could ensure fulfilling all selected forest functions at 100%. Similarly to [26], stronger conflicts were found between SEV representing production function and non-production indicators (Shannon index or H/D ratio, Figures 4 and 6). Furthermore, we revealed the Pareto efficiency front for stability of the stand against wind and its trade-offs with diversity and SEV. This analysis could optimally allocate the examined harvest-regeneration systems in the multi-criteria decision space and simultaneously reveal their effects on decision criteria. The Pareto curves deliver useful information for decision-makers to trade forest values in an optimal way. The proposed Pareto-optimisation approach can be adopted for other forest services based on the outcomes from SIBYLA or other simulation tools.
The analysis of our results showed that continuous-cover forestry approaches were identified as those that can best fulfil multiple criteria. Such management has been found to be more efficient in fulfilling ecosystem services than clearcutting [46,47]. In addition, continuous cover forestry was also found to increase carbon stocks in secondary forests under climate change scenarios [4]. However, our multi-criteria analysis revealed that the fulfilment of the production function was significantly reduced if variants based on single tree selection were applied to the stand of our interest, although their application could achieve high levels of both non-production functions (Table 3).
From the point of view of stability, the variant based on target diameter with long rotation and regeneration periods of 160 and 40 years (STD_160_40) was found to be optimal for the stand of interest in a multi-purpose decision space, although other optimal variants based on single tree selection could also fulfil stability at high levels of its absolute maximum ( Figures 5 and 6, Table 3). Selective approaches were found to have a positive effect on mechanical stability of trees [48,49], small-scale stand heterogeneity [50], and better acclimation to wind disturbances [51]. The remaining trees in the stand can use greater growing space, which is reflected in their fast volume increment [52]. Thick trees are considered more profitable, and the remaining prospective trees represent production capital [53].
By applying this regeneration system for a longer time, stands become gradually uneven-aged, which are considered more stable than even-aged forests [54].
Palmer [55] stated that a target diameter approach is better than selective thinning also from the point of the valourisation from cutting, because in the case of a target diameter approach the trees suitable for the wood-processing industry are cut, while in the latter case the selection depends also on tree positions [56]. The positive economic effect was however not confirmed by our results, which showed that the optimal variant based on target diameter (STD_160_40) fulfilled SEV only to 9%, while the optimal single tree-cutting variants (STC_65_1, STC_70_1, STC_75_1) fulfilled SEV to 11 or 12% (Table 3, Figure 6). Considering the criteria that we analysed in our study, namely height-diameter ratio and Shannon index, single tree selection variants seem to be more suitable to increase stand stability as a whole, because they maximise tree species diversity and keep mechanical stability at a high level (Table 3). This is in accordance with the so-called "diversity-stability theory" [49], according to which mixed stands are usually less affected by windthrows or other natural disturbances than pure stands [57,58].
It is very important that complex forest management problems, where multiple decision-makers and stakeholders are involved, are addressed by an objective approach [42]. Although we did not directly involve stakeholders in the presented study, the three selected criteria represent contradictory interests possibly prioritised by different groups. With the applied approach we can analyse the requirements of several decision makers at the same time by accounting for a whole range of weights of individual management objectives ( Figure 5). The visualisation of the whole decision space ( Figure 6) enables users to perform a thorough analysis of all examined management variants and to select solutions that would best meet the requirements interactively on the computer screen. These opportunities may facilitate communication, cooperation, negotiation, and information-sharing between stakeholders with various preferences and interests [12]. This will enable decision makers to evaluate individual management objectives, and the trade-offs between them, to specify the preferences or targets for different criteria, and to reach necessary compromises when making final decisions. The presented results demonstrated the use and the potential of the suggested decision analysis tool for the selection of an optimal harvest-regeneration system in a typical secondary forest stand of central Europe with unnatural tree species composition dominated by Norway spruce on the basis of different requirements of decision-makers. The outputs showed that the combination of the SIBYLA growth simulator [27] and the optimisation tool OPTIMUS provided us with results that were comparable with those from the complete DSSs (e.g., Sadflor [59], Anafore [60], Heureka [61]). Such a combination of simulators with optimisation or multi-criteria decision analysis tools has been identified as a promising trend to support forest management decision making because of increasing computing capacity [62], and the variety of possible simulated outputs [63]. The ability and flexibility of the SIBYLA simulator to model various management scenarios defined by a forest manager including different types of thinning and harvest regeneration systems from simple clearcutting systems up to the selection of individual trees [32] enabled us to investigate the impact of 450 different harvest-regeneration variants ( Table 2) on the selected objectives. The simulator is free software known to foresters in Central Europe, which has been designed with regard to the needs of practical forestry, i.e., it is a user-friendly software allowing a user to define its forest stand based on the readily available data from e.g., forest management plans. Since SIBYLA can produce a large range of different production, biodiversity, and economic outputs at a tree and stand levels [27], its application is not limited to the three stated goals of forest management planning. Moreover, the model is sensitive to climate [64], CO 2 and nitrogen following the model of ecological classification by [65], which are pre-requisites for evaluating alternative management strategies under climate change conditions [66]. The applicability of SIBYLA for forest growth modelling in Central Europe has already been documented in multiple works [22,[67][68][69][70] and its development is still in progress.
The presented application of the suggested decision aid tool showed that it can be used as an efficient support when selecting an optimal harvest-regeneration system for a particular stand, as it can examine and evaluate the impact of a great number of possible management alternatives according to multiple criteria. Moreover, this approach is more objective than the currently applied planning system in Slovakia, which is based on the classification of a stand into one of predefined so-called forest management models that determine main management goals, and specify methods, rules and instructions, on how to achieve the defined goals within a rotation period. Hence, the management prescribed for a particular stand greatly depends on this categorisation and is independent of the objectives of stakeholders [22]. In addition, the creation of forest management plans and management strategies can be supported by multiple visualisation possibilities of the forest state and growth in the SIBYLA simulator and of the fulfilment of the predefined management objectives in a multi-dimensional decision space. The optimal variant can be selected directly on the computer screen based on their preferences from the set of feasible goals [21]. In addition, by accounting for the uncertainty when selecting optimal variants, alternative insignificant variants can be provided (Supplementary Materials) and subsequently examined from the point of their practical applicability in specific terrain conditions. The whole approach of the proposed decision aid tool is flexible and enables stakeholders to change their management goals at any time of the rotation period. This is of particular importance in secondary forest stands, which are susceptible to environmental changes [4], and need to be managed adaptively to ensure their long-term sustainability and forest functions [66]. Such a solution using ready-to-use tools may solve the most urgent problems in forest management planning in the context of ongoing climate change and disturbances and serve as a starting point or as a basis for more sophisticated solutions in the future.
Since the provision of many ecosystem services requires planning at a large scale, it is also possible to apply the tool also at a level of landscape. Since SIBYLA is a model simulating individual forest stands, a specific landscape can be specified as a set of representative forest stands, for which we can propose optimal management variants with regard to multiple management objectives. Such an approach will allow the trade-off between the ecosystem services at both stand and landscape levels to be efficiently analysed. Combining the fulfilment of multiple ecosystem services is an integrative approach of forest management, which prevails in Central Europe [71] and has been shown to be more effective than segregation [72].
It should be noted that the presented results are affected by the interest rate used for the calculation of NPV [40,73]. The interest rate is actually the weight that determines how future costs and revenues flow compared to current costs and revenues until the valuation. An interest rate that is too low leads to economically inefficient intentions, while an interest rate that is too high leads to the rejection of effective intentions [74]. The determination of the interest rate is based on the determination of the nominal interest rate most often determined on the basis of the rates stated in the loan agreements. The real interest rate is calculated as the nominal interest rate, which is reduced by the value of the weakening of the borrowed or deposited amount during the loan or deposit period (inflation). There is also a risk in determining the interest rate. The rate of return increased by the occurrence of incidental cutting is greater than the non-risk rate of return. The considered risk is significantly affected by the decline in the yield value of forest land and the financial profitability of forest management [75]. Therefore, in our future work we would like to analyse the impact of different interest rates on the selection of optimal management variants.

Conclusions
The presented decision analysis tool based on the SIBYLA growth simulator and the optimisation tool seems to be useful for forest management planning of forest stands in central Europe that can help a decision maker examine a great number of possible harvest regeneration variants and their ability to fulfil a wide range of forest functions. The growth simulator contains tools which enable the application and simulation of a variety of harvest regeneration variants, and the interactive decision maps can be used by decision makers to examine the relationships between the required objectives and to find the optimal harvest-regeneration variant that fulfils the preferences of or compromises