Optimizing Fuel Treatments Allocation to Protect the Wildland–Urban Interface from Large-Scale Wildfires in Greece

Abstract

A crucial risk governance priority of the Greek forest managers is to reduce damages in the wildland–urban interface (WUI) by controlling wildfire behavior through fuel management practices. To support decisions for where management should be applied and how, this study experimented with new methods for fuel treatments allocation over a typical Mediterranean fire-prone landscape in the peninsula of Kassandra (an area of 350 km²), northern Greece. The Minimum Travel Time (MTT) fire simulation algorithm and the Treatment Optimization Model were used to produce eight spatial exclusionary and non-exclusionary datasets that were used as criteria for the spatial optimization of fuel management interventions. We used the Multicriteria Decisions Analysis method with Geographical Information Systems to cartographically intersect the criteria to produce two priority maps for two forest management scenarios (i.e., a control and a realistic one). The results revealed that 48 km² of the study area was characterized as high-priority locations in the control scenario (i.e., with equally weighted management priorities), while 60 km² was assigned to the high-priority class in the realistic scenario (i.e., with different weighted management priorities). Further analysis showed a substantial variation in treatment priority among the four major forest land cover types (broadleaves, sparse Mediterranean shrublands, conifers, and dense Mediterranean shrublands), revealing that the latter two had the highest selection values. Our methodological framework has already been operationally used by the Greek Forest Service branch of Kassandra to decide the most effective landscape fuel treatment allocation.

1. Introduction

The increase in continuity and accumulation of forest fuels is attributed mostly to the prevailing policy of abandoning forest management, especially in low-elevation Mediterranean conifer and maquis stands, to save money in favor of fire suppression. That leaves a limited budget available for fire prevention efforts [1,2,3,4]. Moreover, residential and forest flammable vegetation mixing zones, known as the wildland–urban interface (WUI), have expanded rapidly in recent decades, resulting in a dramatic increase in the number of people and structures affected by wildfires worldwide [5,6,7]. Understanding where the highest WUI risk emanates and what measures are necessary to make the WUI less exposed and vulnerable is an important challenge for natural resources managers. This is a very well-understood concept in the US, where wildland fuel treatments to reduce wildfire risk in adjacent communities and houses are a high priority for federal land management agencies, such as the US Forest Service [8].

Fuel treatments are forestry measures that change surface and canopy fuel characteristics. They aim at mitigating fire behavior, thus increasing fire suppression effectiveness. For each case, the most appropriate fuel treatment technique depends on factors such as the type and size of the forest stands, fuel characteristics, topography, available funding, land ownership, relevant management experience, traditions, and legislation [9]. According to the Fire Protection Plans for Forests and Forest Ecosystems specifications [10], the term “cleanings” is used to describe a forest cultivation practice to achieve disruption of forest fuel connectivity. By applying this method, the shrub understory is partially removed and the trees are pruned at a height ranging from one third to half of their total height. A minimum distance of 2 m from the ground is retained to create a vertical fire-break zone between the surface and the forest canopy. The dominant, stunted, dry trees that are affected by insects and fungi are logged, and the canopy cover is maintained such that growth of the shrub understory is not favored. The retrieved large-diameter biomass is attributed to cover the needs of local residents for heating, while small-diameter biomass is chipped and scattered in the forest floor [11]. When the treatments are made in strips or zones, their width is defined by minimum ecological and fire safety boundaries. Where treatments are applied in the whole stand or area, as general forestry cultivation or management controls, they are regulated by relevant laws to ensure that economic or other management purposes (hunting, soil protection, etc.) are achieved.

Previous studies have widely analyzed the complexity of designing fuel treatment plans for changing forest structure and reducing wildfire behavior [12,13,14,15,16,17,18]. Finney et al. [19] developed an optimization algorithm identifying fuel management priority areas that can achieve the largest delay in fire spread rates. In other studies, treatment strategies were evaluated by fire spread simulators, followed by the calculation of fire exposure parameters that were subsequently used in tradeoff analyses for large-scale restoration programs [20,21,22,23,24,25,26]. This approach has been successfully implemented in the Mediterranean Basin over the last few years. In Sardinia, Italy, Salis et al. [27] defined three fuel treatment strategies and objectives, and analyzed how these different strategies affected wildfire exposure. In Catalonia, Spain, Alcasena et al. [2] identified suitable strategic locations on forest lands for prescribed fire treatments to achieve three competing objectives, and calculated production possibility frontiers (PPFs) to explore the trade-offs among the different objectives. In Macedonia, Greece, Palaiologou et al. [28] developed three alternative forest management scenarios to evaluate the effectiveness of five management priorities.

The need for scientifically based decisions, especially in cases of complex problems that require the study of many conflicting factors, was highlighted not only within the framework of wildfires but also regarding other socioecological challenges [29,30,31]. Thus, Multicriteria Decisions Analysis (MCDA) is a typical procedure that provides the ability to identify areas based on a variety of attributes and criteria each area should meet [32]. When the Geographical Information System (GIS) and MCDA are combined, a significant relevance multicriteria spatial analysis is achieved, known as GIS-based multicriteria decision analysis (GIS-MCDA) [33,34,35,36].

In this study, we applied the GIS-MCDA method for fuel treatment spatial optimization over a typical Mediterranean wildfire-prone landscape. Our goal was to assist the responsible managers of the local Forest Service branch, which has the main responsibility for forest fire prevention, in their decisions to prioritize forested areas for receiving fuel treatments targeting the protection of human life and property. To our knowledge, this work is the first application of fuel treatment spatial optimization in Greece that combines fire simulation modeling, GIS-MCDA, factor weights, and forest management methods, and one of the first in the Mediterranean Basin.

Results from the analysis presented in this paper on spatial optimization of the areas that must be managed were used to draw up the local forest management plan for the following years. The plan has already been funded and applied in the Kassandra peninsula (Chalkidiki, northern Greece), through the implementation of the ”AntiNERO” fire prevention management program of Greece during 2022. Since this is a three year forest management plan, its effectiveness can only be assessed and recorded when a few years pass, and if these treatments experience a wildfire event during their lifetime. Interesting findings were acquired from the first-year verification process, discussed in the following chapters of this work.

2. Materials and Methods

2.1. Study Area Description

Kassandra covers an area of 350 km² and is the western peninsula of Chalkidiki Prefecture, northern Greece, located about 60 km south of Thessaloniki (Figure 1a). The topography is relatively smooth, apart from some isolated northwest-facing slopes that are steeper. The maximum altitude is approximately 320 m above sea level. Its climate is characterized as Mediterranean type-temperate, with mild winters and dry hot summers. The mean annual precipitation is approximately 580 mm, and the mean temperature is ~17 °C, with the dry period lasting from May to September with a mean maximum temperature of ~30 °C.

Bachantourian et al. [37] estimated the land cover of the Kassandra peninsula using ground-survey acquired data, mostly dominated by Aleppo pine (Pinus halepensis Mill.) (44%) and agricultural areas (44.5%), comprised of olive trees (Olea europaea L.) and non-irrigated croplands that occupy a unique large continuous feature area in the north and the central part of the peninsula (Figure 1b, in yellow). Moreover, approximately 10% is urban cores and WUI areas with hundreds of dispersed recreational houses and hotels intermingled with the wildland forest vegetation (Figure 1b). The natural attractions of the study area, a combination of pine forests and sandy beaches, led to land use changes in favor of residential development for tourism and recreational uses during the last 40–50 years.

In the last 40 years, wildfires (>0.1 km²) burned 29% of the total area, while some regions were impacted by multiple events. Nowadays, most of these formerly forested areas have been degraded to other vegetation ecosystems as the dense and mature conifer forest has been replaced by a dense understory of evergreen sclerophyllous shrubs (Mediterranean maquis) consisting of Pistacia lentiscus L., Quercus coccifera L., and Arbutus unedo L., exceeding 2 m in height in some stands [37,38,39].

Until 2021, forest management plans in Kassandra were limited to logging permissions that the Greek Forest Service (GFS) issued to residents to cover their household needs, estimating that approximately 3500–4000 m³ of dry timber were removed annually from the forests by this procedure. Alternatively, in cases of exceptional and unexpected natural disasters, e.g., a strong summer windstorm in 2019 or the 2006 wildfire, “emergency harvesting” is applied and forest management permits are given to Forest Logging Cooperatives under the supervision of GFS [40]. However, the recent devastating fires motivated the Greek State to holistically restructure the policy of wildfire confrontation, adapting its fire management planning to more prevention-oriented approaches [28]. Thus, a forestry fire protection program called “AntiNERO” was funded with a budget of 72 million euros for prevention projects inside public forests across the country, from which 1.4 million euros were allocated in Kasssandra. The project was implemented from June to November of 2022, and 120 ha of public forests received fuel treatments, while 650 km of the forest road network and five water tanks for firefighting use were maintained.

2.2. Land Cover Groups, Classes, and Fuel Models

In our previous work [37], we developed a fine-resolution land cover dataset, hereafter referred to as the “Kassandra Dataset”, after combining ground surveys with orthophotos. This dataset was used to categorize land cover types into wildfire fuel models suitable for fire behavior modeling. Initially, the process analyzed vegetation data from field data collected on 154 circular sampling plots (0.1 ha each), allocated inside forested areas throughout the Kassandra peninsula. According to the principles and methods of forestry science [41], dendrometry characteristics were measured and estimated, i.e., stand height (SH), diameter at breast height (DBH), crown base height (CBH), canopy cover (CC), and minimum and maximum diameter below the crown; while a detailed inventory of the forest composition (i.e., dominant species, coverage percentages) and physical conditions (i.e., weather damages, diseases) of each sampling plot were recorded. The final land cover type product came up as the combination of the ground surveys dataset with the official digital vector layer from the Chalkidiki Forestry Cadastre [42], literature data and reports, i.e., forestry management plans and fire protection plans, and photointerpretation of high-resolution (0.25 m × 0.25 m) Large Scale Orthophotography (LSO 25/2015) in natural color. For our research, we also used vector layers with the subsidized olive-producing areas [43], the dominant wildfires for the period 1980–2021, and the forest stands [37] to create fuel model maps. Moreover, the discontinuous urban developed areas, i.e., villages, dispersed country houses, and hotel units [44], and Kassandra’s road network [45] were used to calculate proximity.

By considering the predicted fire behavior that each one of the 12 land cover classes of the Kassandra Dataset can potentially produce, we assigned for each land cover type one or more fuel models from the 40 standard fuel models of Scott and Burgan [46]. The fuel models represent sets of parameters that define dead and live fuel properties to various fire behavior modeling algorithms [47]. The selected fuel models for each land cover type were then inserted in the Behave Plus 6 software [48] to test which of their simulated fire characteristics were closer to the behavior of a real fire. As the fuel model map was created before the fuel treatments that were implemented during the forestry fire protection program AntiNERO 2022, the fuel type changes were not included in this map.

The complete list of the Kassandra Dataset with five class groups, 12 land cover classes, and 11 fuel models, one assigned to each class, and their corresponding percentages across the study region are shown in

Table 1. Kassandra Dataset land cover groups and classes, matched with the appropriate fuel model and their corresponding percentages. LC: Land Cover.

LC Group	LC Classes	Fuel Model Code	Fuel Model Name	LC Class PERCENTAGE	Km2	LC Group Percentage
Artificial surfaces	Urban/Suburban developed	91	NB1	2.91	10.19	9.61
	WUI	102	GR2	6.70	23.51
Agricultural areas	Croplands	101	GR1	27.65	96.92	44.43
	Olive groves	161	TU1	16.78	58.83
Forest	Shrublands/moderate load	142	SH2	0.64	2.24	43.96
	Sclerophyllous vegetation/Mediterranean maquis	147	SH7	8.13	28.48
	Coniferous forest/Treated	161	TU1	1.84	6.47
	Coniferous forest/Dwarf conifer	164	TU4	9.10	31.90
	Coniferous forest	165	TU5	23.96	83.99
	Broad-leaved forest	182	TL2	0.29	1.00
Wetlands	Open water	98	NB8	0.66	2.32	0.66
Barren	Bare ground	99	NB9	1.34	4.69	1.34
	Total			100.00	350.55	100.00

.

2.3. Input Data for Wildfire Simulations

Wildfire simulation modeling was conducted using the minimum travel time (MTT) fire spread algorithm and the treatment optimization model (TOM), as proposed by Finney [19], to obtain two fire-related metrics that can provide insights on where fuel management should be applied, i.e., the fire transmission to houses (FTH) to locate areas where ignitions can potentially affect a high number of houses and the fuel treatment grid (FTG) to identify areas where treatments allocated at the landscape level can most effectively retard fire growth rates.

A gridded landscape file (called LCP) with spatial raster data and a high special resolution (30 m × 30 m) was created to deliver the required spatial inputs in the wildfire simulation model of FlamMap [49]. The elevation, aspect, and slope required layers were retrieved from the digital elevation model (DEM) of high-resolution LSO 25/2015 [50], while values for the stand height (SH), the crown base height (CBH), and the canopy cover (CC) were retrieved from the Kassandra Dataset. Canopy bulk density for Aleppo pine was first estimated at the tree level, followed by a reduction to stand level [51]. The total weight of the available crown biomass (Crown Fuel Weight) was calculated using the allometric equations of Mitsopoulos and Dimitrakopoulos for the Aleppo pine [52]. Surface fuel quantity and arrangement were described with fuel models (see Section 2.2).

We used meteorological data from the prevailing weather conditions during a large wildfire occurrence in the study area on 21 August 2006, regarding the air temperature, the relative air humidity, the respective wind speed at 6.1 m high, and the wind direction. Dead fuel moisture and foliar moisture content were set as corresponding to very dry canopy fuel conditions, using the BehavePlus D1L1 scenario (very low dead, fully cured herb) (

Table 2. Fire weather data corresponding to the prevailing weather conditions during the largest wildfire event in the study area (21 August 2006), and other parameters used for wildfire simulations.

Parameter	Value	Units
1 h fuel moisture (0–0.64 cm)	3	%
10 h fuel moisture (0.65–2.5 cm)	4	%
100 h fuel moisture (2.6–7.5 cm)	5	%
Live herbaceous fuel moisture (LH)	30	%
Live woody fuel moisture (LW)	60	%
Wind direction	90	(0)
Wind speed	62	Km × h−1
Foliar moisture content	70	%
Spot Probability	0.25
Resolution of calculation	30	m

). Regarding spot probability, as FlamMap only considers the volume but not the type (i.e., tree species) of canopy fuels, we used a universal value for the forest type of the study area that produces the most and more flammable embers, i.e., Pinus halepensis Mill.

A total of 10,000 randomly allocated ignitions across the study area were simulated with the ΜΤΤ algorithm of FlamMap 6.1 [53] to generate fire behavior metrics such as fire rate of spread, burn probability, and conditional flame length that consequently were included as GIS-MCDA criteria. This high number of simulations ensured that all pixels with burnable fuels were burned at least once. For the simulations with the use of TOM, equal parameterization values with MTT were used. TOM requires as input the ideal landscape that describes the proposed management implementation to each pixel of the landscape in which treatments can potentially be applied, and not only where they are planned to be applied. The TOM output treatment opportunity grid (TOG) compares the rate of fire spread (ROS) between the real and the ideal landscape. The results are rendered on a spatial raster file where each pixel receives values of +1, 0, and −1. The positive sign reflects areas where ROS decreased after the treatment implementation, zero reflects those with no difference, and the negative sign reflects the areas where ROS increased after the treatments. The TOM selected cells with a value of +1 for treatment. Cells with zero and negative values were not treated. The MTT calculates the fastest travel routes of the fire in the field while TOM proposes (without applying it) the special optimization for management between those areas that showed a decrease in the intensity of fire behavior after the treatment implementation. During the simulations, no suppression by firefighting forces (air or ground) was taken into consideration.

2.4. Suitability and the Analytical Hierarchical Process

Eight spatial and legislative criteria were defined and categorized into two types, i.e., constraints (exclusionary criteria) and factors (non-exclusionary criteria). The exclusionary criteria limit the analysis to specific geographical areas, excluding those that may be considered unsuitable for management. By using the Raster Calculator of QGIS, for each exclusionary criterion, a Suitability layer, i.e., a binary file, was created to divide the study area into two classes. The class with value = 0 represents areas that are excluded from management, while a value = 1 represents areas included in the spatial optimization. The mathematical formulation for selecting areas using only exclusionary criteria is:

$$S_e={\displaystyle \prod }_{j=1}^K{b}_j$$

(1)

where S_e = exclusionary criteria suitability (0 or 1); b_j = criteria score of constraint j; K = number of constraining criteria; Π = product.

The non-exclusionary criteria that determined the degree of suitability of the various geographical areas were scored according to their distance from the point of interest and they were qualified on a 0–100 value scale, with zero representing areas with the lowest priority for management and 100 the areas with highest priority for management. When the factors are combined by applying a weight to each and a summation of the results is followed, the suitability can be estimated using the following formula [54]:

$$S_{ne}={\displaystyle \sum }_{i=0}^N{w}_{i }{x}_{i }$$

(2)

where: S_ne = non-exclusionary criteria suitability; w_i = weight of factor i; x_i = criterion score of factor i; N = number of factors.

For the cases where both types of criteria apply, the mathematical formulation for assigning the overall suitability (S_all) is:

$$S_{all}={\displaystyle \sum }_{i=0}^N{w}_{i }{x}_{i }\times {\displaystyle \prod }_{j=1}^K{b}_j$$

(3)

The spatial optimization of the management interventions was carried out by cartographic intersection of the data layers (criteria), converted to raster and resampled to the same spatial resolution (30 m × 30 m). The involved criteria were homogenized and reclassified from their various value scales to a 0-100 scale. According to their priorities/objective’s relative importance using the Analytic Hierarchy Process (AHP) method [32], weights were assigned to four of those criteria (see Section 2.6). AHP was proposed by Satty [55] and is based on the principle that during decision making, the experience and knowledge of people are equally important to the available data. Combining the products of spatial constraints with those of simulations by multi-objective analysis and cartographic overlaying, the most suitable areas for management were identified, and the treatment allocation maps were produced. The research was completed in the field where the spatial optimization maps were reviewed by more than 70 ground surveys and recordings. Thus, the treatment units selected in the office were confirmed or rejected in the field before finalizing the management plans.

2.5. Criteria Description and Application

Eight criteria were defined and used, based on expert choice, the relevant literature, or because they are enforced by the Greek and European legislation. The factors cover climatological, technical, location, cost, ecological, and social dimensions that impact site selection for receiving fuel treatments. Four of these were of constraint type (exclusionary criteria): (1) slope (SLP), (2) land ownership (OWN), (3) protected natural areas (PNA), and (4) fuel treatment grid (FTG); while four were non-exclusionary criteria: (5) suitable fuel models (SFMs), (6) proximity to the road network (PRN), (7) proximity to urban development areas (PUD), and (8) fire transmission to houses (FTH).

2.5.1. Slope of the Land Surface

Surface slope is a crucial factor that determines which areas are appropriate for fuel management, as steep slopes are naturally prohibitive for applying fuel treatments. Slope was derived in percent from the digital elevation model (DEM), as resulting from the LSO25/2015. The generated raster was resampled from the 0.25 cell size to the spatial resolution of 30 m, to be compatible with the other thematic criteria. A binary file was used, setting the limits of safe fuel treatment application. Based on forest managers’ experience, we indicated with value = 0 (i.e., exclusion) the areas with surface slope > 100%, while with value = 1 those areas with surface slope < 100% (Figure A1a).

2.5.2. Ownership of the Land

As our scope was to allocate treatments only in the forested lands, we used Kassandra’s dataset land cover to flag and remove artificial surfaces, agriculture areas, wetlands, and barren classes (value = 0), while the polygons of the forest land cover group (Table 1) received a value = 1, i.e., qualified for spatial optimization. The Greek legislation allows the Greek Forest Service to apply fuel treatment plans exclusively in public forests, so different forest ownerships, i.e., private, managed by the local government, and forests belonging to the Church, were excluded from receiving fuel treatments (Figure A1b).

2.5.3. Protected Natural Areas

In the northern part of the peninsula, an area of wetlands and forests that covers 40 ha is included in the category “Special Protection Areas” (SPA) and belongs to the European Natura 2000 Network with the code GR1270013 (H.P. 37338/1807/E.103/6.9.2010, Government Gazette 1495B′). According to the Council Directive 92/43/ECC of 21 May 1992 on habitats (European Union Law), approved management plans in Natura 2000 areas presuppose that the interventions do not have negative impacts on the ecological integrity of the site. In many cases, the management of the protected areas requires relevant resolutions by the Greek Ministry of the Environment with implementations of the corresponding compensatory measures. As the negative consequences of the interventions were not identified by this research, we excluded all Natura 2000 sites from receiving fuel treatments (Figure A1c).

2.5.4. Fuel Treatment Grid

The Fuel Treatment Grid is an output produced by the TOM of FlamMap 6.1 that identifies areas with priority for fuel management. It indicates the areas with the highest delay in fire spread after the fuel treatments implementation [19]. The ideal landscape was generated as a result of management without restrictions, through a set of rules that changed the values of fuel model (FM), canopy cover (CC), crown base height (CBH), and crown bulk density (CBD) (

Table 3. Set of rules applied to define how theoretical fuel treatments were applied for the generation of the ideal landscape.

Rule	Target Theme	Abbreviation	Rule Description
1	Crown Base Height	CBH	SET CBH × 2 IF FM = (164 OR 165) AND CBH ≥ 3 AND CBH ≤ 6
2	Crown Base Height	CBH	SET CBH × 3 IF FM = (164 OR 165) AND CBH ≤ 2.9
3	Crown Base Height	CBH	SET CBH × 2 IF FM = (161) AND CBD ≥ 0.06
4	Canopy Cover	CC	SET CC −20% IF CC ≥ 70%
5	Fuel Model	FM	SET FM = 141 IF FM = (142 OR 147)
6	Fuel Model	FM	SET FM = 161 IF FM = (164 OR 165)
7	Crow Bulk Density	CBD	CBD × 0.7 IF FM = (161 OR 164 OR 165 OR 182) AND CBD ≥ 0.12

).

Thus, dense conifer forests with understory (TU5 from Scott and Burgan fuel models) and dwarf conifers (TU4) in the current landscape were changed after the theoretical application of fuel treatments to low-loading conifer forests (TU1) in the ideal landscape, while moderate and heavy shrub fuel loading was changed to low shrub fuel loading (SH1). For dense conifer stands (CC ≥ 70%), canopy cover was reduced by 20%, CBH was increased by two times for trees with CBH between 3 and 6 m, while CBH was tripled for trees with CBH values ≤2.9 m. For the fuel model TU1, CBH was increased by two times for areas with CBD ≥ 0.06 kg/m³. With rule #3, the olive tree areas were separated from coniferous forests, as both were assigned with the fuel model TU1 and were excluded from theoretical fuel management in the ideal landscape. Finally, for conifers and broadleaved forests, CBD was 30% off for areas with CBD ≥ 0.12 kg/m³. The fuel treatment grid received the value of 1 for the areas proposed for treatment, while those that were rejected from the spatial optimization received the value of zero (Figure 2 and Figure A1d).

Regarding TOM execution, along the eastern coast of the peninsula, an ignition line perpendicular to the main wind direction was defined to affect the entire study area. The treatment fraction was set to 0.20, meaning that the algorithm would find 20% of the landscape that creates the highest difference in fire behavior. Fuel treatment areas were selected from those pixels that decrease fire behavior in the unmanaged landscape.

2.5.5. Suitable Fuel Models

Not all fuel models have the same priority for treatment, even though they might be suitable for receiving fuel treatments. Therefore, we reclassified the fuel model map into five classes to reflect fuel model suitability with a real number in the range from 0 to 100: “very low” (0), “low” (10), “moderate” (30), “high” (50), and “very high” (100). Areas with the fuel models TU5, TU4, and TL2 indicated very high priority for treatment and received a value of 100. Areas with the FM SH7 (shrublands—very high load) were classified in the high-priority class (value = 50), areas with FM TU1 (sparse or managed forest) received the moderate value (30), the FM SH2 (shrublands—medium load) received the low value (10), and the remaining fuel models that represent the non-forest lands were characterized as very low suitability (0) for fuel management (Figure A1e).

2.5.6. Proximity to the Road Network

Fuel management priority increases inversely proportional to the proximity to the road network, due to the high risk for fire initiation on either side of the roads and the ease of access to the sites, in implementing the treatments and retrieving and transporting the logging products. As a result, the fuel treatment implementation cost is lower closer to the road network. The road network was converted to a proximity raster layer that indicates the distance of the centroid of each pixel from its nearest feature (point of line vertex) using the Euclidean distance algorithm of QGIS, reclassified on five classes to reflect low to high suitability within the value range of 0-100. Thus, the “very high” class (value = 100) was attributed to pixels located within a distance of 50 m from the roads, the “high” class (value = 50) was assigned to pixels within 50–100 m, the “moderate” class (value = 30) was assigned for the pixels located within 100–500 m, and the “low” class (value = 10) was assigned for pixels at a distance >500 m from the roads. The pixels on the road edges were assigned with the value of zero (Figure A1f).

2.5.7. Proximity to Urban Development Areas

The fuel treatments close to residential and urban development areas are preferred when the fuel management goal is to protect communities and human-made assets. The two classes from the Kassandra dataset that represent artificial surfaces, i.e., Urban/Suburban developed and WUI areas, were joined to generate an ”urban development areas” layer that includes villages, dispersed country houses, hotel units, hotel resorts, and camping sites. The weighted average distance that each pixel had from a residential area was used for the treatment optimization, and pixels were classified into five classes: “very high” (value = 100), “high” (value = 50), “moderate” (value = 30), and “low” (value = 10), for pixels located within the distances of 500 m, 1000 m, 2000 m, and >2000 m from residential areas, respectively. The urban development areas received the value of zero and were excluded (Figure A1g).

2.5.8. Fire Transmission to Houses

The MTT was used to produce the non-exclusionary criteria of fire transmission to houses, which accounts for the number of houses exposed to fires ignited at a particular location, threatening even a single house [2,22,56]. The goal was to understand where dangerous fires initiate and, there, apply fuel treatments. The fire perimeters modeling grid layer output and the CSV file with the ignition points’ coordinates were joined and connected in QGIS [57]. The produced vector layer was intersected with the point data shapefile of the discontinuous urban developed area (a total of 25,632 records) [37], and the total number of the affected structures by each ignition point was estimated. By applying the Inverse Distance Weighting (IDW) interpolation process [58] on all simulated fires with exposed residential houses >0, a continuous surface with 30 m resolution was created (Figure 3), and the intermediate values for each pixel were inferred. Areas where fires started and were predicted to expose 1 to 50 structures were classified with “low” threat (value = 10). Areas that were predicted to expose 51 to 150 structures were classified with “moderate” threat (value = 30), those with a predicted exposure between 151 and 350 were classified with “high” threat (value = 50), and for predicted structure exposure > 350, the areas were classified with “very high” threat (value = 100) and received the maximum priority for management (Figure A1h).

2.6. Assigning Weights to Criteria

The overall score of each criterion is influenced by the weighting factor assigned to it and indicates the criterion’s importance; thus, it has a significant impact on the result. There are several methods for assigning factor weights, with the simplest being the application of identical coefficients in all criteria, dividing the unit (which is the sum of the coefficients) by the number of criteria [31]. However, AHP assigns different factor weights to the criteria, an approach that is considered more realistic. In this research, weights were assigned only to factors [59,60] based on the pairwise comparisons method. The resulting combinations depend on the criteria’s number and are given by the mathematical formulation:

$$\frac{n\left(n-1\right)}{2}=\frac{4\left(4-1\right)}{2}=6$$

(4)

The result of these comparisons is a pairwise comparison matrix (

Table 4. Pairwise comparisons of the importance of the four different criteria. Suitable Fuel Models (SFMs); Proximity to the Road Network (PRN); Proximity to Urban Development areas (PUD); Fire Transmission to Houses (FTH).

	SFM	PRN	PUD	FTH	Factor Weights
SFM	1.00	5.00	3.00	1.00	0.41
PRN	0.20	1.00	0.20	0.20	0.07
PUD	0.33	5.00	1.00	1.00	0.21
FTH	1.00	5.00	1.00	1.00	0.31

Consistency Ratio = 0.04.

), where the Fundamental Scale was used with the number 9 meaning absolute supremacy and the number 1 representing criteria with equal importance [55]. When the criterion of the line is more important than the comparable criterion of the column, it receives a number of 3, 5, 7, and 9 that yields in the corresponding order an increasing value. The comparison matrixes of the AHP’s process are available in Tables S1–S3 in Supplementary Materials S1. The factor weights were calculated by averaging each row of the normalized matrix, and an acceptable value for Consistency Ratio (CR) was produced. When CR < 0.10, the estimates are considered within normal limits of consistency as the array is better than a random array and it is, therefore, accepted; otherwise, subjective preferences should be reevaluated. CR is given by the mathematical formulation:

$$\mathrm{CR}=\frac{\mathrm{CI}}{\mathrm{RI}}$$

(5)

where CI = Consistency Index; RI= Random Consistency Index.

In this study, two scenarios were applied that differ in terms of the criteria’s weighting factors. In the first scenario, the weighting factors are the same (0.25) for all criteria (control scenario), whereas in the second scenario, as the landscape managers can have different priorities, different factor weights were assigned. We considered their contribution to the achievement of a defined goal based on the author’s personal experience regarding fuel treatments’ application and the significant ability to objectively assign weights to the criteria. Thus, applying the implemented technique of factor pairwise comparison, the highest weights were attributed to the SFM and FTH criteria, while the PUD and PNA followed with lower weights (

Table 5. Calculated factor weights for the two scenarios after applying the weighted sum method for the four different criteria. Suitable Fuel Models (SFMs); Proximity to the Road Network (PRN); Proximity to Urban Development areas (PUD); Fire Transmission to Houses (FTH).

	SFM	PRN	PUD	FTH
Scenario 1	0.25	0.25	0.25	0.25
Scenario 2	0.41	0.07	0.21	0.31

).

3. Results

3.1. Spatial Optimization Model

The output maps resulted from the application of GIS-MCDA on four of the criteria in a spatial resolution of 30 m × 30 m. The cartographic background was created by considering the exclusion criteria and removing all excluded pixels from receiving fuel treatments (Figure 4). An expected high percentage of the total area (80.94%), corresponding to 281 km², was excluded, while the remaining 66 km² (approximately 19%) qualified for fuel management. Subsequently, the appropriate weighting factor was assigned to the remaining four eligibility criteria, i.e., SFM, PNA, PUD, and FTH, that were consequently added together. The procedure was performed for both scenarios, resulting in two intermediate suitability maps, i.e., without excluding any pixel, whose results were divided into four classes, from low to very high suitability on a scale of 0–100; i.e., “low suitability” (0–10), “medium suitability” (11–30), “high suitability” (31–50), and “very high suitability” (51–100) (Figure 5a,b).

Table 6. The area of each suitability class and its proportion of the total study area for the two scenarios.

Suitability Class	Scenario 1		Scenario 2
	Area (km2)	Proportion %	Area (km2)	Proportion %
Low	4.64	1.33	13.57	3.90
Medium	61.62	17.71	85.44	24.56
High	141.08	40.55	95.34	27.40
Very high	140.62	40.41	153.46	44.14

shows the area for each group and its proportion of the total study area that resulted from each scenario. Comparison of the two scenarios revealed that the dominant differences were found for the class of “high suitability”. Thus, a significant decrease of 45.74 km² between the first and second scenario occurred for this class, which was evenly allocated and distributed in the other three suitability classes by increasing their values. The total fuel treatments’ suitability with a detailed description of the pixel sum, the area (km²) of each class, and the percentages for the two scenarios are reported in Tables S4 and S5 in Supplementary Materials S2.

The two final overall priority maps (Figure 6) resulted from the multiplication of the exclusion area map (Figure 4) with the corresponding intermediate suitability layer of each scenario (Figure 5), identifying areas that maximize the levels of achievement of multiple criteria to protect the residentially developed areas and receiving a range of values from 0 to 100. This rating classified the study area into four classes: “excluded for treatments” (value = 0), “low treatment priority” (values between 1 and 30), “medium treatment priority” (values between 31 and 50), and “high treatment priority” (values between 51 and 100) (Figure 6a,b).

The two overall priority plans have concentrated high-priority areas (in dark green) in central and eastern parts of the peninsula, where the main WUI areas are located. Approximately 19% of the study area can receive fuel treatments for each scenario as the percentage of the excluded areas remained unchanged in both scenarios (

Table 7. Fuel treatment spatial optimization results for the two scenarios.

Priority Class	Scenario 1		Scenario 2
	Area (km2)	Proportion %	Area (km2)	Proportion %
Excluded	281.36	80.94	281.36	80.94
Low	0.50	0.15	0.38	0.11
Medium	18.03	5.19	6.17	1.78
High	47.70	13.72	59.69	17.17

). A total of 47.70 km² (13.72%) in scenario 1 and 59.69 km² (17.17%) in scenario 2 are of high priority for receiving fuel treatments, revealing that the realistic scenario 2 increases the cost of the fuel treatment’s plan implementation. Medium and low classes have very low proportions of priority compared to the high class, proving the high risk of the majority of areas for major wildfire events. The detailed description of the pixel sum, the area (km²) of each class, and the percentages for the overall priority results for both scenarios are reported in Tables S6 and S7 in the Supplementary Materials S2.

In Figure 7, boxplots represent the average overall priority values indicating the different ranges and distribution between the two scenarios. The treatment unit’s average bulk of values is higher for Scenario 2, meaning that by using unequal weighting, more pixels are allocated in higher-priority classes. The two distributions are considered symmetrical as only one outlier (grey dot) is recorded below the 10th percent and above the 90th percent.

3.2. Fire Transmission to the Houses and Fuel Treatment Grid

The FTH showed interesting spatial patterns with a strong linkage to the major tourism areas that have high structure density. The results showed that some fire events could potentially affect up to 2500 structures. The higher values of FTH (>350 structures) observed at the northern part of Kassandra, i.e., Nea Potidea, and at the central part, i.e., Kalandra (Figure 3), occurred in crop agricultural lands that were excluded from receiving fuel treatments; thus, in those areas, we cannot do anything to lower structure exposure. However, areas covered by conifer forests such as Sani, Chanioti, and Pefkohori (Figure 3) have greater chances for the reduction in their fire exposure as they are located in high-fuel-treatment-priority areas (Figure 6). The lowest exposure (<50 structures) was found at the northeastern and central parts of the peninsula, in extended rural areas with low structure density. FTH is not necessarily related to fire size, as the highest values of FTH correspond to simulated fires between 400 and 600 ha, while the largest simulated fires (>3000 ha) expose less than 1350 structures (Figure 8).

Average FTH values inside the pixels selected for fuel treatments showed different distributions between the four major forest types. The values of FTH for all the four major forest types fell within the lower end of the range, and extreme values were recorded outside the whiskers, which means that the four distributions were not symmetrical (Figure 9a). Thus, most values were not close to the average value of the distribution, so the average number of structures exposed cannot be predicted considering the fire ignition’s land cover type. However, the four distributions differed significantly. A wildfire could potentially expose more structures when initiated from the SH7 fuel model, as this fuel type has the highest median and the widest range of values (50% of the observations affect between 148 and 438 houses). On the other hand, the narrowest value range was observed in broadleaved forests (FM TL2) and sparse shrublands (FM SH2). Coniferous land cover (with FM TU1, TU4, and TU5) had the lowest median value range but produced the most upper extreme values, although this is rather due to the variable environmental conditions encountered. These extreme upper values of the coniferous land cover distribution are mostly due to the expanded area (i.e., 122 km² vs. the other three fuel models where their sum area does not exceed 30 km²), and secondly due to its wide spatial distribution both near and far away from the WUI.

Green areas in Figure 2 showed TOM simulations results, indicating areas where fuel treatments can reduce wildfire behavior in the six fuel models (TU1, TU4, TU5, SH2, SH7, and TL2). Overall, the central and the eastern parts of Kassandra have the largest and most continuous areas for successful fuel treatments. The percentages of the FTG values showed substantial variation between the four forested land cover types (Figure 9b), i.e., 53.40% of the coniferous forests area can potentially reduce fire behavior after the implementation of fuel treatments, with the corresponding percentage for Mediterranean maquis as 42.56%. On the other hand, treatments allocated in sparse shrublands (FM SH2) have a very low potential for a successfully wildfire containment (3.17%), with the lowest values of FTG observed in broadleaves (FM TL2) with 0.26%, meaning that even if those areas were managed, they would not achieve a substantial reduction in fire behavior, due to their low flammability.

4. Discussion and Conclusions

This is the first application of small-scale fire simulation modeling to prioritize fuel treatments with an integrated forest management plan for the protection of residential areas in Greece. The significant point of our work is the use of a high-resolution ground dataset created for fire risk reduction efforts. The Kassandra Dataset was produced through a labor-intensive and costly process with detailed land cover information that cannot be easily acquired or produced for large landscapes [37]. Similar surveys need to be conducted even for smaller areas in many different WUI zones of Greece because they act as case studies for the informed spatial optimization of fuel treatment management plans. To the best of our knowledge, this work is the first application that considers concurrently so many selection criteria relevant to environmental, socioeconomic, and planning parameters, and one that developed a realistic weight assigning process as it was generated from experts’ opinions and persons who are actively involved in the fuel planning management of the study area. Moreover, a second scenario with equal weighting factors was applied to represent the control scenario.

Key findings of this study include: (a) 19.06% of the Kassandra peninsula can receive treatments for both scenarios, with 13.72% in high-priority areas for scenario 1 and 17.15% for scenario 2; (b) fuel treatments should focus on central and eastern parts of Kassandra where the main WUI areas are located; (c) fuel treatments inside coniferous forest areas have the best prospects to mitigate fire behavior; (d) wildfires initiated from dense Mediterranean shrublands could expose more structures compared to the other land cover types; (e) the fire size does not necessarily relate to structure exposure.

According to our results, in Kassandra’s Forest Service managed lands, with a 10-year fuel treatment plan, an area of 500 ha in high-priority-class pixels should be treated annually. This illustrates that a dynamic workforce and a high budget are required for the scenario’s implementation. However, there is always the possibility of a 10-year period fuel treatment plan to be potentially reviewed or cancelled by future wildfires that may burn areas prior to their treatment [28,61]. With the current fuel treatment rate and their associated costs, as calculated by the actual implementation of the “AntiNero” program in Kassandra, 50 years and a budget of 35 million euros (almost 6000 euros per ha) are required to treat only the high-priority areas identified in Scenario 2. This fact alone highlights how unfeasible, excessive, overpriced, and bureaucratic fuel management is in Greece, especially when compared with published figures of the relevant costs in various countries across the globe [62,63]

Results from this spatial analysis framework could help to provide effective protection of the residential areas by treating fewer hectares, focusing on areas where the results of the simulations show high fire exposure to structures. The main output, that used the fuel treatment grid as a basis for fuel treatment units’ selection, is the Optimized Fuel Treatment Allocation map (Figure 6) that was produced after applying the GIS-MCDA method. This dataset was appropriately used to operationally inform the decisions to allocate fuel treatments in a relatively large area, i.e., the Kassandra peninsula.

As, in the past years, a large-scale wildfire burned the study area resulting in the destruction of residential structures that, in turn, endangered human lives, it is important to identify and monitor high-wildfire-risk areas to pre-organize evacuation plans, to avoid confusion and prevent injuries. The results of this methodological framework, and especially the simulation output of the FTH factor that specifies where fires that affect developed areas are initiated, has already been used by the Greek Fire Service branch of Kassandra by adjusting their annual vehicle allocation for early fire detection and pre-suppression planning. In an alternative approach, more or different criteria could be involved in this study even if it focuses on the same management priority, i.e., population safety in the region. Thus, instead of the FTG simulations output, which is an exclusion criterion that rejects or selects areas for management, other fire indicators such as burn probability or conditional flame length could be used to account for fire intensity or chances that a pixel will eventually be burned by a future fire. Moreover, additional placement criteria may be explored aiming at the reduction in the implementation cost or to avert potential new ignitions. For example, the locations of the national electricity network or the installed solar photovoltaic panels (PV) networks could be involved as points of high ignition risk. Statistically, about 30% of fire ignitions in Kassandra come from short circuits of the electricity supply pole cables [64], while, during the last decade, the wide installation of PV is responsible for fire events that were caused by the dust density around the PV array, the ambient temperature, or the material structure of the PV [65].

The final overall priority maps were verified with more than 70 field surveys by the land managers of the local Forest Service, to confirm or reject the areas that were selected. Thus, our results informed the official final decisions of where fuel treatment should be allocated in the previously mentioned long-term fire protection management masterplan of “AntiNERO”. The process of field verification revealed several issues for consideration. In many cases, the fragmented landscape of Kassandra has forested areas that intermingle with fenced residential areas. This intermix restricts access to forested lands to receive fuel management, leading to the rejection of several initially selected projects. Moreover, the field verification revealed a significant deviation between actual slope measurements and the raster layer of the slope, although it was created by high-resolution large-scale orthophotography’s DEM (2 m × 2 m). Thus, for areas that were rejected in the office due to their steep inclination, the field survey revealed that they could be managed and vice versa.

In the present study, a high-resolution and ground mapping of the land use and vegetation cover was used. However, some limitations related to the fire modeling technique utilized must be considered, e.g., the use of fuel models that were not created for the study area could potentially affect the accuracy of the results, although we ensured that the most appropriate fuel models that can capture the actual fire behavior were selected and used. The use of a meteorological scenario is another limitation of the study. However, we were mainly interested in accounting for the most severe weather conditions that occurred during a wildfire event. At the end, the study also includes the hierarchical classification of the managed areas, considering economical and temporal restrictions.

The fuel treatment plan implementation in the field also highlighted the complexity of forest ownerships in Kassandra, a fact that prevented its uninterrupted implementation in the landscape. Where ownerships are smaller or are with multiple owners or with competing interests, as in many local cases, all these impede fuel management. An effective treatment plan’s implementation presupposes that forest managers can bring together most of these actors and convince them, as they cannot legally force it, to activate a cooperation of all the involved stakeholders toward the common goal of protecting all lands from wildfires. Finally, from a different point of view, intense removal of biomass could subsequently induce depletion of nutrients from the forest [66,67], considered as the fuel treatments’ greatest disadvantage. Therefore, there is a great need for future research to examine fuel treatment plans according to their costs and their ecological impacts to achieve optimal balances and the best results for the Greek forested ecosystems. Nevertheless, in some places throughout the country, it was observed that prejudice and negative reactions to the implementation of preventive treatments by some people or groups occurred, although all regulations and State-imposed instructions to implement the officially approved projects were met by the responsible agencies. This fact may be attributed to the lack of systematic education of the Greek society that could help them to obtain a fundamental perception of the science around wildfires and the proper means of their prevention and values-at-risk or WUI protection [68].

This study addresses the need for a planning protocol or framework that can inform the optimal fuel treatment allocation in a typical Mediterranean area with a dense wildland–urban interface. Our main goal was our implemented and tested framework to become a useful tool for local forest management planning that can be applied to prioritize fuel treatments in target locations (e.g., the WUI or protected areas). Different forest management scenarios and trade-offs can be considered within this framework, such as fire resilience, ecological restoration, or facilitating fire suppression, and even identifying areas that can maximize the levels of achievement of multiple goals [2,28]. In other words, achieving one management goal, e.g., changing the stand structure for reducing wildfire hazard, could, at the same time, achieve ecological targets, timber retrieval, or the aesthetic quality upgrade of the managed forest [56]. Management priorities should be defined considering each region’s special characteristics, i.e., the spatial configuration of fuels, ecological, legal, and social constrains, the available budget, and the location of communities and values-at-risk [12,13].

Additional studies need to examine the performance of appropriate spatial fuel composition and density, and the effect of different governance policies in reducing fire potential. More specifically, emphasis should be placed on a combination of thinning, mechanical treatments for managing or arranging the fuel, and understory burning to reduce surface and canopy fuels to control fire sizes and intensities [69]. The combination of mechanical treatments (mostly thinning, pile burning, and grazing) is a promising choice. Prescribed fire can further reduce implementation costs, to avoid hazardous fuel accumulations. Logging residuals should be either lopped and scattered or piled and burned later when conditions allow it. In addition, pruning of trees that are left after thinning is suggested to break the continuity of the ladder fuels [70].

As fuel treatments in earlier years were limited or had relatively few treatable areas, the forests fuel loads and continuity have extremely increased, and as the prescribed burning is still not legal in Greece, the implementation’s cost for a fuel treatment plan is significantly higher compared to a similar application in other countries where prescribed fire is used instead of the costly mechanical treatments [62,63]. To face the ever-growing wildfire problem in Greece, we proposed to implement landscape-scale fuel management based on defined goals and priorities, as described in this work. Without proper fuel management to mitigate the spread and intensity of large wildfires before they occur, the effectiveness of firefighting is also dramatically reduced [71]. Proactive management of forest fuels has proven to be an effective method for stopping fuel continuity and reducing fuel quantity, lowering the rate of spread and the intensity of fires, and increasing the resilience of forest ecosystems [24].

It is well recognized that the fire-prone Mediterranean forested ecosystems can attain a significant decrease in basic fire behavior attributes when they are properly managed. However, applying optimizations to a fuel treatment planning is an extremely complex problem. This research provides useful guidelines to forest managers interested in applying a similar decision-making process, enabling them to formulate and change their preferences or decide on their final management goals. The applied method offers the most effective utilization for fire management investments aiming at mitigating fire behavior that, in turn, can positively affect wildfire suppression around or inside residential areas. This study is not a hypothetical case study, but it has been tested and applied to guide and resolve real fuel management questions and issues in Kassandra, Greece. The results revealed profound differences on where fuel treatment projects should be applied between the different forest types and risk areas. We intend to monitor and assess the effectiveness of this fuel management plan’s spatial allocation and performance a few years after its application.