Multi-criteria analysis to compare multiple risks associated with management alternatives in planted forests

Aim of study: Adaptation of silviculture in planted forest may help to mitigate damage due to biotic and abiotic hazards. However, compromises have to be found because it is not possible to minimize the risk from all hazards through application of a single forest management approach. The objective of this study was to improve a multi-criteria risk analysis (MCRA) method that makes it possible to rank forest management alternatives (FMAs) according to multiple risks. Material and Methods: We defined eight FMAs for maritime pine forests in France, Spain and Portugal. We used as the definition of risk the combination of hazard, susceptibility and exposure. Hazard level was estimated using archive data on occurrence and severity of damaging agents over the last few decades. Forest susceptibility to hazards was evaluated by experts who scored the effect on stand resistance of eleven silvicultural operations characterizing each FMA. Exposure was estimated as value at stake, which combined forest standing volume, simulated with forest growth models, and wood prices.Main Results: Using the PROMETHEE algorithm, we found that the overall ranking of FMAs was consistent across all countries, with short rotation plantations to produce pulpwood or energy wood were the least at risk. The ranking was mainly driven by forest values at stake. We found that by improving the accuracy of forest values exposed to damage, based on growth models and representative wood prices, the MCRA outcomes were more useful and realistic.Research highlights: Our methodology provides a relevant framework to design FMAs that would minimize risks while maintaining income.Keywords: Pinus pinaster; vulnerability; hazards; growth modelling; expert assessment; wood price; southwestern Europe.


Introduction
Planted forests represent 7% of the global forest area and their surface is increasing worldwide (Payn et al., 2015). They are likely to play a significant role in the near future to meet the increasing demand for wood products and therefore decrease the logging pressure on natural forests (ICPF, 2013;Pirard et al., 2016). Meanwhile, planted forests are also expected to provide many other ecosystem services such as regulating services (e.g. carbon sequestration, soil conservation and watershed resources) but also biodiversity conservation at the landscape scale (Brockerhoff et al., 2013;Thompson et al., 2014). However, forests are increasingly affected by biotic and abiotic hazards (Lindner et al., 2010;Seidl et al., 2011a;Albert et al., 2015). In many parts of the world, pest and diseases are expanding their range and increasingly causing damage thus compromising the delivery of forest ecosystem services (Schelhaas et al., 2003;Boyd et al., 2013;Seidl et al., 2014). In recent decades, forest disturbance regimes have intensified markedly in Europe, resulting in much higher levels of damage from hazards such as windstorms, insect outbreaks and wildfires (Lindner et al., 2010;van Lierop et al., 2015). Moreover, climatic changes were identified as key drivers of this increase (Seidl et al., 2011c;Boyd et al., 2013) mainly through change in hazard occurrence, duration, severity and intensity (e.g. extreme climatic events such as drought, Allen et al., 2010) and forest susceptibility (e.g. water stress reducing tree resistance to pest, Jactel et al., 2012b). However, changes in forest management can also increase the level of damage (Jactel et al., 2009). For example, increasing both the age of the final cut and the standing volume increase storm vulnerability (Gardiner & Birot, 2013). Besides, reduced tree diversity can result in higher susceptibility to biotic and abiotic disturbances (Jactel et al., 2017). Because 99.9% of planted forests are monoculture plantations (Nichols et al., 2006), producing ca. 50% of industrial round wood (Payn et al., 2015), it is essential that new methods for evaluating the risk to planted forests are developed and implemented.
Risk depends on the combination of three components (Kron, 2005;Jactel et al., 2009): 1) the occurrence and severity of the hazard (e.g. wind, fire, pathogen, pest) which is the cause of damage, 2) the susceptibility of the system (e.g. trees and forest stand) to the hazard which determines the amount of damage when it occurs and, 3) the value exposed to damage, i.e. the socio-economic impact that depends on the value at stake. There is a multitude of biotic (e.g. insect pests, fungal diseases, pathogenic nematodes, mammal grazers) and abiotic (e.g. drought, frost, storm, fire) hazards in forests. In addition, they do not operate independently but often interact or generate cascade effects, leading to synergistic adverse effects on tree growth or survival. For example, drought often increases tree sensitivity to infection by pest and diseases (Jactel et al., 2012b) and windstorm triggers bark beetle outbreaks (Temperli et al., 2013). There is thus a need to consider forest risks in a holistic perspective, i.e. to operate a paradigmatic shift from individual to multiple risk analyses.
While several reviews demonstrated the possibility to reduce stand susceptibility through changes in silvicultural treatments (e.g. Fettig et al., 2007;Jactel et al., 2009;Klapwijk et al., 2016;Jactel et al., 2017), they also pointed out that a given particular operation (i.e. a cutting or any tending operation such as pruning or bush cleaning…) can have multiple, and sometimes, contradictory effects on stand susceptibility to different damaging agents. For example, thinning can improve individual tree vigor, through reduced competition, thus increasing individual tree resistance to bark beetles (Fettig et al., 2007), whereas thinned stands can be more heavily attacked by tree defoliators (e.g. pine processionary moth, Régolini et al., 2014). It follows that no single forest management alternative (FMA) can simultaneously minimize all risks of damage and that compromises have to be made.
Several empirical or process-based models have been developed to link forest stand management with susceptibility to specific damaging agents and infer on associated risks but only a few can deal with two combined hazards, e.g. wind and bark beetles (Temperli et al., 2013), drought and bark beetles (Lexer et al., 2002) or fire and bark beetles (Simard et al., 2011). However, none, to our knowledge, are generic enough to incorporate multiple hazards in many different forest conditions. There is thus a need for a flexible and robust methodological framework enabling to rank different forest management alternatives according to their relative exposure to multiple hazards.
Multi-criteria decision analysis (MCDA) methods provide a suitable framework and effective techniques for finding compromise solutions (Ishizaka & Nemery, 2013). They are flexible with input data (able to combine qualitative and quantitative information), propose mathematical algorithms to solve complex questions of data aggregation while being user friendly, thus enabling decisions makers to participate in the process (Ishizaka & Nemery, 2013). One of the most widely used MCDA method is PROMETHEE (Preference Ranking Organisation Method for Enrichment Evaluation). It is an outranking method, typical of the European MCDA school (Brans et al., 1986). It has been successfully used to address decision making problems for risk management (e.g. Linkov et al., 2006;Ruzante et al., 2010;Su & Tung, 2013;Liu et al., 2016;Doyi et al., 2018) and also forest planning (e.g. Huth et al., 2005;Ananda & Herath, 2009;Ghaffariyan et al., 2013;De Barros et al., 2014;Segura et al., 2014;Schuler et al., 2017) or both (Seidl et al., 2011b). In 2012, the method was adapted to forest risk assessment and the term MCRA was coined for "Multi-criteria Risk Analysis", which enables ranking forest management 3 Multi-criteria risk analysis in European maritime pine planted forests alternatives according to their vulnerability to multiple hazards (Jactel et al., 2012a).
We focused on forest risks for maritime pine (Pinus pinaster Ait), a native and major timber tree species of southern Europe, to exemplify our approach (Abad Viñas et al., 2016). The study area encompassed three regions: Aquitaine (France), North Portugal and Galicia (Spain), where maritime pine is the main tree species in productive forests and accounts for ca. 47, 45 and 44% of the forested area, respectively (Sanz et al., 2006;AFN, 2010;IFN, 2010). Maritime pine is also one of the main harvested wood material in the three studied regions (MAGRAMA, 2011;Dias & Arroja, 2012;DRAAF Aquitaine, 2015). Maritime pine forests thus play a major role in the regional economies. Recently, large disturbances occurred in the three regions: in 2009, in Aquitaine, the storm Klaus blew down approximately 43 million m 3 of wood, which represented an economic damage of approx. 2 billion Euros (Hanewinkel & Peyron, 2013); wildfires that occurred in Portugal during 2003 (312,430 ha damaged area) and 2005 (274,783 ha) were responsible for economic losses of over 3 billion Euros (ICNF, 2013); in 1999 the pine wood nematode was discovered in Portugal and five million trees had been felled by 2008 (Rodrigues, 2008). In 2006, forest fires affected an area in Galicia of 86,000 ha (around 60% covered by trees), being responsible for economic losses of 248-336 million Euros (Loureiro & Barrio, 2009). In 2017, major fires burnt another 46,000 ha in Portugal (with 65 fatalities) (Nunes et al., 2019). These large-scale events raised a lot of concern among foresters about the sustainability of wood production in cultivated maritime pine forests. Such a context therefore offered a good opportunity to engage with forest stakeholders in comparing forest management alternatives (FMAs) regarding multiple forest risks.
The main objective of this paper is to move a step forward and present an enhanced version of MCRA that addresses previous limitations. In particular, the study aimed at improving two key elements in the multi-criteria risk analysis framework, i.e. at proposing and testing: i) more accurate estimates of hazards levels, combining occurrence and severity and based on archive data instead of expert opinion, ii) more meaningful values at stake, based on economic valuation of forests using forest growth models and real wood prices, again instead of a simple score provided by experts.

Methods
We first present the Forest Management Alternatives (FMAs) designed by forest experts according to local silvicultural practices in the three selected regions. Then we assess the three components of the risk used in this paper as "hazard level", "forest susceptibility" and "wood value at stake" (Fig. 1). They were combined to rank the forest management alternatives (FMAs) according to associated multiple forest risks in a given region, using the outranking method PROMETHEE (Mareschal, 2013). PROMETHEE was developed by Brans (1982), further extended by Brans & Vincke (1985) and Brans & Mareschal (1994).

Forest management alternatives (FMAs)
Eleven silvicultural treatments were defined: site preparation, fertilisation, regeneration type, tree genetic material, stand structure, stand composition, cleaning of understorey vegetation, stand clearing, thinning, pruning and final harvesting (Table 1). For each treatment, several options were considered. The combination of these options enabled the construction of eight contrasting FMAs exploring a wide range of possibilities (e.g. with large variation in rotation length, site management techniques and tree genetic materials) and with diverse objectives in terms of quantity of timber and forest biomass production. The FMAs are explained in the next paragraph and described in Table 1 under the following simplified names: M1-High quality, M2-Standard classic, M3-Low investment, M4-Short-term with subsidies, M5-Low density without thinning, M6-Half-dedicated to biomass, M7-Biomass, M8-No management.
Some FMAs are currently applied in the case studies regions: for example, the "Standard classic" M2 is very common in Aquitaine and in Galicia and aims at producing timber in the medium-term. "High quality" M1 is representative of the objective of producing high quality timber on the long-term. The "Low investment" M3 is designed to consider stands in which forest owners invest as little money as possible before getting a return on investment with the first thinning. The "Short-term with subsidies" M4 is designed considering subsidies for planting, leading to active silviculture until the first thinning and then applying no other treatment until the final cut at 25 years, which is close to the official Galician Regional Government models PP2 (DOG Diario Oficial de Galicia, 2014). The "Low density without thinning" M5 is aimed at limiting thinnings and associated costs by planting at final density and harvest earlier. Others FMAs are more innovative, like the "Biomass" M7 aiming at producing energy wood and pulp-wood rapidly and the "Half-dedicated to biomass" M6 in which half of the stand is early thinned for energy wood and the other half is kept longer for timber production. They are considered possible options to restore quickly wood resources for the local forest-based sector in Aquitaine (Lesgourgues & Drouineau, 2009) and Galicia following the Klaus storm in 2009. Finally, "No management" M8 is a FMA with natural regeneration and no further management. It has no wood Forest Systems August 2020 • Volume 29 • Issue 2 • e004 production objective. However, harvesting of high-grade timber may occur in small private forest properties.

Hazard level assessment
We let the list of hazards vary with regions for the sake of realism but we set as constraint to comprise both abiotic and biotic hazards. Therefore, at least two abiotic and two biotic hazards were chosen in each region according to their local relevance based on local occurrence data for the last 10 to 40 years ( The invasive pine wood nematode (PWN, Bursaphelenchus xylophilus Steiner and Buhrern) and drought, even though locally relevant in the studied area (southwestern Europe) were not included in the analysis. PWN, a quarantine organism within the European Union, was not retained because it requires specific phytosanitary measures (mainly sanitary cuts) that were not included in the FMAs analysed here. Drought symptoms have been shown lagging one or several years behind the climatic event (Bréda & Badeau, 2008) and cannot be easily disentangled from damage caused by biotic hazards, making drought severity difficult to estimate properly.
After, hazard selection, hazard level was defined as a combination of mean occurrence and severity; it was defined for each region based on archive data (Fig. 1). We estimated hazard occurrence as the mean percent of forest area or volume damaged each year in the region (P1 Hi ). It was averaged across the longest possible period with reliable data, which was often of several decades. However, hazard occurrence is not enough to account for the relevance of damage. For example in Aquitaine, 16% of the surface of maritime pine forest have been impacted per year by the pine processionary moth since 1981 (Data from the French Forest Health Department) whereas only 1.5% of the volume has been impacted per year by storms since 1976 (communication from CNPF, Figure 1. Conceptual flow diagram of multi-criteria risk analysis (MCRA) methodology for a given regional case study.
Multi-criteria risk analysis in European maritime pine planted forests 2014). On the other hand, damage caused by the Klaus storm (2009) had much larger economic consequences on the forestry sector than chronic defoliation by pine processionary moth. To account for this issue, we estimated hazard severity as the percent of losses induces per area or volume affected (P2 Hi ) (see more details in Table 2). We used various estimates depending on hazards: losses due to occasional damage such as those produced by storms, fires or bark beetles are usually well assessed and a direct loss on wood price is generally observed, e.g. 86.5% for storms and fire on average ( Table  2). On the contrary, for some other hazards the consequences are often not discernible on the short-term, but become visible after a certain amount of time or a certain threshold of damage. They operate in a similar way as the creeping environmental problems defined by Glantz (1994) as "incremental, low-grade, cumulative and longterm biospheric changes induced by humans". In this category, we included herbivore damage (game) (Landmann et al., 2015) and insect defoliators like the pine processionary moth. It was more difficult to estimate hazard severity for these problems and we mainly relied on relative growth loss estimates (Ballon & Hamard, 2003;Gatto et al., 2009). To estimate P1 Hi and P2 Hi , we used scientific or grey literature, including regional reports or studies in local languages (Appendix A1 [suppl.]). We multiplied the two values to provide an estimate of absolute hazard level in each region. It was further divided by the sum of all hazards levels in a given region in order to produce 'relative level' per hazard and region (eq. 1). These values were used to weight criteria in the MCRA (Fig. 1).
Where W Hi is the relative weight of hazard Hi; P1 Hi is the mean percent of forest surface or volume damaged by hazard Hi each year, P2 Hi is the loss in wood biomass or price in a given region due to hazard Hi (n=6 for Aquitaine, n=6 for Portugal and n=4 for Galicia).
For example (see in Table 2), the average volume affected by storms per year was 1.5% in Aquitaine (P1 storm =0.015). Then, the loss due to storms was estimated to be 86.5% of wood price (P2 storm = 0.865), using the example of what happened after the Klaus storm (Lesgourgues & Chantre, 2009). Therefore storm level in Aquitaine was estimated as P1 storm ×P2 storm =0.0130. The sum of P1 Hi ×P2 Hi for all hazards in Aquitaine was 0.0746. The relative weight of storm in Aquitaine was thus 0.013/0.0746=0.174 (17.4%). The calculation of hazard relative weights in Aquitaine is fully documented in Tables 2 and 3 and details for Portugal and Galicia are provided in Table 3 and Appendix A1 [suppl.].
Last, to account for the fact that the number of hazards varied between regions and might bias the comparison of MCRA outcomes across regions, we re-did the analyses with a selected set of four hazards by region, each of them The defoliator is the pine processionary moth Thaumetopoea pityocampa, the root rot fungus is Heterobasidion annosum. "Bark beetle" refers to Ips sexdentatus. P1 H is hazard occurrence estimated as mean area or volume affected by a hazard each year. P2 H is % of losses induced per area or volume affected. The W H is the product of P1 H by P2 H divided by the sum of all hazard levels in a given region giving the relative losses per year induced by a given hazard.

7
Multi-criteria risk analysis in European maritime pine planted forests being present in at least two regions (Table 4). This change had no major impact on the conclusions.

Forest susceptibility assessment
Following the method developed by Jactel et al. (2012a), we used expert panels to evaluate maritime pine stand susceptibility to multiple hazards.
Panels: In each region, a panel of experts was gathered during a one-day workshop. It was highly recommended to have several experts per hazard but one expert could be consulted about different hazards. Experts were highly qualified researchers in forest ecology or forest managers with expertise in targeted hazards and maritime pine silviculture, belonging to the leading universities, research institutes, advisory companies or risk management authorities of each region. The number of experts per panel was between 10 and 15 people, because bigger groups are difficult to moderate and the number of available experts in each region was anyway limited. The number of experts participating in the panels and the expertise covered are detailed in Appendix A2 [suppl.].
Scoring method: Regional experts had to reach a consensus before giving a score from 0 to 1 for the effect of each silvicultural treatment (Table 1) on stand susceptibility to individual hazards. For each silvicultural treatment, they identified the "reference" option triggering the lowest amount of damage for a given hazard and gave it a score of zero. Then the experts gave a score of 0.25, 0.50, 0.75, or 1.0 respectively, to any other option that would increase a little amount, a moderate amount, a large amount or an extreme amount the damage due to a given hazard as compared to that of the reference option. We note that this was a semi-quantitative approach by which experts produced scores according to a relative scale from non-susceptible (0) to most susceptible (1) to a given hazard. Still, in multi-criteria decision analysis it is not the absolute numbers (net flows) that matters to qualify a given option (here a FMA) but its ranking, i.e. relative position with  respect to the others (i.e. whether it should be preferred or not to the others). At the end, for each of the eight FMAs, scores were summed across the eleven silvicultural treatments (Tk) that characterized a FMA (Mj), for each individual hazard considered (Hi), in order to get a general score of susceptibility (S) (Fig. 1)

Assessment of wood value at stake
The score of forest value at stake, that was usually provided by experts in Jactel et al. (2012a), was here replaced by a quantitative estimate of the wood standing value of a virtual forest that is composed of N stands of one hectare corresponding to each year of a rotation cycle of N years. The wood standing value was estimated for each stand (i.e., the entire virtual forest) and then divided by the number of stands (equal to the number of years in the rotation cycle). This provided the value at stake of a one hectare of this "virtual forest" where all ages and corresponding wood values would be represented, irrespective of the rotation length (age of trees at final harvesting time). This method of calculation enabled testing for the impact of all hazards on all possible ages of a stand managed with a given FMA. It should be noted that age effects were taken into account when assessing the susceptibility of stands to different hazards.
To evaluate the forest growing stock, each FMA, except M8, was simulated using a maritime pine growth model adapted to regional conditions. The models used were PBRAVO (Páscoa, 1987) in Portugal, GesMO (González González et al., 2012) in Galicia and the Lemoine model in Aquitaine (Lemoine, 1995;Meredieu & Labbé, 2006). Each model had its own equations (see the FORMODELS database http://www.iefc.net/formodels_database_forest_ modeles_liste/ for details). To start a simulation, all the models need the values of initial tree density, stand basal area at an initial age of 10 -12 years old and a site index (dominant height at a reference age; one average value per region). Fertilisation and improved genetic materials were taken into account by using higher site indices. The estimation of wood volume produced in M8 (No management) was based on data from the study of Sánchez-Orois & Rodríguez-Soalleiro (2002) for Galicia and from selected national forest inventory plots of maritime pine stands with irregular structure in Aquitaine (IGN, 2012). As we missed similar data for Portugal, we used the Aquitaine data to simulate M8 in Portugal after checking the plausibility of productivity figures.
Outputs of growth models were stand characteristics at every age including tree density, basal area, mean height, mean diameter, and total volume. Therefore growing stock in cubic meters was available for all stands composing the "virtual forest" for each FMA. Then current wood prices (S4) were used to convert the cubic meters into Euros and these values were averaged across stands to provide a wood standing value (E M ) for the "virtual forest" corresponding to each tested FMA. Wood prices were extracted from regional data based on wood sales in Aquitaine (2014), information from the Centro PINUS in Portugal (2016) and information from the "Asociación Forestal de Galicia" and auctions of the Galician Regional Forest Administration for Galicia (2011)(2012)(2013)(2014)(2015)(2016).
We defined vulnerability as the product of stand susceptibility to a particular hazard by its exposed value (Jactel et al., 2012a) (Fig. 1). To estimate the vulnerability of a FMA to a particular hazard (V M,H ), we multiplied the wood standing value in Euros (E M ) by the score of forest susceptibility (S M,H ) (eq. 3) = × (eq. 3)

Data integration into MCRA tool and sensitivity tests
We used Visual PROMETHEE ©, a multi-criteria decision analysis software (Mareschal, 2013) to process the data. The software can deal with quantitative and qualitative input data and implements the PROME-THEE algorithm and GAIA multi-criteria decision tools (Mareschal, 2013).
We ranked the FMAs (Mj) according to multiple criteria, using as criterion the maritime pine forest vulnerability to a given hazard (V Mj,Hi ). We used the relative level of hazard to weight the criteria (W Hi ). One multi-criteria decision matrix was thus built per region and the 3 matrices were incorporated in Visual PROMETHEE © (Appendix A5 [suppl.]).
We performed PROMETHEE II analyses (Brans et al., 1984;Brans et al., 1986) to make a complete ranking of FMAs according to multiple risks (see Jactel et al., 2012a for details). As main output, we used the net flow  , which is the difference between outgoing  + and incoming  -flows. The outgoing flow + estimates how far a given FMA outranks other FMAs and the incoming flow -, estimates how far it is outranked by other FMAs.
The calculation of  needs the choice of a preference function, which translates the difference in criterion value between two FMAs to a preference value ranging from 0 to 1, thus enabling the combination of several criteria.
Multi-criteria risk analysis in European maritime pine planted forests Each preference function is associated with a preference threshold (which is the smallest difference that is considered sufficient to generate a true preference). As in the Jactel et al. (2012a) study, we chose the V-shape preference function for all criteria because preference was likely to decrease proportionally to both susceptibility and value at stake. The distribution of vulnerability per FMA and per hazard, which led us to choose the V-shape model, is available in Appendix A6 [suppl.]. For each hazard and region, we sorted the vulnerability values in ascending order, calculated increment in vulnerability values between two successive FMAs and averaged increment values. Then we averaged the mean increment values across all hazards by region, which we used as preference threshold. The preference thresholds were 4238 for Aquitaine, 2987 for Portugal and 1419 for Galicia.
Because the main objective of the present study was to identify which FMAs were more or less at risk compared to others, we decided to set the decision rule to "minimize" all criteria, i.e. to minimize the product of exposure, vulnerability and hazard for any damaging agent considered. The classification would then rank first the FMAs least at risk. We also used the multiple scenarios tool in Visual PROMETHEE© to combine the three multi-criteria regional decision matrices and then provide an overall complete ranking of the FMAs irrespective of the regional case studies.
We then performed several types of sensitivity tests to assess the relative importance of the three risk components, hazard level, forest susceptibility and wood value at stake, in the ranking of FMAs. To rank FMAs irrespective of hazard agents, we set all weights to be equal. Then to rank FMAs irrespective of wood values at stake, we used only the susceptibility scores (S) instead of vulne-rability (V) in the multi-criteria decision matrix. We also combined the two by setting all weights equal and ranking FMAs only according to their susceptibility. Finally, we tested the impact of the number of hazards taken into account by doing the analysis only with the four hazards selected in Table 4.

Wood values at stake
The forest standing values (E) greatly varied between FMAs and to a lesser extent between regions for a given alternative (Fig. 2). Standing values increased with the duration of forest rotation (e.g. from M4-M7 to M1-M3) with the exception of M8 where the low quality of forest products explained a low value at stake. On average, standing values were generally higher in Aquitaine than in the other two regions due to higher wood prices in this region (Appendix A4 [suppl.]).

Complete ranking of FMAs according to associated risks
The ranking of FMAs was remarkably consistent across regions (Fig. 3) except for the "No management" (M8) and the "Low density without thinning" (M5) alternatives. The ranking of these two FMAs greatly varied between regions, being better ranked in Aquitaine than in Galicia. This is likely due to difference in hazard relative weights and susceptibility scores, as values at stake were quite similar (Fig. 2). The most preferred FMA was

Forest Systems
August 2020 • Volume 29 • Issue 2 • e004 always the "Biomass" (M7), followed by the "Short rotation" (M4). These two FMAs had in common short rotations (from 9 to 25 years) and produce small diameter timber or wood biomass. The least preferred FMAs (at higher risk) were the "High quality" (M1) and "Low investment" (M3), which were characterized by long rotations (ca. 55-60 years) and high expectation in terms of forest products (large stem wood) (Fig. 3 and Fig. 4a). "Standard" (M2) and "Half dedicated" (M6) were intermediate, being of medium rotation length (35-40 years) with an objective of stem wood production. Using only four common hazards agents (Fig. 3) did not change the ranking of FMAs, except for M8 in Aquitaine and M8 and M3 in Portugal, which were ranked lower than in the analysis with all hazards.

Ranking of FMAs to assess the relative importance of the three risk components
Setting all weights to equal, i.e. considering all hazards to be of equal damaging effect within a given region, hardly changed the ranking of FMAs (Fig. 4b).
By contrast, considering only the susceptibility of FMAs (i.e. value at stake not taken into account), drastically changed their ranking ( Fig. 4c and Fig. 4d). The most preferred FMAs (least susceptible) were the standard alternatives (M1, M2 and to a certain extent M4) which were characterized by the use of genetically improved material, planted at low density and with several thinning operations. The least preferred were the biomass alternatives (M6 and M7) which were characterized by high density planting and no or few thinning operations. Using hazard levels to weight susceptibility criteria hardly changed the ranking of FMAs (Fig. 4c vs. Fig. 4d).

Discussion
The methodological framework that we presented in this study represents a step forward for multiple risk analyses in forests. In particular, we improved the MCRA proposed by Jactel et al. (2012a) by providing more objective and accurate estimate of two components of risk, namely hazard level and the value at stake of forest (exposure value to damage).
The relevance of the improved hazards level evaluation appears more clearly when comparing the relative importance of the four main hazards agents shared by the three regions. Table 4 shows that storm importance was considered higher than fire in Aquitaine where landscape and climate is more favourable to fight against fires, whereas it is the other way round in Portugal. In Galicia, both storm and fire are the major concerns. This is consistent with what was recently observed in the regions: Aquitaine faced two major storms in the past 15 years (Colin et al., 2010), whereas Portugal and Spain are among the Southern European countries with the highest percentage of burnt forest every year since 1980 (San-Miguel-Ayanz et al., 2016). Moreover, Galicia was one of the two Spanish regions most affected by the Klaus storm (see the FORESTORMS database http://www.iefc.net/storm/ for details) and it ranks first in Spain for the surfaces affected by fires during the period 2001(ADCIF, 2012. Bark beetles were considered more important in Aquitaine than in Portugal, probably because past important storms in Aquitaine triggered two major bark beetles outbreaks (Aumonier & Maugard, 2006;DSF, 2010). Game and Fusarium circinatum are not of major concern in any of the regions. Sensitivity tests with the entire dataset showed that hazard levels (used as weights in the MCRA) had a negligible effect on the ranking of FMAs in terms of associated risks. This may be because we considered a high number of hazards agents together, which may have levelled out their individual effect, resulting in something similar to the equal weight option. This was confirmed by the ranking of FMAs based on their susceptibility to a number of hazards reduced to four (Fig. 3). In such situation, two FMAs were ranked worse, M8 in Aquitaine and M3 and M8 in Portugal. They have in common low susceptibility scores for the two hazards removed from the analysis (the insect defoliator and the root rot fungus for Aquitaine and the insect defoliator and the torrential rain for Portugal, respectively; see S3). This suggests that the number and the identity of hazards retained for the multi-criteria analysis might influence the overall risk ranking and that forest managers have to be consulted to set priorities based on their risk perception.
The sensitivity tests showed that when setting weights equal, i.e. considering equal hazard levels, the ranking of FMAs hardly changed (Fig. 4b) whereas when risk was reduced to susceptibility (i.e. without taking into account value at stake, Fig. 4c and Fig. 4d) the ranking of FMAs completely changed (Fig. 4c). Wood values at stake thus emerged as main drivers of risk ranking. The FMAs most at risk were those with highest values at stake, e.g. M1, M2 and M3 (Fig. 2). They were also the management alternatives with the longest rotation length (Table 1). This could be explained as the result of a longer period of time during which forest damage can occur, such as those due to occasional hazards (e.g. storms : Gardiner & Birot, 2013) or those that accumulate such as creeping hazards (e.g. insect defoliation or drought). Reciprocally, forest damage can be avoided if clearcutting was applied in the stand before the hazard occurs (Gardiner & Quine, 2000;Dhôte, 2005;Jactel et al., 2012b). Forest growing stocks are also higher in longer rotations, thus increasing the volume of damage (e.g. volume of windthrow) if affected by hazards (e.g. storms, Seidl et al., 2011a). Because forest values at stake were calculated by combining wood volume and wood price, the quality of products (timber vs. biomass energy wood) had also a great impact on standing value (Fig. 2). One might argue that when experts scored the silvicultural treatments for their susceptibility to hazards, they may have unconsciously overestimated those associated with FMAs designed to produce high value products. This was not the case, as demonstrated by our sensitivity analysis where values at stake were set equal: the FMAs most at risk (e.g. M1 and M3, Fig. 4a) were ranked as least susceptible (Fig. 4c).
Because new hazards may emerge in the future such as invasive pest species (Kenis et al., 2008;Seebens et al., 2017), or the occurrence of existing ones might change due to landscape or climate evolution, it remains important to adapt forest management in order to reduce stand susceptibility per se. The sensitivity test with no value at stake and hazard weights set equal revealed that the least resistant maritime pine plantations to a large array of hazards would be those produced by M6 or M7 (Fig. 4d). Their main differences with other FMAs are their high stand density, absence of thinning and their short rotation length. These silvicultural treatments were considered by experts to sharply increase stand susceptibility to fire due to their effects on fuel accumulation (Fernandes & Rigolot, 2007). High stand density is also known to reduce individual tree vigour, which may result in higher sensitivity to bark beetle attacks (Fettig et al., 2007). Reduced distance between neighbouring trees in dense pine plantations can also facilitate the transmission of root rot fungi through root contact (Garbelotto & Gonthier, 2013).
There are four main limitations to our study. First, the combination of heterogeneous sources and the lack of accurate and direct figures may have limited the accuracy of hazard assessment. It remains very difficult to obtain quantitative estimates of hazards occurrence due to the lack of a European-scale monitoring network where causes of damage are identified as well as their economic consequences. Still, since MCRA uses weights (relative levels of hazards) for criteria (FMA vulnerabilities to hazards), the use of absolute quantitative data is not imperative and qualitative data are as relevant as quantitative data, as long as scores are properly used to produce consistent relative ranking positions. Second, FMA susceptibility was estimated with the help of experts due to the lack of experimental studies on this topic, and yet only one panel per region could be gathered. As risk assessment also depends on risk aversion (Hanewinkel et al., 2011), it is difficult for experts to provide purely objective estimates. It would be advisable in the future to organize several panels of experts (Jactel et al., 2012a) and use structured forecasting approaches (e.g. the Delphi method: MacMillan & Marshall, 2006;Salas-Garita & Soliño, 2019) to elicit expert opinions so as to come to a consensus on risk rating. We can also use quantitative technics such as Saaty's pairwise comparison (see in Diaz-Balteiro et al., 2016) provided that the number of experts are large enough to do so. Third, we could only estimate stand value at stake for even-aged stands using available growth models for maritime pine, preventing us from testing other FMAs based on uneven-aged forestry even though we assume that this type of structure would perform better in case of hazard. Fourth, the susceptibility of a given stand might depend on that of neighbouring stands, particularly for hazards with contagion dynamics, like storms or fire. It would be then advisable to upscale forest risk analysis to the landscape level.

Conclusions
Our MRCA methodology provides a solid starting point for more detailed studies incorporating the concept of multiple risks in forest planning when no modelling tools are available.
First, we highlighted the great importance of forest value at stake for ranking FMAs according to vulnerability to multiple hazards. We showed that short rotation forestry is less prone to risk due to a reduced probability of being hit by large disturbances with long return periods (e.g. storm, fire, bark beetles) but, above all, to the lower economic value of standing growing stock (and also the lower operational costs for reforestation). However, wood biomass production is also probably less profitable for forest owners due to reduced revenues (Schwarzbauer & Stern, 2010) so forest owners have to combine risk aversion and gain expectation. In addition, managers may consider paying insurance to cover losses due to forest damage. The next step would be thus to combine risk and profitability as criteria to rank FMAs and thus propose more realistic advice to forest managers.
Second, we restricted our analysis to pure stands, focused on wood production and without any interaction with the other components of the landscape (hazard propagation, benefits of surrounding patches of diverse forests...). However there are other silvicultural options that could reduce forest risks while maintaining high productivity, in particular the use of tree species mixtures (Zhang et al., 2012;Gamfeldt et al., 2013;Jactel et al., 2017;Jactel et al., 2018). Finally, to compare more complete sets of forest management alternatives, more criteria related to forest value could be included in the decision-making process such as the provision of ecosystem goods (de-Miguel et al., 2014) and services (Gamfeldt et al., 2013;Branco et al., 2015), which are increasingly taken into account when assessing forest profitability (Deal et al., 2012). Forest price volatility could be also considered as additional hazard.
A way forward to test more combinations of silviculture treatments and to identify FMAs reducing risks while increasing the provision of services will be to design forest growth models  able to model stand evolution under several types of forest management and hazard dynamics. This will require knowing trends in hazard occurrence in a given region and identifying linking functions between stand features (including species choice and composition), resistance to hazards (Pukkala et al., 2005;Seidl et al., 2007) and delivery of services (Lindroth et al., 2009). Model outcomes could be then used as criteria in multi-criteria decision analyses where forest stakeholders would be invited to weight the criteria according to the ecosystem services they expect from forests, with the objective to identify tradeoffs and synergies.