Spatial Optimization and Tradeoffs of Alternative Forest Management Scenarios in Macedonia, Greece

: Managing forests has been demonstrated to be an efﬁcient strategy for fragmenting fuels and for reducing ﬁre spread rates and severity. However, large-scale analyses to examine operational aspects of implementing different forest management scenarios to meet ﬁre governance objectives are nonexistent for many Mediterranean countries. In this study we described an optimization framework to build forest management scenarios that leverages ﬁre simulation, forest management, and tradeoff analyses for forest areas in Macedonia, Greece. We demonstrated the framework to evaluate ﬁve forest management priorities aimed at (1) protection of developed areas, (2) optimized commercial timber harvests, (3) protection of ecosystem services, (4) ﬁre resilience, and (5) reducing suppression difﬁculty. Results revealed that by managing approximately 33,000 ha across all lands in different allocations of 100 projects, the area that accounted for 16% of the wildﬁre exposure to developed areas was treated while harvesting 2.5% of total wood volume. The treatments also reduced fuels on the area that are responsible for 3% of the potential ﬁre impacts to sites with important ecosystem services, while suppression difﬁculty and wildﬁre transmission to protected areas attainment was 4.5% and 16%, respectively. We also tested the performance of multiple forest district management priorities when applying a proposed four-year fuel treatment plan that targeted achieving high levels of attainment by treating less area but strategically selected lands. Sharp management tradeoffs were observed among all management priorities, especially for harvest production compared with suppression difﬁculty, the protection of developed areas, and wildﬁre exposure to protected areas.


Introduction
The recent high death toll from wildfires burning into developed areas in Greece (i.e., in 2007 and 2018) highlights the importance of forest management for wildfire risk mitigation for community protection. This is in sharp contrast to the 1950s-1970s when timber production and management with traditional agroforestry were the highest priorities [1][2][3]. Recent extreme fire events have motivated managers and policymakers to further re-examine fire policies in order to develop comprehensive and strategic forest management programs, allowing them to adopt the best compromise between fire control and fuel treatment approaches while considering the fact that full suppression of wildfires is not a feasible and reasonable strategy in the long term [4,5]. One of the most important challenges for forest management agencies is to consider and decide among competing landscape management priorities, and such a challenge creates a complex spatial tradeoff problem for managers that seek to identify optimal arrangements within economic and other constraints [6]. Developing a broader mix of forest management priorities that are tailored to particular landscapes based on fire regimes, human values, and land use could potentially highlight where alternative and integrated strategies provide a long-term solution to better coexist with fire [5]. Management goals can include, among others, fire resiliency in forested landscapes, the restoration of the cultural fire regime, a safe and efficient fire response, and fire-adapted developed areas (see the cohesive strategy used in the United States [7]).
In Greece, 6.5 million hectares are either characterized as forests (30% of the total Greek territory) or forested lands (20%). During 1990-2020, forests increased by 18.3%, but timber retrieval fell from around 2.5 million m 3 in the 1990s to 1.4 million m 3 in 2016, with the decrease mostly attributed to fuelwood production. In total, timber stock increased from 156 million m 3 in 1990 to 185 million m 3 in 2013, with a stable average growing stock volume per hectare of 47 m 3 [8]. A combination of an increase in forest cover, a decrease in annual timber retrieved, and the subsequent accumulation of dead and live vegetation indicates a need for an increase in the area treated annually to target reducing the spread rate and intensity of fires that can potentially harm developed areas, cultural monuments, ecological assets, and ecosystem services. The problem is intensified by the decrease in employment in the forest management sector, which was reduced by half between 2009 and 2015 (i.e., around 23 thousand people) [3,8]. However, this needed increase in forest management is unlikely to come from the few hundred private forest owners that manage 8% of all forests in Greece, as these are primarily coppice enterprises producing mostly fuelwood with very low profitability and limited part-time employment, mostly in rural areas [3]. From the above, it is clear that the Greek state should be more actively engaged in forest management and promote policies which can eventually increase the annual area that is either treated for fire risk reduction or managed for achieving ecological or restoration objectives, while enhancing cooperation between the major land tenures (private, municipal, state, and church).
The existing laws, regulations, and constitution articles to sustain forests and promote biodiversity conservation make it more difficult to accelerate and expand the rate and scale of forest management. For example, all Greek forests and forested areas are protected by the Greek Constitution (Article 24) [9], which dictates that forest land use changes are prohibited unless they are required for public interest. Deforestation activities are limited and permitted only in specific cases for the public interest and benefit (e.g., construction of roads, railways, high tension lines), following the direct administrative procedures under the provision of Greek laws (Legislative Decree No. 86/1969 [10]; Laws No. 998/1979 [11] and No. 1734/1987 [12]). The concept of sustainable yield (Administrative Regulation No. 120094/499/1937 [13]) includes assumptions regarding forest protection from fire and other disturbances, ensuring regeneration after timber retrieval or disturbance, secure forest growth by replacing mature and over-mature forests with slow or no net growth rates, and optimal timber retrieval that allows the forest to evolve and supply at least the same or greater timber volume in the future. Furthermore, the Administrative Regulation No. 10223/958/1953 [14] provides guidelines for the implementation of forest management plans for state and private forests, and the Legislative Decree No. 86/1969 [10] sets very strict regulations regarding forest management for both public and private forests.
For more than three decades, Law No. 998/1979 [11] has regulated the protection of the country's forest and other wooded lands, providing the legislative framework to apply forest management following specific rules and guidelines for practices driven by the fundamental principle and predominant goal of preserving and promoting the "sustainability" of forests in terms of their provision of products, growing stock, and services. Any temporary loss of forest cover is not considered deforestation and is declared reforested to recover its former state (article 61 Legislative Decree No. 86/1969 [10], articles 37, 38, 46, 47 of Law No. 998/1979 [11]). Over the years, there have been a number of biodiversity conservation reforms (European Union Directives 92/43/EEC [15] and 2009/147/EC [16]), including imposing restrictions on clearcutting; prohibiting felling of trees with nests; calling for the retention of overmatured, malformed, dead trees at In the second scenario, given the assumption that we follow a business-as-usual approach to manage only inside forest districts, we investigate whether an accelerated fouryear treatment plan with an annual pace (i.e., instead of the established 10-year treatment plan) can be applied to achieve a substantial attainment for each of the five management priorities. This scenario is more localized and attempts to identify where the best projects can be created across all forest districts and what the expected management outcomes are for each management priority.
In the third scenario, given the assumption that we can manage the entire available fortreatment area of each forest district, we investigate which of the four non-timber priorities has the best joint production with timber volume retrieval and where these treatments are located.

Overview
In brief, we partitioned the study area into hexagon stands (n = 81,000), 37,000 of which were within the 290 forest districts, where potential fuel treatment projects were tested and evaluated. Stands were populated with spatial data derived from multiple sources including fire simulations and other assessments (e.g., growing stock volume, management status, socioeconomic attributes) to map and quantify five management priorities (see Table A1 in Appendix B). We then used a scenario planning model [21] to maximize each of the five management priorities under the three different modelling scenarios. These scenarios were designed to explore different strategic approaches to allocate fuel treatment projects within the study area in terms of efficiency and tradeoffs among the priorities. We assumed that if a stand were treated, fire exposure is then mitigated.

Study Area
Our study area covers 34,000 km 2 across three regions of Greece (i.e., Western Macedonia; Central Macedonia; Eastern Macedonia and Thrace), further divided into 13 prefectures encompassing seven national forests and parks and 290 forest districts ( Figure 1A). Forest administration and the application of forest management plans are achieved by 25 Greek Forest Service local branches, excluding the autonomous monastic region of the Mount Athos peninsula, which is a separate entity that is further divided across 20 monasteries; these monasteries independently manage their lands based on forest management plans that have been approved by the Greek Forest Service. The population of the study area is over 2.5 million people residing in more than 2000 cities, communities, and villages.
The Greek Macedonian landscape is characterized by large plains irrigated by long rivers, coastal areas with low elevation vegetation types, mountainous areas with productive conifer and broadleaf forests mixed with grasslands, and alpine zones. Currently, most low elevation forests are high in density with ladder fuel structures and have a high potential for stand-replacing, high severity fire events. The climatic factors in the study area controlling large fire weather conditions present substantial region-wide differences [22].
According to the 2018 Corine Land Cover data [23], private agricultural lands cover 30% of the study area, followed by broadleaved forests (20%), transitional woodland-shrub (9%), pastures (8%), and land principally occupied by agriculture with significant areas of natural vegetation (8%) ( Figure 1A). Artificial surfaces and water or wetlands make up to 3% and 2.5% of the total study area, respectively. Sclerophyllous vegetation or shrublands and sparsely vegetated areas cover 7.5% and 2% of the study area, respectively. Almost 6% of the study area is covered by conifer species, half of which are low elevation (Pinus halepensis, Pinus brutia, Pinus pinea, and Juniper spp.), and the rest dominate mountainous regions that include mostly species such as Abies spp., Picea abies, Pinus nigra, Pinus sylvestris, and Pinus heldreichii. Finally, mixed conifer-broadleaf forests cover 5% of the study area, mostly comprising Fagus sylvatica mixed with Abies spp. or Pinus nigra or Quercus petraea, Q. frainetto or Q. pubescens.

Forest Districts and Stands
Although forest districts have defined stand boundaries, as indicated within forest management plans, lands outside forest districts do not have such information. To account for this, we partitioned the entire study area into individual forest stands by overlaying a hexnet layer, while considering the intersections with managed forest boundaries, major land cover types, developed area boundaries, and protection status ( Figure 1B). Stands functioned as spatial decision units to simulate management [29]. First, we created a hexnet for the entire study area with 53,000 hexcells (stands), all with an area of 65 ha. Stand size became variable, and the stand number increased due to the intersections, which partitioned the original hexnet into 81,000 stands with sizes ranging from 10 to 65 ha (mean of 42 ha). This better captures the spatial gradients in treatment priority metrics across the landscape.
Forest management plans usually disregard degraded stands that have existed for one or more decades and are not receiving any official management. These stands are mostly grazing landscapes or places that have a protective role in soil stabilization and erosion aversion, or they have low timber stock volume. Protected area cores or virgin forests are also left unmanaged, targeting either conservation or biodiversity. To account for this, all stands that fell inside one of the protected areas (national forest cores, aesthetic forests, national parks, and nature protection areas) or were not forests or forested areas (e.g., grasslands) were flagged as administratively unavailable (nonmanageable). In addition, we flagged as nonmanageable those stands with more than one-third of their area either treated or disturbed by wildfire, insects, or disease during the 2010-2019 period [30]. All stands that were coniferous, broadleaved, or mixed forests, have sclerophyllous vegetation and transitional woodland-shrub [23], and with an average growing stock volume >35 m 3 /ha (see next sections for additional details and [31]), were available for treatment. Stands were then assigned to the 290 pre-existing forest districts (average size = 5185 ha, minimum = 514 ha, maximum = 66,750 ha) ( Figure 2A), with half of all stands not belonging to any forest district.
x FOR PEER REVIEW 9 of 30 especially in forested stands with low canopy base heights and adequate canopy fuels, can reduce crown fire severity potential and increase future fire resistance. FConstMTT produced gridded text files with information about the burn probability and fire intensity for each point on the landscape. Fire intensity estimates were calculated using the 20 Fire Intensity Level (FIL) classes with 0.5 m increments, up to the maximum value of 9.75 m. To quantify wildfire hazard, we used conditional flame length, which is the probability-weighted flame length given a fire occurs (Equation (1)): where BPi is the probability of a fire at the i-th flame length category (the FIL value), and FLi is the flame length midpoint of the i-th FIL category (e.g., 0.25 m, 0.75 m, etc.). We calculated the CFL value in meters at the stand level by averaging all pixels whose centers intersected the stand. Overall, the CFL is very high (>3 m) on 14% of the study area, denoting an increased complexity for fire suppression efforts that require indirect firefighting tactics ( Figure 2E). The majority of the study area (77%) has the potential to experience low flame lengths (>0-2 m), meaning that ground crews aided by machinery can confront such a fire event. Finally, 9% of the study area has the potential to experience fires with moderate CFL (2-3 m), requiring aerial means in collaboration with ground forces for successful fire containment. CFL was used to prioritize stands for treatment to reduce suppression difficulty.

Stochastic Fire Behaviour Modelling
Stochastic fire behavior simulations were conducted with the command line version of the minimum travel time (MTT) fire spread algorithm proposed by Finney [32], as implemented in FConstMTT [33]. The necessary data layers describing the fuels and topography of the study area were assembled in a gridded landscape file with a 100 m spatial resolution binary file. The required layers include elevation, aspect, and slope to describe the topography, fuel models to describe surface fuel quantity and arrangement, canopy cover, canopy bulk density, canopy base height and canopy height of the wooded areas to simulate crown fires, all of which were retrieved from the Copernicus database, or the Greek Forest Service, or other existing data sources and previous research efforts (see [34]).
The study area was divided into 11 fire regime macro-areas that capture the changing fire activity gradients based on historical fire activity and the climatic stratification of the environment as described by Metzger [35]. Analyzing fire activity separately on these macro-areas facilitated the segmentation of the study area into 30 smaller zones (min. at 30,000 ha; max. at 280,000 ha; mean at 140,000 ha) with differing wildfire season duration (i.e., annual period concentrating 90% of the burned area from fires >50 ha) and very particular burn patterns associated with local fire-weather variability [36].
For each zone, we assigned at least one representative remote automatic weather station (RAWS) with at least 10 years of recorded data (see [34]). For each RAWS, we estimated the frequency of the dominant wind direction and used it as the base of each simulation scenario, with the frequency translated into a scenario selection probability. The average wind gust speed for each RAWS was used as the base wind speed for all simulations since large human-caused wildfires coincide with higher wind speeds [37]. Lastly, we set a base simulation duration of 300 min and a spotting probability of 30%. All weather estimates covered only the months of July and August.
In total, 10,000 fires were simulated for each zone, i.e., 300,000 for the entire study area. During fire modelling, each fire was independently modelled, and weather conditions were held constant, while fire suppression efforts were not considered due to relatively limited containment capabilities during extreme fire events. To calibrate the surface fire spread model, we matched historical with simulated large fire size (>50 ha) distributions in every macro-area. The FConstMTT simulations generated (1) fire perimeters, (2) annual burn probabilities (BP), and (3) conditional flame length (CFL). The annual BP is the ratio of the number of times a pixel burned to the total number of fires simulated, and CFL measures the expected flame length in meters given a pixel burn.

Protection of Developed Areas
Protecting communities and structures in the wildland-urban interface (WUI) is the primary goal of both fire suppression and annual fuel treatment programs that involve the application of actions such as fuel breaks, fuel removal, and to a lesser extent, thinning in the periphery of developed areas and along major roads [20,38]. Recently, funding for such treatments has increased (see Introduction) with the protection of developed areas being the highest management priority for the Greek Forest Service, in collaboration with local authorities.
We geometrically intersected all simulated fire perimeters with a polygon layer depicting the major land tenures and communities of the study area. To estimate the predicted exposure to developed areas, we assumed that the structures reported in Hellenic Statistical Authority's 2011 census data for each community polygon are spatially distributed equally, and the percentage of burned area from each fire within each polygon was translated into the annualized number of structures affected.
The simulated fires were used to delineate the fireshed around developed areas where large fires are likely to ignite and/or spread, and were attributed with the conditional total structures affected, assuming a fire occurred at the location. We applied the ArcGIS Inverse Distance Weighting (IDW) process on all simulated fire ignition locations with a structure exposure estimate >0 with a 1 km search radius; we then interpolated the estimated total structure exposure of each simulated fire ignition, subsequently smoothing it using a focal mean process. The extent of the fireshed for the approximately 2200 developed areas is 9250 km 2 , i.e., 27% of the total study area ( Figure 2B). Finally, for each stand we estimated the average value of the fireshed of all pixels intersecting it; this value represents the number of structures exposed to wildfire from ignitions within the stand. Each stand can have many ignitions that affect multiple developed areas. Structure exposure was the metric used to prioritize protection of developed areas.

Forest Production
Prioritizing forest production identifies the locations where treatments can generate revenue that in turn can be used to subsidize noneconomic fuels treatments [6]. In addition, rural communities can also participate in the projects by retrieving fuelwood for household purposes. Accurate estimates of the timber volume available for harvesting, including the type and quality of the retrieved stocks, are only available for forest districts. Some examples of the growing stock of common tree species found in the area are as follows: P. halepensis and P. brutia have a growing stock of about 26 m 3 /ha, P. nigra has about 54 m 3 /ha, P. sylvestris has about 123 m 3 /ha, and Picea abies has about 342 m 3 /ha [39]. One major drawback is that most of these plans were either not available to us or outdated (i.e., more than a decade old), and thus additional processing was required in order to approximate the level of current timber volume in each stand.
We retrieved the raster layer (≈100 m cell size) of growing stock volume (GCV) for the year 2010 [31,40]; this layer mapped the volume (m 3 /ha) of all living trees more than 10 cm in diameter at breast height, measured over bark from ground or stump height, excluding smaller branches, twigs, foliage, flowers, seeds, stump, and roots ( Figure 2C). The GSV estimates were obtained from spaceborne SAR (ALOS PALSAR, Envisat ASAR), optical (Landsat-7), LiDAR (ICESAT), and auxiliary datasets with multiple estimation procedures [40]. To account for forest losses from 2010 onward, we used the Global Forest Change dataset [30], which is a time-series analysis of Landsat images that describes a stand-replacement disturbance, i.e., a change from a forest to nonforest state, during the period 2000-2019. These losses can result from infestations, weather incidents, wildfires or timber harvesting, and silviculture, and all pixels flagged as "forest loss" after 2010 were used to zero out the values of GCV pixels. For each stand we assigned the average value of all pixels whose centers intersected the stand, using zonal statistics in ArcGIS, and we then multiplied this value (m 3 /ha) with the stand area (ha) to retrieve the stand GSV. This value represents the volume of all trees of the stand, regardless of their commercial value, their stand structure, and condition, and it was the metric used to prioritize forest production.

Ecosystem Services
Prioritizing the protection of ecosystem services from wildfire through fuel treatments can provide manifold improvements to current forest management efforts in order to mitigate a loss of revenue in fire-prone forested and rural areas. These areas are of critical importance to the well-being and financial stability of local populations. Since most of the country's forests are low productivity, the forestry sector alone contributes less than 0.2% to Greece's gross domestic product. Recent studies used additional metrics to holistically estimate the value of the forest ecosystem functions in Greece, which was approximated at 56 billion EUR (2 billion EUR per year) [39,41]. Recreation is the most important economic function (27% of the total value), followed by grazing (23%), biodiversity (21%), and soil protection (18.5%), while timber value contributes only 5% to the total land value [41]. Other economic functions include carbon sequestration, hunting, and nonwood forest products.
We retrieved an official spatial estimation of the economic value of land in Greece, ordered by the Ministry of Environment and the Greek Forest Service and produced by the government research organization ELGO-Dimitra [42]. To create this product, economic functions and values of forested lands, including pastures and shrublands, were separately evaluated. In addition, negative externalities from wildfires, soil erosion, and floods were included in the estimation. The outcome was a 100 m cell size raster dataset in which the total land value was estimated at euros per one tenth of a hectare (the Greek unit of area measure, stremma, i.e., 1000 m 2 ) ( Figure 2D). We first estimated the average value at the stand level from all pixels whose center intersected the stand and then multiplied it by the area of the stand adjusted to stremma. The higher the stand land value is, the more important it is for us to apply fuel treatments on this stand to protect its functions from losses due to potential wildfire. The stand land value was used to prioritize the protection of ecosystem services.

Suppression Difficulty
Flame length is the distance measured from the average flame tip to the middle of the flaming zone at the base of the fire; it is measured on a slant when the flames are tilted due to effects of wind and slope and serves as an indicator of fireline intensity. Despite serious considerations on how flame length can be accurately measured in the field [43], it is very apparent to fireline personnel, and so it can readily convey a sense of fire intensity that it is worth featuring as a primary variable representing the potential wildfire hazard on the landscape [44]. Higher flame length values represent fuels with a higher probability of experiencing torching, crowning, and other forms of extreme fire behavior under conducive weather conditions. Treating stands with highest flame length values, especially in forested stands with low canopy base heights and adequate canopy fuels, can reduce crown fire severity potential and increase future fire resistance.
FConstMTT produced gridded text files with information about the burn probability and fire intensity for each point on the landscape. Fire intensity estimates were calculated using the 20 Fire Intensity Level (FIL) classes with 0.5 m increments, up to the maximum value of 9.75 m. To quantify wildfire hazard, we used conditional flame length, which is the probability-weighted flame length given a fire occurs (Equation (1)): where BP i is the probability of a fire at the i-th flame length category (the FIL value), and FL i is the flame length midpoint of the i-th FIL category (e.g., 0.25 m, 0.75 m, etc.). We calculated the CFL value in meters at the stand level by averaging all pixels whose centers intersected the stand. Overall, the CFL is very high (>3 m) on 14% of the study area, denoting an increased complexity for fire suppression efforts that require indirect firefighting tactics ( Figure 2E). The majority of the study area (77%) has the potential to experience low flame lengths (>0-2 m), meaning that ground crews aided by machinery can confront such a fire event. Finally, 9% of the study area has the potential to experience fires with moderate CFL (2-3 m), requiring aerial means in collaboration with ground forces for successful fire containment. CFL was used to prioritize stands for treatment to reduce suppression difficulty.

Wildfire Transmission to Protected Areas
If an area is characterized with a protection status in Greek or European Union legislation (Wetlands-Ramsar; Specially Protected Areas under the Barcelona Convention; Biogenetic Inventories; Biosphere Reserves; World Heritage sites; NATURA 2000; National Parks and Forests), limitations on human interventions are imposed, including fuel and forest management. Restrictions can include limiting the quantity and quality of retrieved timber, allowing for only specific mechanical treatment methods to be applied, controlling the opening of new roads, and forbidding the use of machinery, equipment, and vehicles. The role of protected areas is important in maintaining ecosystem services, ecological balance, aesthetic quality, and biogenetic reserves. A wildfire can impact the functioning of protected areas, and years can be required to restore them to their prefire state.
Extensive parts of the study area have some type of protection status, but this does not necessarily mean that fuel and forest management is restricted. There are more than 200 protected areas that cover approximately one quarter (900,000 ha) of the study area. We kept the 15 most important protected areas that include National Forest cores, aesthetic forests, National Parks, and nature protection areas where management and human intervention are limited to a higher degree or totally prohibited, covering approximately 190,000 ha.
Similar to the protection of developed areas, our goal was to understand where fires initiate and burn into these areas. To achieve this, we geometrically intersected all simulated fire perimeters with a polygon layer depicting the major land tenures and protected areas. This enabled the estimation of outgoing fire from each ignition that burned into protected area polygons. Then, we applied the IDW process with a 1 km search radius on all ignitions that affected protected areas, with annual hectares burned >0. Finally, for each stand we estimated the average value of all pixels intersecting it ( Figure 2F). Each stand can have many ignitions that affect multiple protected areas. Higher stand level protected area exposure denotes a higher priority for receiving fuel treatments.

Multicriteria Spatial Prioritization and Optimization
To design and simulate landscape forest management prioritization scenarios and tradeoffs, we used the multicriteria, hierarchical planning model ForSys [29], which was previously known as the Landscape Treatment Designer [6,21,45,46]. The model identifies the treatment units that maximize the attainment levels for multiple priorities with varying priority weights, subject to treatment constraints (e.g., budget or treatment limits), treatment thresholds (e.g., forest stands susceptible to high severity prescribed fire) and legallymandated exclusions to treatment (e.g., excluding protected areas). ForSys has been applied in several prior case studies in the Mediterranean [47][48][49] and elsewhere [29,46,50,51].
ForSys uses spatial information in the form of shapefiles that describe stands (hexcells) that can be also grouped into predefined operational planning areas (i.e., forest districts in our study). The ForSys algorithm uses a relatively simple search heuristic that tests each polygon as a seed to build a restoration project in the surrounding landscape; it absorbs adjacent stands based on their potential contribution to the priority value and treats those that meet the treatment thresholds [46]. The user supplies a forest management scenario in terms of priorities, treatment constraints, and treatment thresholds, and the program identifies one or more project areas within the landscape that maximize the priority [46]: This is subject to: where C is a global constraint on the investment level per project area (e.g., area treated), Z is a vector of binary variables indicating whether the j-th stand is treated (e.g., Z j = 1 for treated stands and 0 for untreated stands), N ij is the contribution to objective i in stand j if treated, and A is the area of the j-th treated stand. W i is a weighting coefficient that can be used to emphasize one objective versus another. The flowchart of program control for ForSys showing the decision framework for simulations that use the aggregation option can be found in Appendix A ( Figure A6). The fundamental idea is that if a substantial portion of a forest district is treated on an annual basis, it has a positive effect on reducing fire intensity for the entire forest district, while low treatment intensities or scattered projects have little or no influence on limiting large fire spread [48,52]. Another assumption is that each restoration priority would be addressed by thinning for overstocked stands and surface fuel reduction with mastication or fuel removal. Stands were added to each project based on the mean value of the priority for that stand, and when the area treatment constraint was met (projects of 100, 250, and 500 ha, depending on the modelled scenario), the total priority value was summed. ForSys runs were applied in stands with merchantable volume greater than 35 m 3 /ha and were set to be aggregated, meaning that there were spatial constraints in having the selected stands with common boundaries.

Linking Stands with Fire Behavior Modelling Outputs and Other Forestry Metrics
Each stand was populated with the average values or binary flags from variables derived from fire behavior modelling or other ecological and timber-related data sources (e.g., protection status, available growing stock volume, disturbances, etc.). Due to the various stand shapes and sizes, each variable's raw value had to be converted to stand value with Equation (4): where S ij is the stand value for stand i and for variable j, C ij is the raw averaged or summed value for stand i and for variable j, and A i is the area (hectares) of stand i. In addition, each stand had a unique identifier, centroid coordinates, and an area field. Then, the percentage contribution with respect to the total problem of all treatment units (i.e., PCP) was calculated to standardize reporting across all priorities and to assess the attainment degree of all treated units with Equation (5) where PCP ij is the percent contribution to the total problem (sum of the data subset K) of stand i and for variable j. These transformations were required for the subsequent scenario planning analysis. Finally, to assess the importance of each stand that defined whether it will be selected or not in a project, we estimated the percentage difference from the maximum value of all stands of the hexnet for each priority's stand value with Equation (6): where SPM ij is the percentage difference from the maximum value of data subset K (landbase) of stand S i and for variable j. PCP and SPM for each of the five priority values were calculated for all stands of each landbase hexnet.

Modelled Treatment Scenarios for Each Management Priority
We simulated three forest management scenarios that explore all five management priorities: (1) prioritize forest management projects across the entire landscape; (2) develop a four-year fuel treatment plan restricted to the 290 forest districts aiming to treat the best projects until one-third of the attainment achieved by all simulated project is met; and (3) understand and quantify potential tradeoffs and joint production between forest production and the other four priorities inside forest districts. The overview of the key run parameters used in ForSys, the postprocessing approach, and results for each prioritization scenario are in Table A1 in Appendix B. Similar to our previous studies [29,51], treatment within all scenarios was a combination of thinning (either commercial or noncommercial depending on stand conditions) and broadcast/pile burning based on silvicultural prescriptions. We assumed that if a stand were treated, fire exposure is then mitigated. The current forest management program through approved forest management plans treats on average 41,000 thousand ha to retrieve 530,000 m 3 of timber annually within the study area during the past ten years (Tables A2 and A3), with the largest treated area occurring in 2019 with approximately 46,000 ha and the lowest in 2011 with approximately 30,000 ha.

Landscape-Scale Project Prioritization (Scenario 1: All Lands)
Our first scenario generates 100 fuel treatment projects, each with a maximum treated area of 500 ha (total treated area of approximately 50,000 ha); this scenario prioritizes each of the five management priorities independently. Projects were prioritized by ranking the attainment achieved for each priority run and mapped. We assessed cumulative attainment by increasing the area treated for each priority relative to the other four priorities. We also identified which stands were selected in the projects of more than one priority, and we estimated the cumulative attainment of all five priorities if those stands are treated.

Forest District Project Prioritization (Scenario 2: Annual Treatments)
In our second scenario, we simulated fuel management scenarios to find the best five 100 ha projects inside each of the 290 forest districts; this was performed for each of the five management priorities. Approximately 1300 projects were created and ranked for each priority. Then, we set a goal to select and treat within four years the best projects of each priority until one-third of the attainment achieved by all simulated projects was reached. In each of the four years, the percentage attainment was approximately the same, but the amount of treated area varied. Thus, the projects were selected each year based on their attainment rank.
To achieve the annual attainment goal, projects could be located across all forest districts or not, and this depends entirely on the project attainment. For example, if the goal was to apply fuel treatments to treat the total developed areas exposure by 16% over the entire four-year period, during the first year it is possible that 4% could be treated with only five projects on three forest districts-which are considered the best of all available projects-while in the fourth year for the same attainment, it is necessary to treat ten projects across five forest districts. Finally, we estimated the attainment achieved for each of the five priorities at the end of the four-year treatment plan if the entire area that contributed to all projects from all priorities were treated.

Forest District Tradeoffs (Scenario 3: Tradeoffs)
For the third scenario, we simulated 250-ha projects until the entire forest district was treated (for all 290 districts) by prioritizing paired combinations of forest production (i.e., growing stock volume) and each of the other four management priorities by changing the relative weights of each objective. We used different sets of integer weights in all combinations from 0 to 4 in increments of 1, resulting in 14 outcomes for each pair. For instance, weights of 1 and 0 for priority A and B respectively represent the maximum production for priority A, whereas weights of 2 and 2 for each priority represent a mixed production for both [46]. For each pair of priorities, we ranked each forest district by joint attainment of both priorities, and we used the 15 highest ranked forest districts to create production possibility frontier graphs between growing stock volume attainment and the attainment of each of the other four priorities. In these graphs, lines represent the forest districts and the tradeoffs between the two priorities, with individual points on each line representing different project solutions within that district. The more right a point on the line is located on x-axis, the greater the attainment in growing stock volume, while the higher a point is on the y-axis, the greater the attainment of the paired priority. The shape of the line is indicative of the relationship between the two priorities, with convex shapes representing decreasing opportunity cost; in other words, as we move along the curve with increasing growing stock attainment, less attainment is achieved for the other priority.

Scenario 1: All Lands
For each priority we ranked and mapped the best 100 projects (approximately 500ha each) based on attainment and assessed the cumulative attainment and area treated for each objective (key results can be found in Table A1, Appendix B). The line with the highest attainment was always not surprisingly the selected priority, except for suppression difficulty where treated projects resulted in higher attainment for wildfire transmission compared to protected areas. We also filtered the stands from the five prioritization runs to identify those that were selected more than once to create new projects that can achieve multiple objectives. Although this process does not guarantee spatial colocation of the new projects, it resulted in more than 100 spatially aggregated projects (mean area: 330 ha; max. area: 1485 ha) that include 94% of all stands with >1 priority met. In total, 5050 stands from all five prioritization runs were identified, covering 210,000 ha ( Figure 3A; purple, blue, and red). We found 722 stands (30,200 ha) that were selected for two priorities ( Figure 3A; blue), and 72 stands (2350 ha) that were selected for three priorities ( Figure 3A; red). If management is applied on these ≈800 stands ( Figure 3A; upper left graph), the attainment is higher for the protection of developed areas (16%) if 15,800 ha from the total of 32,550 ha are treated, followed by the wildfire transmission to protected areas priority (16%) if 19,000 ha are treated. Attainment was found to be lower for suppression difficulty (4.5% if 22,200 ha are treated), and approximately 2.5% for forest production and 3% for ecosystem services, if all 800 stands were treated.
The cumulative attainment of the 100 projects for the protection of developed areas priority was 47% if we treat 43,000 ha ( Figure 3B). This level of attainment is much higher compared to any attainment achieved by the other four prioritization runs, suggesting that small parts of the landscape account for most of developed areas fire exposure, making it easier to manage less land in order to achieve high attainment. In addition, on the priority projects for the protection of developed areas, we can achieve a 10.5% attainment for wildfire transmission to protected areas, 6.5% attainment for suppression difficulty, 2.3% attainment for forest production, and 3.2% attainment for ecosystem services. The best projects are in clusters close to the capital city of Thessaloniki and in coastal areas of the Prefectures of Kavala, Chalkidiki, Pieria, and Thassos Island.
The cumulative attainment of the 100 projects for the forest production priority was 7.4% if we treat 45,000 ha ( Figure 4A). Regarding the other four priorities, we can achieve a 3.9% attainment for ecosystem services, 0.8% attainment for suppression difficulty, 0.4% attainment for wildfire transmission to protected areas, and 0.1% attainment for developed areas exposure. The very low attainment for the other priorities suggests that joint production is not possible on the selected stands. It is notable that the best projects are clustered in three prefectures (i.e., prefectures of Pella, Florina, and Imathia) where productive high elevation conifer and broadleaved forests prevail.
The cumulative attainment of the 100 projects for the ecosystem services priority was 5.4% if we treat 45,300 ha ( Figure 4B). Regarding the other four priorities, we can achieve a 4.9% attainment for forest production, 3% attainment for developed areas exposure, 2.5% attainment for suppression difficulty, and 1.4% attainment for wildfire transmission to protected areas. The attainment of other priorities is substantially higher compared to that achieved when forest production was prioritized, suggesting that joint production is possible on some projects, although not to a substantially high degree. The best projects are clustered in the first peninsula of Chalkidiki, in the western part of the study area (prefectures of Kastoria and Kozani), and at the prefecture of Pella. selected for two priorities ( Figure 3A; blue), and 72 stands (2350 ha) that were selected for three priorities (Figure 3A; red). If management is applied on these ≈800 stands ( Figure  3A; upper left graph), the attainment is higher for the protection of developed areas (16%) if 15,800 ha from the total of 32,550 ha are treated, followed by the wildfire transmission to protected areas priority (16%) if 19,000 ha are treated. Attainment was found to be lower for suppression difficulty (4.5% if 22,200 ha are treated), and approximately 2.5% for forest production and 3% for ecosystem services, if all 800 stands were treated. The cumulative attainment of the 100 projects for the protection of developed areas priority was 47% if we treat 43,000 ha ( Figure 3B). This level of attainment is much higher compared to any attainment achieved by the other four prioritization runs, suggesting that small parts of the landscape account for most of developed areas fire exposure, making it easier to manage less land in order to achieve high attainment. In addition, on the priority projects for the protection of developed areas, we can achieve a 10.5% attainment for wildfire transmission to protected areas, 6.5% attainment for suppression The cumulative attainment of the 100 projects for the suppression difficulty priority was 10% if we treat 45,000 ha ( Figure 5A). This attainment level suggests that reducing wildfire hazard-and thereafter suppression difficulty-is not an easy task and requires more extensive management to make a substantial difference for firefighting forces. Regarding the other four priorities, we can achieve 14.5% attainment for wildfire transmission to protected area, 7% attainment for developed areas exposure, 3.3% attainment for ecosystem services, and a 2.2% attainment for forest production. The best projects are in the prefectures of Chalkidiki, Kozani, and Kavala. Finally, the cumulative attainment of the 100 projects for the wildfire transmission to protected areas priority was 36% if we treat 44,300 ha ( Figure 5B). This is the second highest attainment achieved from the five priorities. Regarding the other four priorities, there are opportunities for joint production since we can achieve a 13.5% attainment for developed areas exposure, 7.5% attainment for suppression difficulty, 3% attainment for ecosystem services, and a 2.5% attainment for forest production. The selected projects are very scattered across the study area, and we do not notice the clusters formed by the projects of the other four priorities, with the highest ranked projects located in the prefectures of Kavala, Serres, Chalkidiki, and Thessaloniki.

Scenario 2: Annual Treatments
We found that for the protection of developed areas priority, approximately 15% of the total attainment can be achieved in four years with 34 projects treating 2195 ha ( Figure  A1 and Table A1 in Appendix B). For the forest production priority, a 5% attainment can be achieved from 289 projects treating 27,144 ha in four years ( Figure A2). For the ecosystem services priority, a potential attainment of 4% can be achieved with 281 projects applied on 27,250 ha in four years ( Figure A3). To achieve a 6% attainment for the suppression difficulty priority, 181 projects were needed treating approximately 17,000 ha in four years ( Figure A4). Finally, an attainment of 16% is possible for the wildfire transmission to protected areas priority if treatments were applied to 70 projects on 6500 ha in four years ( Figure A5). For all the five priorities, the number of projects implemented increased over time to achieve the same level of attainment compared to year one.
When the abovementioned projects were combined to remove duplicate stands (i.e., stands that were selected by the prioritization process for more than one priority), the total area treated during the four-year treatment plan was 77,350 ha with 855 projects on lands with sufficient growing stock volume ( Figure 6). Most of these projects are in the mountainous zones of the prefectures of Imathia and Pieria, in Kastoria and Chalkidiki. If the entire area shown in red in Figure 6 was treated, we expect that an attainment of 24% can be achieved for the protection of developed areas and wildfire transmission to protected areas priorities, 11% for the suppression difficulty priority, 9% for the forest production priority, and 8.5% for the ecosystem services priority. Compared to the attainment percentages achieved from the four-year management plan of each priority

Scenario 2: Annual Treatments
We found that for the protection of developed areas priority, approximately 15% of the total attainment can be achieved in four years with 34 projects treating 2195 ha ( Figure A1 and Table A1 in Appendix B). For the forest production priority, a 5% attainment can be achieved from 289 projects treating 27,144 ha in four years ( Figure A2). For the ecosystem services priority, a potential attainment of 4% can be achieved with 281 projects applied on 27,250 ha in four years ( Figure A3). To achieve a 6% attainment for the suppression difficulty priority, 181 projects were needed treating approximately 17,000 ha in four years ( Figure A4). Finally, an attainment of 16% is possible for the wildfire transmission to protected areas priority if treatments were applied to 70 projects on 6500 ha in four years ( Figure A5). For all the five priorities, the number of projects implemented increased over time to achieve the same level of attainment compared to year one.
When the abovementioned projects were combined to remove duplicate stands (i.e., stands that were selected by the prioritization process for more than one priority), the total area treated during the four-year treatment plan was 77,350 ha with 855 projects on lands with sufficient growing stock volume ( Figure 6). Most of these projects are in the mountainous zones of the prefectures of Imathia and Pieria, in Kastoria and Chalkidiki. If the entire area shown in red in Figure 6 was treated, we expect that an attainment of 24% can be achieved for the protection of developed areas and wildfire transmission to protected areas priorities, 11% for the suppression difficulty priority, 9% for the forest production priority, and 8.5% for the ecosystem services priority. Compared to the attainment percentages achieved from the four-year management plan of each priority (see previous subsection), the attainment of each priority from the combined projects ( Figure 6) was higher due to the larger area treated (77,350 ha). (see previous subsection), the attainment of each priority from the combined projects ( Figure 6) was higher due to the larger area treated (77,350 ha). Figure 6. The best 100-hectare projects for each priority with a four-year fuel treatment plan with an annual treatment pace, combined to remove repeated inclusion of the same stands selected by two or more priorities. In total, 855 projects resulted in 77,350 ha of treated area. Attainment for each priority, if all stands that appear on the map are treated, is shown at the upper left corner.

Scenario 3: Tradeoffs
In Figure 7 we mapped the 15 forest districts with the best joint production for each pair of priorities. Each point located along each line represents a 250-hectare project that resulted from the ForSys simulations with varying weight combinations between the two priorities for a particular forest district. It should be noted that the forest districts not shown could have any shape-some with straight diagonal lines indicate that increasing the attainment of one priority results in a proportional increase in the attainment of the other priority. Figure 6. The best 100-hectare projects for each priority with a four-year fuel treatment plan with an annual treatment pace, combined to remove repeated inclusion of the same stands selected by two or more priorities. In total, 855 projects resulted in 77,350 ha of treated area. Attainment for each priority, if all stands that appear on the map are treated, is shown at the upper left corner.

Scenario 3: Tradeoffs
In Figure 7 we mapped the 15 forest districts with the best joint production for each pair of priorities. Each point located along each line represents a 250-hectare project that resulted from the ForSys simulations with varying weight combinations between the two priorities for a particular forest district. It should be noted that the forest districts not shown could have any shape-some with straight diagonal lines indicate that increasing the attainment of one priority results in a proportional increase in the attainment of the other priority.
The sharpest tradeoffs with concave lines were produced for the forest production and suppression difficulty priorities ( Figure 7A), with elongated curve-like shapes starting high on the y-axis (suppression difficulty priority) and ending with low suppression difficulty but with very high forest production attainment percentages. Most of these forest districts are in Chalkidiki and in the eastern parts of the study area ( Figure 7A). Forest district 58 is on the frontier for both priorities, and yet the district is tiny, so for very little area treated we can get high attainment. The relationship between forest production and wildfire transmission to protected areas produced a variety of line shapes ( Figure 7B), i.e., two are represented with linear shapes and four with convex shapes, while the remaining are concave but with less sharp tradeoffs (shorter lines). Four forest districts are in the western parts of the study area, but most are in the prefectures of Serres and Kavala ( Figure 7B). The lines between forest production and protection of developed area are mostly concave, indicating sharp tradeoffs ( Figure 7C) on the forest districts of the central and eastern parts of the study area (prefectures of Thessaloniki, Serres, and Kilkis). Some forest districts, such as districts 456, 36, 5 and 1, appear in all four pairs of priorities, and some such as districts 456 and 1 are on the frontier for more than one pair of priorities. Finally, regarding forest production and ecosystem services, except for one forest district in Chalkidiki, all the forest districts have flat lines, indicating that joint production of the two priorities was limited since ecosystem service attainment was constant with few opportunities to increase attainment regardless of forest production ( Figure 7D). The sharpest tradeoffs with concave lines were produced for the forest production and suppression difficulty priorities ( Figure 7A), with elongated curve-like shapes

Discussion
This work is the first application of large-scale fire simulation modelling to evaluate forest management planning for wildfire exposure mitigation in Greece, and it is the first to introduce a prioritization scheme for forest management priorities through scenario planning. Specifically, it provides measurable outcomes on how focusing on one forest management priority results in tradeoffs in the others, and it identifies fuel treatment projects to achieve multiple fire management outcomes [6,[53][54][55]. Recent studies in the United States showed that treating less than 30% of landscapes can have substantial effects on wildfire spread and intensity; in other words, not every stand needs to be treated [29].
Overall, across all scenarios, the attainment of ecosystem services and forest production were correlated, suggesting that the same projects have high values for both priorities (i.e., in places where a large amount of growing stock volume is obtained, the value of land is also high). This is not surprising given the land value dataset ( Figure 2D) includes timber value in the estimation process. For Scenario 1, more efficient forest management planning is possible if stands with high attainment for two or more priorities are selected, minimizing the treated area and maximizing the expected positive outcome. Across the landscape there are more opportunities to treat the sources of exposure to developed areas and protected areas than other priorities. Protection of developed area attainment can be very high by treating fewer hectares in projects closer to urban areas, where simulation results predicted high structure exposure. As a result, with an overall low investment, a big change can be achieved in reducing fire exposure to developed areas. The prioritization of protection of developed area projects can also achieve high attainment for the wildfire transmission to protected areas and for suppression difficulty, however comes at a cost to forest production and ecosystem services attainment. Prioritizing forest production sacrifices attainment of the other priorities that are less than 1% of the total attainment, except for ecosystem services. Lastly, when we compare the modelled area selected for treatments ( Figure 3A) with the reported annual treated area by prefecture (Table A2), we notice that Thessaloniki and Kavala-although both with historically low area treated-were selected for a large number of modelled projects. This is probably a result of the projects that are not timber retrieval oriented, such as developed areas protection or suppression difficulty.
For Scenario 2, looking for stands that address all five priorities allows for the treatment of areas that are sources of fire to developed or protected areas-and to a lesser extent for ecosystem services, forest production or suppression difficulty. This is indicative of the cost to prioritizing developed or protected areas. Compared to the reported annual treated area (Table A2), results revealed that 6 out of the 13 prefectures that historically received high management are also those with the highest modelled treatments, although ranked differently. Our hypothetical fuel treatment plan for Scenario 2 was designed to be implemented within four years. Treating projects with the four-year management plan can be effective to rapidly reach some goals of attainment for each priority. Focusing on the first two years, fewer projects with less area are required for management to achieve the same attainment percentage compared to the following two years. Joining all projects for the five priorities (see Figure 6) proved ineffective compared to the all lands results from Scenario 1 (see Figure 3A), since we almost doubled the combined total area treated for all priorities but this did not result in doubling the attainment achieved. Overall, both scenarios suggest that it is better to start the application of treatments on the best projects up to a certain threshold at which attainment ceases to increase substantially, rather than treating all modelled projects (low return of investment).
Scenario 3 highlights the tradeoffs between forest production and the other three priorities inside forest districts. Sharp project tradeoffs between forest production and suppression difficulty were revealed. This also shows the forest districts where there are joint opportunities for attaining timber and another priority, with most of them located in the northeast. However, we must note that this scenario in not a pure tradeoff analysis since we applied a timber volume threshold, which removed stands that could achieve very high attainment, e.g., for the protection of developed lands but had low timber values.
Recent studies emphasized the need to incorporate wildfire uncertainty into effective fuel reduction plans, prompting discussions on whether future wildfires burn the project areas prioritized for fuel treatments prior to their scheduled implementation, and whether this would be executed at a fast pace [29]. For example, despite the widespread implementation of fuel treatments in the United States, their benefits are diminished over time on large landscapes due to the low probability that treated areas will be burned by a subsequent fire within a treatment's lifespan (temporal mismatch) [56]. Our prioritization process results indicate how much area should be treated and where, but to achieve the estimated levels of attainment, forest management should be effectively applied with vegetation and fuel treatment methods aligned with international experience and practice, while considering the abovementioned temporal mismatch [57,58]. In brief, fuel treatment strategies should employ a combination of surface fuel loading, depth, and continuity reduction treatments (e.g., prescribed burns and mastication), silvicultural practices to change tree crown structure (e.g., thinning and low pruning), and the creation of infrastructure and safety areas to facilitate fire suppression activities (e.g., road networks and water points). Mechanical treatments such as thinning (either commercial or noncommercial depending on stand conditions), mastication, or entire tree harvesting are required in high fuel load conditions or dense forest ecosystems with ladder fuels to reduce canopy bulk density and to mitigate hazard prior to using fire to reduce fuels. Most recent studies agree that thin and burn treatments had positive effects in terms of reducing fire severity, tree mortality, and crown scorch, but burning or thinning alone had either less of an effect or none at all, compared to untreated sites [59].
The temporal mismatch and the appropriate fuel treatment prescriptions for each area should be considered when fire management agencies want to apply fire risk reduction fuel treatments. Significant administrative and application constraints need to be addressed in order to implement the simulated scenarios and increase the chances for reducing fire behavior. In our previous research [38,60], we highlighted that prescribed burning is illegal in Greece and that a policy reform is required to allow for its application. As a result, the investment cost to achieve the expected outcomes is significantly higher compared to a similar application in other countries with a well-established fire use culture. Moreover, the major logistical constraint is how downed logs are moved to the closest timber mill and what infrastructure is required (e.g., opening of new roads). Another consideration is the management of logging residuals to avoid increasing fire risk in cases where no pile burning or mastication follows. Timber harvest volume with silvicultural treatments is determined in the forest management plan with an allowed variation of 10% in the volume of non-firewood products and 20% in the volume of firewood (Legislative Decree No. 86/1969), which is another legislative limitation towards increasing the pace of fuel treatments, especially for forest districts with recently established forest management plans. For reference, in Table A3 we show the reported annual timber harvest volume achieved in each prefecture of the study area for the period 2011-2020.
Once the projects and the stands for fuel management are selected on the larger scale that this analysis was conducted, a smaller scale approach is then required to allocate and arrange treatments on the landscape to maximize their efficiency. Previous studies have demonstrated how this can be applied in Greece [20] and elsewhere [52,[61][62][63][64]. For the case of stands with high attainment for the suppression difficulty priority, our results can be used to design fuel breaks in the form of a network with sufficiently modified fuels that can decrease the potential fire intensity and the rate of spread to a level where suppression resources succeed in containing the fire [64,65].
We identified some limitations in our assumptions and modelling. The ForSys model is data driven, and the inputs used for both fire simulation modelling and scenario planning were a snapshot of the most recent forest and fuel conditions, confined by data availability that is limited for Greece. Results should be interpreted and applied within a short-term timeframe (1 to 4 years). Moreover, management activities were simulated without the effect of potential disturbances that can occur within the application timeframe. As a result, forest growth was not considered in our estimations, and results are applicable until substantial forest growth occurs or until the effect of a disturbance substantially changes stand structure or fuel conditions. This was another reason why a four-year fuel management plan in Scenario 2 was simulated.
The major assumptions we made in scenario planning modelling included the following: (1) the stands selected for forest management would receive an appropriate suite of treatments, including mechanical thinning, fuels mastication, and pile burning of logging residuals; and (2) the posttreatment conditions of the stand would meet priority goals for the various objectives, i.e., reduce transmission of wildfires to developed areas or timber retrieval goals. In addition, our approach for assessing the levels of attainment from applying forest management in the selected project areas are relative rather than absolute values. In terms of forest production, since recent estimates of merchantable timber volume were not available, we used growing stock volume which is a pretty good analog to merchantable timber volume. Some of the selected stands within projects might not be accessible due to lack of roads, which implies the cost of treatments there would be increased.
ForSys has the capacity to solve large landscape problems (>10 7 ha), thereby overcoming computational limitations that existing forest landscape simulation models have [66,67]. On the other hand, ForSys does not have the functionality of a landscape disturbance succession model [66]. ForSys also does not require spatial solver software such as Google CPLEX. Compared to other spatial optimization approaches [68], our approach to applying scenario planning with stochastic fire simulations and the use of ForSys is based on specific objectives and not just criteria derived from expert knowledge; in other words, we decide the objectives to be tested, we set the management constraints (budget, area, thresholds), and then we produce the spatial optimization. Another major difference is that scenario planning with ForSys provides the outline of candidate areas for receiving fuel treatment/forest management projects and does not show how the projects should be translated into fuel management units and their allocation. Recent studies in Northern Mediterranean Sea showcased how localized fuel treatment allocation and assessment of their performance in mitigating fire behavior can be planned [20,64,[69][70][71].
Finally, since increasing fuel loads and continuity represent the main factor responsible for the recent catastrophic wind-driven fire events in Greece, we anticipate that our results can inform fire risk management agencies on the optimization of forest management in order to improve fire suppression efficiency and reduce fire spread to the wildland-urban interface and protected areas. Managing forest fuels has been internationally demonstrated to be an efficient strategy to fragment fuels and reduce fire spread rates and severity. We advise the Greek state to promote legislation that requires new forest management plans to include analyses examining the operational, economic, and logistic aspects of implementing different forest management scenarios that target meeting different priorities (multipurpose forestry), with measurable outcomes not only for timber harvesting targets but for other important priorities as well, such as community protection from wildfires.

Conclusions
Our results revealed the tradeoffs and opportunities among different forest management priorities, and our work produced treatment mosaics that optimize the delineation of forest management units. The applied methods provide the potential to contribute to improving the efficiency of wildfire management investments aimed at creating fire resilient ecosystems, facilitating safe and efficient fire suppression, and safeguarding rural communities from catastrophic wildfires. Our future research is expected to focus on estimating how different allocations of fuel treatments inside the proposed projects can affect fire behavior using simulation modelling. In addition, we intend to assess the efficiency of current forest management plans and fuel treatments outside managed forests by computing whether they spatially coincide with the proposed project locations. In addition, we plan to incorporate economic modelling with scenario planning to include timber haul distances and costs to the closest mill locations, thus allowing us to calculate harvest costs by treatment type and area as well as retrieved timber commercial value.
Funding: This research is co-financed by Greece and the European Union (European Social Fund; ESF) through the Operational Programme "Human Resources Development, Education and Lifelong Learning 2014-2020" in the context of the project "Spatial optimization of alternative forest fuel management scenarios for wildfire risk assessment" (MIS: 5048196). The APC was funded by the USDA Forest Service.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. Figure A1. Annual fuel treatment locations of the best 100-hectare projects for a four-year management plan, as estimated from the prioritization of protection of developed areas. Figure A2. Annual fuel treatment locations of the best 100-hectare projects for a four-year management plan, as estimated from the prioritization of forest production. Figure A2. Annual fuel treatment locations of the best 100-hectare projects for a four-year management plan, as estimated from the prioritization of forest production.  Figure A3. Annual fuel treatment locations of the best 100-hectare projects for a four-year period, as estimated from the prioritization of ecosystem services. Figure A3. Annual fuel treatment locations of the best 100-hectare projects for a four-year period, as estimated from the prioritization of ecosystem services. Figure A3. Annual fuel treatment locations of the best 100-hectare projects for a four-year period, as estimated from the prioritization of ecosystem services. Figure A4. Annual fuel treatment locations of the best 100-hectare projects for a four-year period, as estimated from the prioritization of suppression difficulty. Figure A4. Annual fuel treatment locations of the best 100-hectare projects for a four-year period, as estimated from the prioritization of suppression difficulty.
Forests 2021, 12, x FOR PEER REVIEW 26 of 31 Figure A5. Annual fuel treatment locations of the best 100-hectare projects for a four-year management plan, as estimated from the prioritization of wildfire transmission to protected areas. Figure A5. Annual fuel treatment locations of the best 100-hectare projects for a four-year management plan, as estimated from the prioritization of wildfire transmission to protected areas.    (1) Ranked projects from lowest to highest attainment and mapped them (see Figure 2B-F).
(2) Identified which stands were selected within each priority and examined how many priorities were met by each stand (see Figure 3A); for all stands with ≥2 priorities, we estimated the cumulative attainment for each priority. (1) Find projects until one-third of the attainment achieved by all simulated projects is treated; equal attainment treated for each of the four years; projects can be selected in any forest district, i.e., they do not need to be equally distributed across districts ( Figures A1-A5). (2) Identified which stands were selected within each priority and after removing duplicates, and examined attainment achieved for each priority if all stands were treated (see Figure 6). Selected the best 15 forest districts for each pairwise comparison; showed tradeoffs at forest district level; all points on each curve represent projects (see Figure 7). Attainment was estimated at project level. P1: up to 3% P2: up to 0.25% P3: up to 0.35% P4: up to 0.4% P5: up to 3.5% n/a P1: protection of developed areas; P2: forest production; P3: ecosystem services; P4: suppression difficulty; P5: wildfire transmission to protected areas.