Sustainability in Forest Management Revisited Using Multi-Criteria Decision-Making Techniques

: Since its origins, the idea of sustainability has always been linked to forest management. However, nowadays, sustainable forest management has usually been approached by deﬁning a set of criteria and indicators. This paper aims to address sustainability in forest management including a set of criteria encompassing the most common decisions: whether the stands are even or uneven-aged, and the optimal silviculture that should be applied in each stand. For this purpose, a lexicographic goal programming model with two priority levels has been deﬁned, into which six di ﬀ erent criteria are integrated. Each criterion corresponds to a particular pillar (economic, technical, or environmental). Furthermore, also incorporated into the model are the preferences of diverse stakeholders, both for the criteria considered in the analysis and for the most suitable silvicultural alternatives to be applied in each stand. This methodology has been applied to a case study in Spain, and the results show much more attractive solutions than the current forest management planning, allowing the obtainment of multi-aged systems that could be favourable for other ecosystem services.


Introduction
The concept of sustainability in a forestry perspective has a development history that dates over 300 years back in time [1], without the author himself having supplied any definition of this term [2]. In some respects, this idea at that time was the first conservation principle ever formulated [3]. In recent decades, it has spread in an exceptional manner to other areas [4]. To some authors, the sustainability concept is currently belonging to the field of social ethics [5]. Nowadays, it continues to be a forceful message with a wide representation in diverse disciplines [6]. In synthesis, von Carlowitz's premise promoted the concept of not exhausting resources in order to bequeath them to future generations, safeguarding finite natural resources for future generations [7].
In order to integrate sustainability into forest management, it has been customary to define a multidisciplinary ensemble of criteria and indicators in order to arrive at a consensus as to what sustainable forest management should be like [8], although sometimes it is not easy to transfer this idea to a strategic level [9]. It has been successfully applied on different spatial scales, using more or less aggregated information. Thus, there are indicators for sustainable forest management focused at a national or supranational level [10,11], or at a management unit level [12,13]. Further, in order to include this concept of sustainability in forest management, tools like certification systems aimed at ensuring sustainable management in forests have been developed. Thus, the systems most used nowadays simply require compliance with a series of indicators [14]. However, although it is commonly believed that a certification scheme implies a sustainability attribute, this direct relationship is not always clear [15].
based on a few basic principles in accordance with the type of forest management practiced at the time, which was centered on the "normal forest" or "normality" concept [60], established at the end of the 18th century by German foresters [61]. In fact, it was characterized by being a mono-objective management, focused only on timber production, usually in regular stands, in which one or a few main species prevailed [62], and to which a silviculture with a fixed rotation was applied [63]. In addition, a single decision-maker made the corresponding decisions with the aim of perpetuating the idea of a fully-regulated forest [57,64]. Finally, as a result, and unlike what happens nowadays [65], restrictions to final harvest were very lenient.
The concept of sustainability has been modified and updated due, among other reasons, to the changes in the role played by forests in society [66,67]. Besides that, the demand for goods and services coming from forest systems has been continuously increasing [68], meaning that, today, forest management is much more polyhedral (or multi-faceted, according to [69]) than it was three centuries ago, and includes other ecosystem services not contemplated in the original idea of sustainability, such as in the case of biodiversity [70]. Indeed, some authors have recently gone beyond the theory of von Carlowitz and defined forests as a "multifunctional bioeconomic system" [71]. Thus, numerous authors maintain that the new paradigms associated with sustainability should promote multifunctionality in forest management [72]. At the same time, the idea of sustainability in forest management has been linked to that of resilience [73]. In this regard, although the expression "sustainable" has for many years been on a par with forest management [74], nowadays it has been assumed that the latter should achieve certain goals with respect to economic, social and environmental pillars and, also, balancing present and future demands [75]. Some authors even propose the concept of "ecologically sustainable forest management" to give prominence to a type of forest management that perpetuates ecosystem integrity [76].
The way forest management has evolved has required the reformulation of the original sustainability concept [1] into the idea of a "regulated forest" (the organization of a forest property to provide a sustained yield of timber products forms [57]), by integrating the abovementioned pillars, but taking into account that, on some occasions, the trade-offs between different criteria are not fully known [77]. That is why sustainability should not be defined by only focusing on the conception of its being concentrated on "sustained yield", one of the objectives most sought after in forest management [78], or on other similar concepts [47]. Namely, this idea should be generalized starting from a case in which not all the forest should be associated with the idea of even-aged stands into which different silvicultural alternatives are established [63]; other ecosystem services are integrated, as well as timber production; and the preferences of the different stakeholders in the decision-making process are included. In effect, a multi-aged structure in forest stands is associated with the idea of long-term sustainability [79], incorporating multiple uses into different stands even with variable rotations [32]. Lastly, the inclusion of different silvicultural alternatives in the analysis is basically justified for two reasons. First, because different forest management strategies promote biodiversity [80], and, second, because silviculture as a science has to respond to the demands of society and not only concentrate on timber production [81].

Case Study
The case study corresponds to Valsaín forest, which is in the Central Mountain Range of Spain. The forest area covers about 7206 ha and mainly composed of pure Scots pine (Pinus sylvestris L.) stands, although there are other species less represented in the area such as Quercus pyrenaica L. and Ilex aquifolium L. The first management plan was developed in 1889, and the main objective was timber production. Since 1940, in order to achieve this objective, a uniform shelterwood system with a rotation of 120 years and a 20 years regeneration period has been proposed [52]. However, in recent years, this management has become more multifunctional [25], and the silvicultural treatment was turned into a shelterwood group system [52]. Nowadays, other aspects such as biodiversity conservation are of great importance. Thus, there are several endangered species in the area, like the black vulture (Aegypius monachus (Linnaeus, 1766)). The number of nests of this species is around 131, and harvest limitations have been introduced in order to protect them. In short, around each nest, a restricted area without management of 3.14 ha has been established as a buffer protection. Finally, there are several protection features in this forest, such as the declaration of the Sierra de Guadarrama National Park in 2013, where commercial cuttings have been forbidden, and which affects 3326 ha.

Silvicultural Alternatives
In order to deal with this complexity, seven silvicultural strategies have been proposed to maintain the adaptability of forests to different conditions [81], and to provide different ecosystem services. The set of possible treatments considers the nature of the stands forming the case study, and all of them have been discussed and approved by the forest manager. A priori, the seven treatments (see Figure 1B and Appendix A, Table A1) could be applied to every stand. was turned into a shelterwood group system [52]. Nowadays, other aspects such as biodiversity conservation are of great importance. Thus, there are several endangered species in the area, like the black vulture (Aegypius monachus (Linnaeus, 1766)). The number of nests of this species is around 131, and harvest limitations have been introduced in order to protect them. In short, around each nest, a restricted area without management of 3.14 ha has been established as a buffer protection. Finally, there are several protection features in this forest, such as the declaration of the Sierra de Guadarrama National Park in 2013, where commercial cuttings have been forbidden, and which affects 3326 ha.

Silvicultural Alternatives
In order to deal with this complexity, seven silvicultural strategies have been proposed to maintain the adaptability of forests to different conditions [81], and to provide different ecosystem services. The set of possible treatments considers the nature of the stands forming the case study, and all of them have been discussed and approved by the forest manager. A priori, the seven treatments (see Figure 1B and Appendix A, Table A1) could be applied to every stand.  The first alternative is based on a shelterwood group system. The treatment spans 40 years, cutting every ten years, and it would be the "business as usual" alternative (BAU). After the final cut, 2% of the stand volume remains as dispersed trees for at least one rotation. In order to model this alternative, we have considered the schedule and cutting intensity traditionally used in Valsaín forest [82]. BAU has been divided into two new alternatives: BAU 1 , embracing light thinning, and a more intensive option with strong thinning (BAU 2 ). The third alternative, called Group selection (GS), consists of 2 cuttings with medium patch sizes of 1 hectare. After final cutting, 5% of the stand volume remains as residual trees.
Another two silvicultural alternatives are based on the idea of Green Tree Retention (GTR). Two GTR alternatives with different retention percentages are proposed, before applying BAU 1 or BAU 2 . In short, either 15% of the volume is maintained in an aggregated way as mature forest patches (GTR 1 ), or 30% of it (GTR 2 ). These figures are similar to those proposed for this species in Spain [83]. These patches could simulate the protection buffers discussed above for the protection of the black vulture species.
For the area now included in the National Park, we have proposed an alternative (TNP) similar to BAU, but with a main difference: the final cut is excluded, due to prohibitions in all Spanish National Parks. Moreover, in this case, the rotation has been increased (up to 140 years) in relation to the other treatments, as a conservation measure [84]. Finally, a no management alternative (NOM) has been incorporated in order to facilitate the possibility of creating mature forest stands in the future, given that these stands are allowed to grow naturally, thus forming reserve areas.

Criteria Considered
As mentioned above, the criteria to be introduced should evaluate sustainability based on the multifunctionality intrinsic to forest systems. First, two technical criteria have been defined. The first of these corresponds to a synthetic index that encompasses the classic idea of a normal forest (NF). In short, this criterion embraces tree indicators: one associated with the notion of even-flow policy, another with the idea of regulation or regulated forest [57], and, finally, a third one related to a suitable ending forest inventory, which ensures the perpetuation of the forest. In order to aggregate the above indicators under the normal forest concept, the procedure used in [85] was followed. In addition, a criterion called Ideal Management (IM) has been introduced. It can be quantified by the amount of area (in hectares) to which the silviculture suitable for each stand, according to the criterion optimized in each case, is applied. Therefore, calculating this criterion consists of minimizing the deviations of the silviculture used, compared to that which would be preferred in terms of the characteristics of the stands and in accordance with the preferences reflected by a decision-maker.
On another side, two environmental criteria have been defined. The first one refers to the possibility that, at the end of the planning horizon, there is a multi-aged forest. To enable a comparison of the degree of irregularity between different prescriptions, the Gini index (G) was employed. The latter is frequently used to measure the degree of inequality in different variables, generally with economic data, although applications in many fields have been described [86]. It supplies values in the interval [0 1] and has been defined in such a way that values close to 0 imply the existence of multi-aged stands, whereas values close to 1 ensure an even-aged structure. The second environmental criterion considered was the carbon balance (C), defined as the difference between the estimated gains and losses of carbon due to the growth of trees and their removal (default method, following [87]) throughout the planning horizon. This criterion will be computed in physical units (tons of carbon).
With respect to the production criteria, the Net Present Value (NPV) and the volume (VOL) associated with the final cuttings (cubic meters) over the planning horizon have been considered. Regarding NPV (€), the incomes provided referred only to timber sales. We have considered a fixed maintenance cost (11.6 €/ha), and different harvesting costs in relation to the silvicultural treatment applied. Following data provided from the current forest management plan, the timber price considered was 28 €/m 3 . We have considered a real discount rate of 2%. Concerning the criteria considered, it is appropriate to clarify that VOL, NPV, IM and C are criteria belonging to the type "the more, the better", while NF and G are criteria belonging to the type "the less, the better", following the terms used in [35]. Figure 1B shows the six criteria considered with their acronyms, meanings, and the units in which they are expressed.

Multi-Criteria Model
When analyzing the solutions shown in the Pay-off matrix, if no solution is preferable, it is necessary to look for a satisficing compromise solution, making the harvest scheduling more sustainable. As we have said before, the technique chosen was Lexicographic Goal Programming (LGP), the second most used GP variant in forest management applications [88].
The first issue to decide in this type of model is how to select the number of priority levels that one wants to establish. Starting from the premise that it should be a limited number and that there should not be as many priorities as the number of criteria [89], in this case, two priority levels have been fixed. In the first one, three goals most directly related to silvicultural aspects have been included, while in the second, NF, IM and G are considered. Namely, given that there are no finite trade-offs among goals placed at the two priority levels [88], it is assumed to be a priority to achieve the goals at the first level before optimizing the criteria situated at the second priority level [90]. The target for each goal was established at 70% with respect to the ideal value obtained in the pay-off matrix. This percentage has been employed in the literature [91], and its use avoids the drawback of not obtaining optimal alternatives if the targets are closer to ideal values [92]. Finally, the complete model can be found in Appendix B, and for its resolution the software LINGO 13.0 [93] has been used.

Interactive Process
In this study, we interacted with two decision-makers (DM1 and DM2). First, as mentioned above (Equations (A2) and (A3)), it was considered appropriate to introduce preferential weights for each of the six goals contemplated in the analysis. For this purpose, the interaction was with the forest manager (DM1), and his opinion was asked on the six criteria using a survey based on pairwise comparisons, with Saaty's verbal scale [94] (see Appendix A, Table A3). This way of obtaining the decision-maker's preferences is recommended for integrating them into GP models [95]. These pairwise comparisons are typically used in the forest context [96].
On another side, and with the aim of obtaining the IM criterion, the engineer in charge of carrying out the last forest management planning (DM2) was contacted. Taking into account the silvicultural alternatives considered (see Appendix A, Table A1) the objective was to establish the suitability of each one for being applied to each stand. For this purpose, the 288 stands in the forest have been grouped into seven categories generated as a combination of a series of factors: production, protection, National Park, districts with nests or the presence of mixed stands. Once the categories were defined, in a first contact with DM2, we selected a set of 2-3 preferred alternatives for each stand. Later, after a second interaction with the same expert, some weights were obtained to weight the different silvicultural alternatives for their suitability, attending to the characteristics and nature of the stand in which they were found (see Appendix A, Table A2). To obtain the weights assigned to each silvicultural alternative in each stand, the same procedure as that employed with DM1 was used (see Appendix A, Table A3).

Results
Beginning with the pay-off matrix, Table 1 shows the results of the optimization of the six criteria separately. In Table 1, ideal values (principal diagonal of the matrix) are shown in bold type, while anti-ideal (or least desired ones) are underlined. The main diagonal includes the maximum values that each criterion can reach, known as the ideal point. Other parts of the matrix contain the anti-ideal points, which would be the worst results obtained for each of the criteria. Finally, to increase the informative character of this matrix, several additional rows were included by measuring the area occupied when each criterion is optimized by each of the silvicultural treatments described, as well as the average rotation. The results given in this table show the ideal and anti-ideal values reached by the criteria, according to the optimized criterion (per columns). First, it can be seen that some objectives do not reach ideal values (NF, IM, G). Then, the level of conflict existing between the different criteria, mainly between NF and IM, and between C and G, can be observed, as well as the fact that the results are fairly similar in relation to carbon balance, because the percentage for this criterion varies the least compared to the rest. In synthesis, it can be noted that there is no solution corresponding to the optimization of a single criteria that appears to be sufficiently attractive. In short, the results of this pay-off matrix justify the setting up of an MCDM model in order to achieve solutions that are more balanced from the point of view of sustainability.
However, before revealing the results of those models, the result of the first interaction with DM1 to obtain the preferential weights associated with the criteria is shown ( Table 2). Note that the sum of the weights for each priority level is equal to 1. Finally, Appendix A (Table A2) contains the weights obtained throughout the second interactive process with DM2 for the seven categories, in which the stands and the silvicultural alternatives associated with them have been grouped. Next, Table 3 displays the solutions obtained for the two priority levels proposed. It can be noted that the solutions at the second priority level (U 2 ) are the same as the values found in the preceding Sustainability 2019, 11, 3645 9 of 24 solution (U 1 ), whereas the values of the three criteria included in U 2 improve in comparison with those obtained in U 1 . Further, to facilitate comparisons, the same auxiliary rows have been included. In general, it can also be seen that the areas assigned to each silvicultural alternative are notably modified when moving from U 1 to U 2 .  In the auxiliary rows in Table 3, the area selected by the model for each silvicultural alternative is indicated. Also, this solution of the LGP model can be compared with the value proposed by the DM2 (the area pointed out as being the one preferred for applying each silviculture), information which is available on Appendix A (Table A2). In Figure 2 the difference (in surface) between both areas (applied or selected by the model and preferred) can be noted. However, it should be taken into account that the area preferred exceeds the value of the total forest area, because, for each stand, a set of 2-3 possible silvicultural treatments were proposed.
Next, Table 3 displays the solutions obtained for the two priority levels proposed. It can be noted that the solutions at the second priority level (U2) are the same as the values found in the preceding solution (U1), whereas the values of the three criteria included in U2 improve in comparison with those obtained in U1. Further, to facilitate comparisons, the same auxiliary rows have been included. In general, it can also be seen that the areas assigned to each silvicultural alternative are notably modified when moving from U1 to U2. In the auxiliary rows in Table 3, the area selected by the model for each silvicultural alternative is indicated. Also, this solution of the LGP model can be compared with the value proposed by the DM2 (the area pointed out as being the one preferred for applying each silviculture), information which is available on Appendix A (Table A2). In Figure 2 the difference (in surface) between both areas (applied or selected by the model and preferred) can be noted. However, it should be taken into account that the area preferred exceeds the value of the total forest area, because, for each stand, a set of 2-3 possible silvicultural treatments were proposed.

Discussion
The method here proposed for evaluating sustainability in forest management has the advantage of following the evolution of the criteria over the planning horizon [19], but it could also be adapted when including spatial constraints [97]. In that case, it could be suitable to design a model simultaneously integrating strategic and tactical planning [98]. Furthermore, the methodology proposed could be extended to incorporate the preferences of a group of stakeholders for addressing sustainability [99,100]. Some authors suggest compiling the forest managers' preferences on silvicultural aspects in a goal programming model [101], although the way these individual preferences would be aggregated is still an open question [102]. Finally, although the inclusion of dynamic aspects could trigger highly specific models, limiting their transference to others [103], we do not believe that this could happen in the proposal presented in this study due to its potential adaptation to other works, with other criteria and indicators.
The results presented in Table 3 allow us to observe the differences in relation to the values achieved for the indicators situated at the first priority level, at forest level. Thus, Figure 3 shows how the values of the criteria NF, IM and G obtained in the pay-off matrix are modified (Table 1), with respect to those obtained in the solution of the LGP model (Table 3).

Discussion
The method here proposed for evaluating sustainability in forest management has the advantage of following the evolution of the criteria over the planning horizon [19], but it could also be adapted when including spatial constraints [97]. In that case, it could be suitable to design a model simultaneously integrating strategic and tactical planning [98]. Furthermore, the methodology proposed could be extended to incorporate the preferences of a group of stakeholders for addressing sustainability [99,100]. Some authors suggest compiling the forest managers' preferences on silvicultural aspects in a goal programming model [101], although the way these individual preferences would be aggregated is still an open question [102]. Finally, although the inclusion of dynamic aspects could trigger highly specific models, limiting their transference to others [103], we do not believe that this could happen in the proposal presented in this study due to its potential adaptation to other works, with other criteria and indicators.
The results presented in Table 3 allow us to observe the differences in relation to the values achieved for the indicators situated at the first priority level, at forest level. Thus, Figure 3 shows how the values of the criteria NF, IM and G obtained in the pay-off matrix are modified (Table 1), with respect to those obtained in the solution of the LGP model (Table 3). As could be expected, in the comparison between models (Figure 3), the values of the criteria obtained in the LGP model are lower than the ideal ones reached in the pay-off matrix. These reductions indicate the opportunity cost associated with the simultaneous consideration of the six criteria defined in this model, instead of the values proposed in the individual optimization of each criterion. Thus, taking into account that the NF and G values (located at each end of the graphic) are adimensional and that they represent indicators of the less-the-better type, in the LGP solution their values are 2.7 and 1.6 times worse, respectively, compared to that obtained in the pay-off matrix. Finally, the IM criterion's value of 5358 ha, which represents the area managed by the most suitable treatment in each management unit (according to DM2), was reduced to 3750 ha in the LGP solution; i.e., about 30% of the area under ideal management (IM) is taken advantage of in LGP, using a different silvicultural treatment from the one that a priori was preferred for application in these areas.
From a sustainability perspective, it has been implicitly assumed up to now that the criteria selected functioned within weak sustainability parameters [31]. Namely, a certain degree of As could be expected, in the comparison between models (Figure 3), the values of the criteria obtained in the LGP model are lower than the ideal ones reached in the pay-off matrix. These reductions indicate the opportunity cost associated with the simultaneous consideration of the six criteria defined in this model, instead of the values proposed in the individual optimization of each criterion. Thus, taking into account that the NF and G values (located at each end of the graphic) are adimensional and that they represent indicators of the less-the-better type, in the LGP solution their values are 2.7 and 1.6 times worse, respectively, compared to that obtained in the pay-off matrix. Finally, the IM criterion's value of 5358 ha, which represents the area managed by the most suitable treatment in each management unit (according to DM2), was reduced to 3750 ha in the LGP solution; i.e., about 30% of the area under ideal management (IM) is taken advantage of in LGP, using a different silvicultural treatment from the one that a priori was preferred for application in these areas.
From a sustainability perspective, it has been implicitly assumed up to now that the criteria selected functioned within weak sustainability parameters [31]. Namely, a certain degree of compensation between the different indicators was presumed, although the LGP method implies that there is no finite trade-off between the different indicators of the different criteria considered [89]. However, in the literature, aggregated sustainability indexes based on multicriteria techniques including premises associated with strong sustainability can be found; i.e., it is found that there cannot be compensation between the criteria [104], nor can a certain degree of compensation between criteria even be chosen [105,106].
The methodology proposed here is open to incorporating other silvicultural alternatives. In our case study, GS has been proposed in order to favor the development of multi-aged stands [79], and GTR has been suggested in order to maintain multiple forest values [76], including biodiversity conservation [107], and, especially, to promote silvicultural treatments compatible with the protection of wildlife species [108]. The percentages of mature forest chosen here are those most frequently used in the literature [109,110]. Another silvicultural alternative (TNP) embraces the prohibition of final cuts (this forbiddance is mandatory in all Spanish National Parks). Maintaining this type of protected area is highly beneficial for certain species, and this justifies not employing a conventional silvicultural alternative [111].
On another side, and as mentioned above, other criteria could be integrated into the model, if considered to be appropriate. For example, we are aware that, unlike other studies [112,113], no criterion expressly linked to the social pillar of sustainability has been proposed. Although one study made in the same case study [114] did suggest two recreation indicators, it has not been possible to transfer this idea to our analysis. Along these lines, some studies point to an underrepresentation of social type indicators as being common [103].
Again, and aside from the significance of the Gini coefficient proposed here, which has been applied as an indicator of forest structural heterogeneity [115,116], we have attempted to define our own biodiversity criterion [117] that could be integrated into the LGP model. In fact, one of the objectives in the case study [118] is focussed on the conservation of the black vulture, an umbrella species [119], which is catalogued as "near threatened" at the global level [120]. For this reason, we pursue the possibility that, at the end of the planning horizon, there is a multi-aged forest. The inclusion of this criterion (Gini coefficient) under an environmental pillar is justified when the forest managers aim to promote the existence of uneven-aged stands as a biodiversity conservation measure [63,79].
However, it has not been possible to model the evolution of this species over 100 years, due to the lack of information available for predicting these values, so it has not been included in the analysis. Furthermore, it would have been necessary to determine the impact of the different silvicultural alternatives on the vulture populations, because it has been demonstrated that they have an influence on some requirements of the nesting of this species [121]. The lack of consideration of this criterion illustrates the difficulty in finding criteria compatible with the information assembly over the planning horizon, in spite of the availability of suitable methodologies for their possible integration into the model proposed. However, some studies have been developed under the goal programming technique successfully integrating aspects related to biodiversity [122,123]. Lastly, the methodology proposed here can also be extended, incorporating other silvicultural alternatives centred on other ecosystem services different from timber production [124].
Furthermore, the results of the LGP model permit one to obtain an overall value of sustainability derived from the resolution of that model, making it act as an aggregated index of the sustainability of forest management, assimilating the procedure reported in [125]. Indeed, if the objective is to compare sustainability between similar forests, or between parts of the same one, the solution of the LGP model could allow them to be ranked in terms of greater or lesser aggregate achievement. In [12,13], some examples of this extension are shown. Finally, the opportunity cost in terms of sustainability (gains or losses of sustainability) could also be calculated, by obligating a single silvicultural alternative to be used in some part of the case study. This circumstance can be of interest when, as can be seen in this case, the solutions notably vary with respect to the area that each silvicultural treatment has to cover, and this also happens in other studies [126].
In addition to incorporating new sustainability criteria, one aspect of the methodology proposed here that should be highlighted lies in the combination of flexibility and complexity that can be incorporated, for example, in terms of the silvicultural alternatives considered. The advantage of this flexibility is indispensable for tackling typical contemporary challenges in forest management [127,128]. Since there is no single specific definition of the idea of complexity in silviculture [129], some authors affirm that the impact of certain forest practices on forest system dynamics can compromise their sustainability [81]. Thus, the consideration of diverse silvicultural alternatives for approaching different criteria would lead to what some authors already call "complex adaptive systems" [130]. Lastly, this methodology can be applied at different levels (spatial or temporal) at which it is desired to carry out the planning and silviculture [128], contributing to the flexibility of the proposed methodology.

Conclusions
In this study, a flexible methodology was proposed to address sustainability in forest management. It includes the possibility of incorporating different criteria, or even synthetic indicators, multiple silvicultural alternatives for each stand and the preferential weights of different stakeholders. The solutions of these multi-criteria models allow the decision-maker to select, in terms of sustainability, the optimal silvicultural alternative for each stand along the planning horizon, the degree of multi-aged forest proposed and the fulfilment of the normal forest ideal.
The results obtained in the case study show different solutions from the ideal values shown in the pay-off matrix, with those supplied by the LGP model being much more attractive. Moreover, it can be seen how the areas devoted to each silvicultural alternative are notably modified in each of the solutions obtained. Finally, the flexibility of the model leads us to believe that this methodology could easily be extended to other case studies and planning levels, both temporal and spatial. Thanks are also given to Patricia Riquelme (Sierra de Guadarrama National Park) and Miguel Cabrera, and to Diana Badder for editing the English. Additionally, this manuscript has been proofread by PRS and MDPI. Thanks are given to the Editor and the reviewers for their helpful comments and criticisms, which have greatly improved the presentation and accuracy of the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
being subject to: Starting from v priority levels (Equation (A1)), for each of the priority levels cited, some functions are defined (Equations (A2) and (A3)), composed of the unwanted deviation variables corresponding to the goals (Equation (A4)) situated at those priority levels. t q is the target level for the qth goal, n q and p q are the negative and positive deviations from the target value of the qth goal. α q = w q R q is used if n q is unwanted, otherwise is used α q = 0, β q = w q R q if it is unwanted, otherwise β q = 0. The parameters w q and R q are the weights reflecting preferential and normalising purposes attached to the achievement of the qth goal.

Appendix B.2. Strategic Planning
Following Figure 1A, a strategic forest planning model was designed with a planning horizon of 100 years, divided into 10-year periods, which is a typical timeframe in models with slow-growing species, especially with Scots pine [118]. In short, a timber harvest scheduling model following Model I methodology [121] was built. Initially, Model I requires the definition of the decision variables or prescriptions, which need to encompass all possible management regimes to be considered, taking into account the rotation length and the silvicultural alternatives put into the model, which will remain unaltered throughout the planning horizon [33]. This strategic planning model was fed with the information available in the last forest management plan [82]. In addition, the abovementioned silvicultural alternatives were included in this general model, and we fixed an interval for rotations of between 120 and 180 years. The lowest limit corresponds to the current rotation length for Scots pine in the Sierra de Guadarrama [52], while the upper limit was established at 180 years in order to provide mature forest habitats, which are essential for the survival of several endangered species [84,131]. On the other hand, it is convenient to clarify that, in addition to the endogenous restrictions two exogenous restrictions have been introduced: a minimum harvest area for each silvicultural treatment, and a minimum commercial volume from final cuttings (TNP not included) per period. Finally, considering the rotation length and the 288 stands defined in Valsaín Forest, a set of 11,148 prescriptions was generated.

Constants
T is the length of the planning horizon. t is the time period. L = T/t is the number of periods into which the planning horizon is divided. ts is the time span that defines the final age class. S = T/ts is the number of final age classes, taking into account the range of years in which regeneration can occur.
B is the total forest area. B i is the area in stand i. b BAU 1 is the minimum harvest area when x BAU 1 prescriptions are selected. b BAU 2 is the minimum harvest area when x BAU 2 prescriptions are selected. b GS is the minimum harvest area when x GS prescriptions are selected. b GTR 1 is the minimum harvest area when x GTR 1 prescriptions are selected. b GTR 2 is the minimum harvest area when x GTR 2 prescriptions are selected. H min is the minimum volume harvested per period. mAge and MAge are the minimum and maximum harvest ages during the planning horizon, respectively.
V ijl is the volume harvested per hectare in stand i, prescription j at period l. F i z is the initial forest inventory on z site index. BVN is the Black Vulture Nest protection area. CF 0 is the initial investment and CF t is the cash-flow per year. r is the discount rate. γ is the parameter that expresses the density of wood in tons per m 3 (0.5) and the percentage of carbon content (50%) existing in dry biomass, whose value is obtained through γ = 0.5· 0.5 = 0.25 tC/m 3 .

Index Sets
i is the number of stands. j is the number of prescriptions, defining a complete treatment schedule for each stand, following the Model I structure (Johnson and Scheurman, 1977).
k is the number of initial age classes. l is the index of periods. z is the site index. v is the number of priority levels.

Variables
x ij is the area harvested at prescription j in stand i.
x BAU 1 is the area harvested by BAU 1 at prescription j in stand i.
x BAU 2 is the area harvested by BAU 2 at prescription j in stand i.
x GS is the area harvested by GS at prescription j in stand i.
x GTR 1 is the area harvested by GTR 1 at prescription j in stand i.
x GTR 2 is the area harvested by GTR 2 at prescription j in stand i.
x TNP is the area harvested by TNP at prescription j in stand i.
x NOM is the area selected for NOM at prescription j in stand i. ux i is a binary variable to force decision variable x i to take either a zero value or a value greater than, or equal to, the minimum harvest area designated by parameters b BAU 1 , b BAU 2 , b GS , b GTR 1 and b GTR 2 .
A s is the area belonging to s final age class s at the ending period.
F f z is the ending forest inventory from site index z. SV ij is the total standing volume at prescription j in stand i. P is the cumulative percentage of the number of trees per hectare for each stand i. Q is the cumulative percentage of volume. CC ij is the carbon capture at prescription j in stand i. EC ij are the carbon emissions produced by cuttings at prescription j in stand i. R q is the normalized vector (as the difference between ideal and anti-ideal values) for each criterion q.

Appendix B.4. Definition of LGP Model
Criteria NF is the normal forest as an aggregated normalized function of the even flow harvest volume per period (H l ), the forest ending inventory for each site class (F z ) and the area control for each age class (A s ).
IM is the ideal management to be applied in each stand i. G is the gini index as an inequaluty measure. VOL is the volume harvested at the end of the planning horizon. NPV is the Net Present Value at the end of the planning horizon. C is the carbon balance at the end of the planning horizon.

Achievement function
R VOL + w NPV · n NPV R NPV + w C · n C R C U 3 = w NF · (n NF ) The following constraints were introduced into the LGP model. First, endogenous ones were introduced to ensure that the area chosen by the model cannot exceed the available area in each stand. Second, a minimum harvest area was proposed for each treatment, since, as the treatments propose cutting sequences at several periods (2-4), a minimum area necessary for the intervention (5 ha for BAU 1 and BAU 2 , and 10 ha for GTR 1 and GTR 2 ) for each of them, was considered. In the case of GS, this minimum area refers to the size of the gaps in which the thinnings are applied. Finally, for TPN and NOM, no amount of minimum area was imposed for their application. The last constraint ensures a minimum harvested volume per period, according to the current forest management plan and the timber harvested in past years.