Integrating Remote Sensing Methods and Fire Simulation Models to Estimate Fire Hazard in a South-East Mediterranean Protected Area

: Unlike low intensity ﬁre which promotes landscape heterogeneity and important ecosystem services, large high-intensity wildﬁres constitute a signiﬁcant destructive factor despite the increased amount of resources allocated to ﬁre suppression and the improvement of ﬁreﬁghting tactics and levels of organization. Wildﬁres also a ﬀ ect properties, while an increasing number of fatalities are also associated with wildﬁres. It is now widely accepted that an e ﬀ ective wildﬁre management strategy can no longer rely on ﬁre suppression alone. Scientiﬁc advances on ﬁre behavior simulation and the increasing availability of remote sensing data, along with advanced systems of ﬁre detection can signiﬁcantly reduce ﬁre hazards. In the current study remote sensing data and methods, and ﬁre behavior simulation models are integrated to assess the ﬁre hazard in a protected area of the southeast Mediterranean region and its surroundings. A spatially explicit ﬁre hazard index was generated by combining ﬁre intensity estimations and proxies of ﬁre ignition probability. The results suggest that more than 50% of the study area, and the great majority of the protected area, is facing an extremely high hazard for a high-intensity ﬁre. Pine forest formations, characterized by high ﬂammability, low canopy base height and a dense shrub understory are facing the most critical hazard. The results are discussed in relation to the need for adopting an alternative wildﬁre management strategy.


Introduction
Fire is an environmental factor with a long history on earth, dating back to the first appearance of terrestrial vegetation, and as a result it has played a vital role in the ecology and evolution of species, ecosystems and humans [1][2][3][4]. While low intensity fires, resembling prehuman fire regimes, are important drivers for maintaining landscape and habitat heterogeneity and promote a number of important ecosystem services to humans [5], large high-intensity wildfires constitute a significant destructive factor for many biomes and ecosystems across the world. Man has dramatically changed fire behavior and regime, converting it into a primarily anthropogenic factor [6], with more than 464 million ha affected annually by wildfires [7]. In the Mediterranean Europe, in particular, fire affects an average area of 350.000 ha annually, with an interannual variation ranging from a little more than 100 thousand ha to almost one million ha, following fluctuations in weather patterns [8]. Wildfires do not only threaten the ecosystem's ecological integrity, but in recent years they have been responsible Fire 2020, 3, 31 3 of 23 data at high spatial, spectral and temporal resolution offer a great tool to, thematically and spatially, accurately map the complex vegetation mosaics observed in the Mediterranean landscape [35][36][37]. While moderate resolution Landsat data have been widely used for vegetation and fuel mapping, as have various very high resolution data, the suitability of Sentinel-2 data for fuel mapping in the complex landscape of eastern Mediterranean region remains largely unexplored. A recent study by Stefanidou et al. [38], where various landcover types have been identified using Sentinel-2 data in Greece, achieved promising results and suggested that these data could significantly improve the accuracy of fuel mapping at a regional and national scale. Remote sensing data and methods can serve other aspects of wildfire management also. Sifakis et al. [39] presented an automatic fire detection algorithm based on data obtained by the sensor Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) geostationary satellite (henceforth: MSG-SEVIRI), with a spatial resolution of approximately 4 km and a temporal resolution of 15 min. The algorithm was tested during the megafires of 2007 in Greece. The results obtained show an 82% efficiency in detecting fire spots using the Short-Wave Infrared (SWIR) channel 4 and the Thermal Infrared (TIR) channel 9. A similar work has been developed and tested in Australia by Xu and Zhong [40] employing multispectral Himawari-8 imagery. The results reported demonstrate the high accuracy and time efficiency of the method in detecting hot spots. Such tools are of great value for the monitoring of fire evolution, spread direction and dynamics, especially if one takes into account the very high temporal resolution of the data. The only drawbacks of the presented methods lie on the coarse spatial resolutions of the imagery, which exceeds 2 km and restricts the applicability of the method in directing the suppression forces and support real-time firefighting tactics. The integration of such methods with higher resolution data could be a step forward in utilizing remote sensing technologies for time efficient and spatially accurate fire detection [39,40]. Remote sensing data and methods are also widely employed for post-fire mapping of the burned area, as well as for the post-fire characterization of wildfires in terms of fire severity [41][42][43]. Furthermore, the integration of data acquired by active sensors allows for the mapping of the vertical characteristics of the complex forest structure and can assist the process of characterizing and mapping forest fuels and post-fire damage assessment [44].
In the current study we employ remote sensing data, fire behavior simulation modeling and other fire hazard related components to generate a spatially explicit estimation of fire hazard in a protected area in Cyprus using a Fire Danger Index (FDI) estimation formula. The resulting product is expected to assist in the identification of vulnerable areas and promote the adoption of appropriate management measures for hazard reduction. Furthermore, it will be integrated in an automatic fire detection system to assist the process of decision making and mobilization of resources, once a potential ignition is detected.

Study Area
The study site covers the Western part of the Troodos mountains which constitutes the largest mountain range in Cyprus with geographic coordinates 34 • 53 29.06 N and 32 • 51 54.24 E. Mount Troodos is characterized by its high biodiversity and diversity of endemic plants which resulted in its inclusion in the list of 13 "plant diversity hotspots" in the Mediterranean region. The central part of Mount Troodos has been a designated National Forest Park (NFP) since 1992 and it is included in the NATURA 2000 Network of protected areas (CY5000004). The altitudinal range is from 455 to 1963 m a.s.l. Administratively, it falls within the Prefecture of Limasol and the management of the forested part is conducted by the Forestry Directorate of Cyprus.
Precipitation at the NFP of Troodos ranges from 660 mm at the lowest parts and exceeds 1100 mm at its peak altitude. The first snowfall usually occurs at the beginning of December and the last at the middle of April. According to the weather data of the nearby weather station of Saitas, at an altitude of 640 m, the hottest month is July, with an average daily temperature of 26.5 • C and a monthly average maximum temperature of 38.6 • C. The coldest month is February with an average daily temperature of 8.4 • C and an average monthly minimum temperature of −1.5 • C. The wettest month is December with an average precipitation of 130.1 mm, and the driest is July with 6.6 mm precipitation [45].
The area of the NFP of Troodos is dominated by forest ecosystems of Pinus brutia, P. nigra, Juniperus foetidissima and Quercus alnifolia, in pure or mixed formations [46]. The P. brutia forest starts from the lowest altitudinal zones and extends locally up to c. 1.400 m. The shrub understory consists of Q. alnifolia, Arbutus andrachne, Pistacia terebinthus, Styrax officinalis, Genista fasselata. P. nigra, starts at altitudes of 1200 m in mixed formations with P. brutia and extends up to the top altitude, often in mixture with J. foetidissima. Pure stands with J. foetidissima are mainly restricted to the central part, from about 1.300 m up to the highest peak. On rocky slopes the stands are usually open whereas on better sites they are dense. Riparian broadleaved forests occur from the lowest parts up to an altitude of 1.600 m, and consist of Alnus orientalis, Platanus orientalis, Laurus nobilis and Salix alba, while Nerium oleander is sometimes present. A large part of the study area consists of evergreen broadleaved shrubland formations at varying densities, from very dense and 2-3 m heigh t formations to low density short formations.

Remote Sensing Data and Analysis
For the identification of vegetation formations and the generation of a landcover map of the study area a time series of Sentinel-2 images were employed, and analyzed in an object-oriented image analysis (OBIA) environment, using the software eCognition [47]. The time series compilation of Sentinel-2 images consisted of four images acquired on 31/10/2018, 03/02/2019, 25/03/2019 and 03/07/2019, respectively ( Figure 1). The images were selected to cover an entire annual vegetation cycle because this allows for the capture of the vegetation, both natural and anthropogenic, in the different phenological stages. This is of particular importance for mapping dynamic land cover types, such as low-density shrubs, which have a significant grassland cover, low density pines as well as agricultural land with annual crops. All selected images were cloud free.
Sentinel-2 mission consists of two polar-orbiting satellites launched by European Space Agency [48], carrying the Multispectral Instrument (MSI) sensor and delivering data at spatial resolutions between 10 and 60 m, covering a wide part of the electromagnetic spectrum, from visual to infrared (Table 1), and a temporal resolution of 5 days. The images were processed at Level 2A, so they already provided geometrically and atmospherically corrected at Bottom of Atmosphere (BoA) reflectance. For the image analysis, only the spectral bands at 10 and 20 m spatial resolution were employed. Furthermore, the Aster Global Digital Elevation Model (GDEM) was employed, downloaded by the NASA EarthData database. The sharp relief of the study area results in a zonation in the distribution of vegetation types and the use of a DEM was indispensable to achieve maximum classification accuracy. The landcover mapping was done using the software eCognition v. 9.4, which allows the integration in the same classification project data from different sources and different spatial resolutions, without the need for resampling prior to the classification. However, the spectral bands with a spatial resolution of 20 m were resampled at 10 m, and a band composite of 10 bands for each image date was created at a spatial resolution of 10 m. The resampling does not, of course, improve the spatial resolution of the 20 m resolution bands, but it was done simply for practical purposes and for the better organization of the project. OBIA, as implemented in the software ecognition, has several advantages compared to traditional pixel-based approaches when applied on High Spatial Resolution data. It effectively deals with within-class variability, and since the classification is done not on the pixel unit but on the object unit, various shape, texture and context characteristics can be employed in the classification process [49,50]. Furthermore, OBIA avoids the pixelized (salt and pepper) representation of landcover classes which is often observed in pixel-based approaches, while the final classification product can be integrated directly into a vector-GIS for further analysis [51].    Sentinel-2 mission consists of two polar-orbiting satellites launched by European Space Agency [48], carrying the Multispectral Instrument (MSI) sensor and delivering data at spatial resolutions between 10 and 60 m, covering a wide part of the electromagnetic spectrum, from visual to infrared (Table 1), and a temporal resolution of 5 days. The images were processed at Level 2A, so they already provided geometrically and atmospherically corrected at Bottom of Atmosphere (BoA) reflectance. For the image analysis, only the spectral bands at 10 and 20 m spatial resolution were employed. Furthermore, the Aster Global Digital Elevation Model (GDEM) was employed, downloaded by the NASA EarthData database. The sharp relief of the study area results in a zonation in the distribution of vegetation types and the use of a DEM was indispensable to achieve maximum classification accuracy. The landcover mapping was done using the software eCognition v. 9.4, which allows the integration in the same classification project data from different sources and different spatial resolutions, without the need for resampling prior to the classification. However, the spectral bands with a spatial resolution of 20 m were resampled at 10 m, and a band composite of 10 bands for each image date was created at a spatial resolution of 10 m. The resampling does not, of course, improve the spatial resolution of the 20 m resolution bands, but it was done simply for practical purposes and for the better organization of the project. OBIA, as implemented in the software ecognition, has several advantages compared to traditional pixel-based approaches when applied on High Spatial Resolution data. It effectively deals with within-class variability, and since the classification is done not on the pixel unit but on the object unit, various shape, texture and context characteristics can be employed in the classification process [49,50]. Furthermore, OBIA avoids the pixelized (salt and pepper) representation of landcover classes which is often observed in pixel-based approaches, while the final classification product can be integrated directly into a vector-GIS for further analysis [51].
Apart from the original spectral bands the following spectral indices were calculated and used in the analysis: The classification process was assisted by ground truth data, collected in May 2019, in situ, at 50, 20 × 20 m plots, where the land cover type was recorded together with various vegetation structural characteristics, which were later used to allocate fuel models to the various vegetation types. Training data were also collected by visual inspection of the Very High-Resolution Images available on Google Earth. A scale parameter of 50 was selected, using a trial and error approach and visually inspecting the results, for the creation of homogenous objects based on the spectral characteristics included in the segmentation process. Only the spectral bands of 10 m spatial resolution, of all four images, were used for the segmentation. The first step of the analysis was to identify the areas covered by snow in the two images of February and March 2019, respectively. This was achieved by using the NDSI and an altitudinal threshold of 1200 m. The objects identified as snow were merged together and re-segmented with a scale parameter of 50, using only the spectral bands of 10 m resolution of the July 2019 and October 2018 images, and were classified separately from the rest of the scene using spectral bands and indices from the July 2019 and October 2018 images. The total number of created objects was 23,155 with an average size of 0.86 ha, minimum of 0.01 ha, maximum of 13.85 ha and Standard Deviation of 0.86 ha. These objects represent homogenous areas and constitute the basic unit for classification. Classification Trees (CART), Support Vector Machines (SVM) and Random Forests (RF) classifiers, all embedded in eCognition, were tested for their effectiveness in identifying the selected classes. The evaluation of each classifier's performance was assessed by the efficiency in reproducing the training set. Random Forests [52] classification algorithm was found to be the most accurate and was employed for the classification. The RF classifier is a nonparametric machine learning technique which consists of an ensemble of classification trees. Each tree is built based on a random subset of training data and a random subset of predictor variables. Variables that appear more often in the ensemble are those with the highest predictive power. Similarly, each tree predicts inclusion of a training sample in one class. The class predicted by the final model is the one predicted by the highest number of individual trees. The minimum mapping unit was set at 0.25 ha and various classification refinements were made in order to avoid isolated small polygons across the final classification map. The final classification accuracy was estimated using an error matrix on an independent set of 365 points, located randomly across the study area, where the landcover type was verified on the ground.

Fire Behavior Simulation and Estimation of Fire Hazard
As already mentioned, a survey was conducted on 50 plots across the study area, following a stratified random sampling design. The stratification was done using the Corine Land Cover data, which was the only available land cover dataset prior the commence of the study and the development of the landcover map with the method described above. On each plot the relevant abundance of all tree and shrub species was recorded and used to estimate the dominant and accompanying species. Various structural characteristics were also recorded including breast height diameter (BHD), tree height (used to calculate canopy height), live and dead crown base height, crown diameter and canopy cover. For shrubs, instead of BHD, the diameter at the ground level was recorded. Furthermore, litter depth was measured, and the existence of dead branches and their average diameter was also recorded.
The above measurements were used in allometric equations developed for vegetation types similar to the ones observed in the study area [53][54][55][56][57] to obtain an estimate of the fuel load present on each plot, split into diameter categories as required for the description of fuel models. Based on the fuel estimates and the morphological characteristics of each vegetation type, a fuel model (FM) was assigned to each plot. Two custom FM were built, for the high-and low-density pines, respectively, and five standard FMs, described by Scott and Burgan [58], were employed to capture the variation on the vegetation formations present in the study area. Each landcover type was represented by a single FM and the fuel model raster dataset was created.
Normally, at the final stage of any OBIA classification process, adjacent objects of the same landcover type are merged together before exporting. In this case, this was not done, and the original 23,155 objects were kept separate, since each one of them represented a homogenous area in terms of vegetation compositional and structural characteristics. On each of these objects (polygons) the structural characteristics of Canopy Height (CH), Canopy Cover (CC) and Canopy Base Height (CBH) of the nearest plot of the same landcover type and aspect were assigned, using the point-to-polygon spatial join function in ArcGIS 10.7. This method was preferred to calculating the average values of plots falling on the same landcover type because it allowed for the better representation in the relevant mapping products of the high landscape complexity observed in the study area. As a result, three additional raster datasets were created representing CH, CC and CBH, respectively, for the entire study area, which were necessary to perform fire simulation. Three additional raster datasets were created using the Aster GDEM, a Digital Elevation Model, a Digital Aspect Model (DAM) and a Digital Slope Model (DSM). The latter were resampled at 10 m spatial resolution and spatially matched perfectly to the FM raster dataset, because simulation cannot be executed if the spatial data do not have the exact same number of columns and rows. The resampling of the Aster GDEM derived products does not of course lead to any actual improvement in their spatial resolution, but unfortunately there was not a better DTM available for the study area.
The model Flammap [30] was employed to simulate fire behavior. Flammap is a two-dimensional fire simulation model which estimates fire behavior under uniform weather conditions. It calculates important components of fire behavior including fireline intensity, rate of spread, flame height and crown fire activity. Although flammap does not take into account the ignition points or the variation of weather conditions during a fire event, it is a powerful method typically used to estimate fire risk based on the vegetation and topographic characteristics [59]. In that sense, it is more appropriate for a study like this, compared to other simulators, such as Farsite [60], because it allows for the identification of areas with a high potential of giving a high-intensity fire, once appropriate weather conditions exist.
The aim of the study is to assess the fire hazard in the area, in a spatially explicit manner, and under conditions that favor a high-intensity fire. The largest fires that have recently occurred in the Eastern Mediterranean region (Peloponnese Greece) occurred in the summer of 2007, which was the hottest recorded in the area for almost a century, with three consecutive heatwaves from June to August and wind patterns that favored the spread and intensity of fires [61,62]. These extreme weather conditions are more likely to occur more often in the future than to be an exceptional case of extremely rare occurrence [61,62]. For this reason, the fire behavior simulation was performed under the extreme scenario described by Mitsopoulos et al. [32] in the Eastern Mediterranean. The general orientation of the area is west so a west wind was assumed, which promotes upslope fire, while the fuel moisture parameters were set at 5%, 7%, 9%, 80% and 100% for 1-h, 10-h, 100-h, live-woody and live-herbaceous fuel, respectively. In the absence of in situ measurements of crown bulk density (CBD) it was decided to be kept constant at 0.3 kg/m 3 for all fuels, which is the average value observed in an area of similar orientation, continentality and vegetation composition (forest of profitis Ilias) on the Island of Lesvos by Palaiologou [63].
Fire behavior simulation using flammap generates a number of different parameters of fire behavior including Fireline intensity (often called Byram's intensity), rate of spread, flame length etc. Fireline intensity describes the heat release per meter of the fire front and is calculated by the following equation [64] which is an adjustment of the original equation developed by Byram [65]: where: I = Fireline intensity in Kw/m H = Heat yield in cal/g W = Fuel loading in tonnes/ha R = Rate of spread in m/min Fireline intensity estimations were rescaled to a scale between 0 and 1 to form the Fire Intensity (FI) component of the Fire Danger Index (FDI) formula (6) presented by Xofis et al. [66]. Rate of Spread (ROS) is also an important component of fire behavior, since it determines to a great extent the time required for an ignition to turn into a large size and hard to suppress wildfire. A preliminary analysis in two study areas showed that Fireline intensity and rate of spread are significantly and positively correlated to a degree of around 60%. However, various landcover types, such as grasslands and meadows, have a low amount of released energy, due primarily to the low fuel load, but still have a high ROS [66]. In a complex landscape mosaic, such as the one in the study area, and in the wider Mediterranean region, these landcover types impose a high fire hazard due to the increased possibility of a wildfire to quickly spread into a nearby area of high fuel load and high potential for a high-intensity fire. Given that the purpose of the FDI is not only to identify the most vulnerable to high-intensity wildfire areas, but also to increase the level of organization of suppression forces, the ROS is included as a separate component in the FDI formula, with a relatively low weight, despite the fact that it constitutes a component of the Fireline intensity equation, and in a sense it is already included in the FDI formula. ROS values were also rescaled to a scale between 0 and 1 and formed the ROS component of the FDI formula (6) Fire 2020, 3, 31 9 of 23 The above formula integrates two more important components of fire risk, the Human Index (HI) and the Pyric History Index (PH), which serve as proxies of fire ignition probability. The HI attempts to capture the relevant risk for a fire event associated with anthropogenic activities. As a proxy of anthropogenic activities, the distance to roads was adopted, with values ranging from 0, for areas at a distance of 200 m or more from roads, to 1 for areas at immediate proximity to them. Distance to roads has been reported to be positively associated with ignition frequency [67,68]. The Pyric History does not have a direct relationship with fire hazard, especially if one takes into account that a past wildfire may reduce the fuel load and subsequently the intensity of a fire event. However, it provides an indication of the fire pattern of an area which might indicate trends related to specific land uses and spatial locations. For instance, in the Eastern Mediterranean region, where free range pastoralism is still practiced in the mountainous and semi-mountainous regions, it is very common for deliberate fires to be set by shepherds for the improvement of grazing conditions. This trend results in a clustering of past fire episodes indicating the higher risk for fire in such areas. Although these fires rarely become high-intensity ones, because they primarily occur in autumn and in areas with a low fuel load, the possibility to spread into nearby more flammable ecosystems and vegetation formations always exists. The inclusion of the PH index in the above formula will allow for the identification of areas where a specific spatial pattern of past fires does exist. The PH was calculated using data on ignition events over the last 25 years provided by the local Forest Service. A Kernel Density Estimation Function was applied to convert the point ignition data to a raster file with values ranging from 0, for areas away from a past ignition point, to 1, for areas at immediate proximity to a past ignition. The FDI calculated with the above formula varies between 0 and 1 with higher values indicating a higher fire hazard.

Validation of Fire Hazard Estimation
To validate the results of the estimated fire hazard using the FDI formula, we compared them with the risk estimated using the Conditional Flame Length (CFL) index, which has been widely used as an estimate of fire hazard [28,35,69]. For the CFL estimation a fire simulation under the Minimum Travel Time (MTT) algorithm, embedded in Flammap, with 10.000 randomly distributed fire ignitions and a fire duration of 8 hours was performed. The simulation was done under the exact same burning conditions as the first simulation, performed for the estimation of FI and ROS. The MTT algorithm [29] employs the Huygen's principle to replicate fire growth, where the growth of the fire edge is a vector or wave front. Since the burning conditions are constant, the MTT algorithm reveals the effect of fuel conditions and topography on fire growth [70]. The CFL was estimated using the ArcFuels 10 software, which was embedded as an extension in ArcGIS 10.7, using as input data the flame-length probabilities (FLP) metrics file generated in Flammap. CFL is an estimate of the mean flame length (FL) of all the ignitions (10.000) that burned a particular pixel.
The calculated FDI values were classified into four classes which were defined as very low hazard (0 < FDI ≤ 0.1), low hazard (0.1 < FDI ≤ 0.35), high hazard (0.35 < FDI ≤ 0.7) and extremely high hazard (FDI>0.7). The lowest class incorporates unburnable land cover classes and areas where some vegetations exist but where it is impossible to sustain a high-intensity fire. The upper threshold of the second class was chosen based on the relevant weight of FI and ROS in the FDI. Both these components account cumulatively for the 70% of the FDI. Since areas that score bellow 0.5 in these components carry a lower risk of a high-intensity fire and the weight of those components in the FDI is 0.7, the resulted threshold was set at 0.35. This class incorporates areas where the intensity and rate of spread of a possible fire is low, so even in the case of a fire event the probability of ending up in a high-intensity fire is low. The third class incorporates all these areas that score highly in the components of FI and ROS. The upper class represents the most vulnerable to wildfire areas, since in addition to the favorable fuel characteristics, these areas are also close to roads and/or past fire spots.
The validation was made on 388 points distributed across the study area, under a stratified random sampling design and with a minimum distance of 300 m between them, using the Sampling Design Tool extension in the ArcGIS 10.7 platform. The request was for 500 points, but only 388 could be located under the restrictions of the sampling design and the minimum distance. On each generated point the FDI class was assigned as well as the CFL value that corresponded to the particular location. One-Way Analysis of Variance (ANOVA) and Post-hoc Tukey HSD (Honestly Significant Difference) tests were employed to test for statistical differences in the CFL among the four classes of fire hazard.

Assesment of between Landcover Types Differences in Estimated Fire Hazard
From the 388 points selected for the validation of the estimated fire hazard, 353 points corresponded to burnable landcover types. On each one of these 353 points the estimated FDI value was assigned together with the landcover type, the vegetation structural characteristics (CC, CBH, CH) and topographic data. One-Way ANOVA and Post-hoc Tukey HSD tests were employed to test for significant differences in the estimated fire hazard between landcover types.
For the points corresponding to high and low-density pines landcover types a second step of exploratory analysis was done using Regression Trees (RT), in order to assess the effect of topography and vegetation structural characteristics on the estimated fire hazard. Regression Trees is a data mining method which handles simultaneously continuous and categorical independent variables and they are very effective in handling non-linear and non-hierarchical relationships. Regression trees repeatedly divide the data into two mutually excluded groups using one independent variable at a time, until homogenous groups, in terms of the dependent variable, are achieved or the data cannot be divided any further [71].

Results
Seven classes were identified in the study area, shown in Figure 2 and Table 2. The overall accuracy of the final product was estimated at 93% and the Kappa statistics was calculated at 0.9, indicating an excellent performance by the classifier (Table A1 in Appendix A).
Low-and high-density pines collectively account for more than 50% of the total area, consisting of pure and mixed stands of P. brutia and P. nigra, often with a dense understory of shrubs. Figure 2 demonstrates the zonation in the distribution of landcover types observed in the study area. The lowest altitudes are occupied by agricultural land, and low and mid altitudes by scarce and low shrublands, which locally turn into dense shrublands, and as the altitude increases the pines dominate. The highest altitudes are dominated by J. foetidissima formations often in mixture with pines. Broadleaved forests tend to occur only along the various streams present in the study area.
Following the process described above, on each polygon generated in OBIA which represents a homogenous, in terms of vegetation structure and composition area, one fuel model was assigned, together with the additional structural characteristics (CH, CC, CBH) of the nearest sample plot of the same land cover class and aspect. Two custom fuel models and five standard fuel models were employed, and the corresponding fuel loads are shown in Table 3. For the custom fuel models the properties of FM type, Surface-Area-to-Volume (SAV) ratio, fuel depth, dead fuel extinction moisture and heat content were adopted by Palaiologou [63].   Low-and high-density pines collectively account for more than 50% of the total area, consisting of pure and mixed stands of P. brutia and P. nigra, often with a dense understory of shrubs. Figure 2 demonstrates the zonation in the distribution of landcover types observed in the study area. The lowest altitudes are occupied by agricultural land, and low and mid altitudes by scarce and low shrublands, which locally turn into dense shrublands, and as the altitude increases the pines dominate. The highest altitudes are dominated by J. foetidissima formations often in mixture with pines. Broadleaved forests tend to occur only along the various streams present in the study area. Flammap simulation and the resulted layers of Fireline intensity and rate of spread revealed the diversity of burning conditions in the area (Figure 3a,b). FI exceeds the threshold of 0.5 in more than 42% of the area, while the ROS exceeds the value of 0.5 in more than 53% of the study area.
The highest values of FI seem to coincide with the distribution of pine forest formations. High ROS are also observed in the distribution zone of pines, but the J. foetidissima formations at the highest altitudes also appear to score high in the ROS index. The other two components of the FDI (Figure 3d,c), which are mainly related to the possibility of a fire event, refer to a much smaller proportion of the study area. The integration of these four components using the FDI formula resulted in the final, spatially explicit, index of fire hazard in the study area (Figure 4).    The value of FDI ranges between 0 in the lowland parts of the study area, occupied primarily by non-burnable classes, to 1 in the uplands, where forest formations exist, often in close proximity to roads or past fire spots. The FDI reveals a complex mosaic of areas with different flammability properties and different degree of fire hazard, reflecting the compositional and structural vegetational diversity, as well as the complex relief of the study area. The classification of the FDI into four classes of fire hazard ( Figure 5) indicates that 36% of the study area faces a high fire hazard and an additional 14% an extremely high hazard. In total, 50% of the study area faces a low or very low fire hazard. The value of FDI ranges between 0 in the lowland parts of the study area, occupied primarily by non-burnable classes, to 1 in the uplands, where forest formations exist, often in close proximity to roads or past fire spots. The FDI reveals a complex mosaic of areas with different flammability properties and different degree of fire hazard, reflecting the compositional and structural vegetational diversity, as well as the complex relief of the study area. The classification of the FDI into four classes of fire hazard ( Figure 5) indicates that 36% of the study area faces a high fire hazard and an additional 14% an extremely high hazard. In total, 50% of the study area faces a low or very low fire hazard.
The ANOVA results for the validation of the FDI (Figure 6) show how the four classes of fire hazard, as estimated using the FDI formula have significantly different means in the estimated CFL, which is adopted in the current study as an independent index of fire hazard. All fire hazard classes differ significantly between them, as tested using the post-hoc Tukey HSD test. CFL increases dramatically between the two classes of low and high fire hazard while the increase between the high and extremely high hazard is smaller but still significantly different. Although CFL is impossible to have captured the increased risk of fire hazard imposed by the anthropogenic activities and past fire occurrence, the observed difference between the two high hazard classes is probably the result of the highest consistency, in terms of estimated fire behavior parameters, on the areas belonging to the extremely high hazard class. The ANOVA results for the validation of the FDI ( Figure 6) show how the four classes of fire hazard, as estimated using the FDI formula have significantly different means in the estimated CFL, which is adopted in the current study as an independent index of fire hazard. All fire hazard classes differ significantly between them, as tested using the post-hoc Tukey HSD test. CFL increases dramatically between the two classes of low and high fire hazard while the increase between the high and extremely high hazard is smaller but still significantly different. Although CFL is impossible to have captured the increased risk of fire hazard imposed by the anthropogenic activities and past fire occurrence, the observed difference between the two high hazard classes is probably the result of the highest consistency, in terms of estimated fire behavior parameters, on the areas belonging to the extremely high hazard class.  The ANOVA followed by a Tukey test revealed, as it was expected, that high density pine formations face the highest fire hazard among the land cover types present in the study area ( Figure  7). High fire hazard is also faced by low density pine and dense shrubs, while surprisingly high fire hazard is also faced by the J. foetidissima formations. Scarce and low shrubs on the other hand, are the least prone to high-intensity fires by land cover type, along with broadleaved forests, where no  The ANOVA followed by a Tukey test revealed, as it was expected, that high density pine formations face the highest fire hazard among the land cover types present in the study area (Figure 7). High fire hazard is also faced by low density pine and dense shrubs, while surprisingly high fire hazard is also faced by the J. foetidissima formations. Scarce and low shrubs on the other hand, are the least prone to high-intensity fires by land cover type, along with broadleaved forests, where no significant difference was found between them, despite the higher FDI value observed in broadleaved forests. Figure 6. Conditional Flame Length (CFL) of the four identified fire hazard classes estimated using the FDI formula. Vertical bars denote ± one standard error.
The ANOVA followed by a Tukey test revealed, as it was expected, that high density pine formations face the highest fire hazard among the land cover types present in the study area ( Figure  7). High fire hazard is also faced by low density pine and dense shrubs, while surprisingly high fire hazard is also faced by the J. foetidissima formations. Scarce and low shrubs on the other hand, are the least prone to high-intensity fires by land cover type, along with broadleaved forests, where no significant difference was found between them, despite the higher FDI value observed in broadleaved forests.  The results of the Regression Trees analysis indicate that both, topographic and vegetation structural characteristics affect the degree of fire hazard (Figure 8). Canopy cover, as expected, determined to a great extent the fire behavior, with pine formation of very low density facing a low fire hazard. Aspect and slope are the two topographic characteristics with the greatest effect on fire hazard. Southeast, south and southwest slopes generally result in lower values compared to the rest of the possible aspects. A positive relationship between slope and fire hazard is also observed under every possible aspect. Canopy base height is a very important vegetation structural characteristic in determining the possibility of a high-intensity fire. Even on relatively gentler slopes, a lower CBH, which does not exceed 3.35 m, can lead to a higher fire hazard in pine forests.
determined to a great extent the fire behavior, with pine formation of very low density facing a low fire hazard. Aspect and slope are the two topographic characteristics with the greatest effect on fire hazard. Southeast, south and southwest slopes generally result in lower values compared to the rest of the possible aspects. A positive relationship between slope and fire hazard is also observed under every possible aspect. Canopy base height is a very important vegetation structural characteristic in determining the possibility of a high-intensity fire. Even on relatively gentler slopes, a lower CBH, which does not exceed 3.35m, can lead to a higher fire hazard in pine forests.

Discussion
In this study we integrated remote sensing data and methods, fire behavior simulation, and proxies of ignition risk to assess the hazard imposed by wildfires in a protected area and its surroundings in the South-East Mediterranean. The results and the approach presented here can provide a significant tool towards a more effective fire management strategy. The spatially explicit fire hazard classification allows its integration into an automatic fire detection system, consisting of land optical cameras and thermal cameras loaded on Unmanned Aerial Vehicles (UAVs) and it is

Discussion
In this study we integrated remote sensing data and methods, fire behavior simulation, and proxies of ignition risk to assess the hazard imposed by wildfires in a protected area and its surroundings in the South-East Mediterranean. The results and the approach presented here can provide a significant tool towards a more effective fire management strategy. The spatially explicit fire hazard classification allows its integration into an automatic fire detection system, consisting of land optical cameras and thermal cameras loaded on Unmanned Aerial Vehicles (UAVs) and it is expected to improve its detection accuracy. It can also be integrated in a decision support system to direct the fire fighting forces in a manner that will reduce the time between ignition and first announcement, increasing the efficiency of the initial attack. Satellite-based automatic fire detection systems have already been developed and tested for their efficiency in fire prone areas [39,40]. These systems, despite the coarse spatial resolution of the imagery (i.e., MSG-Seviri and Himawari-8), are of high operational value in wildfire management [72] since they can detect fires as small as 0.1 ha [73]. The only objection to those systems is the difficulty in the precise positioning of the fire event within the large sized pixel. The integration of the spatially explicit fire hazard index developed here into the raw data processing chain may allow for the identification and allocation of new fire spots with higher spatial accuracy. The identification of the vegetation formation that faces the highest fire hazard is also an important outcome of the study and the adopted approach. Starting from the lowlands to the uplands, the fire hazard in the study area increases. The low altitude areas are dominated by the low hazard formations of scarce and low shrubs, which are often intermixed with agricultural land and other anthropogenic activities, including pastoralism. The latter is often associated with frequent fires of a small area and low heat yield, burnt for the improvement of grazing conditions. The intensive human presence in this area has resulted in the degradation of vegetation which currently consists of phryganic-type formations with a scarce presence of grazing resistant shrubs. The fuel load is kept low, and, as a result, the possibility of a high-intensity fire in these areas is low. In the absence of grazing or any other degradation factors the vegetation turns into dense shrublands, dominated by Q. alnifolia which constitute one of the most hazardous vegetation formations in the study area. Dense shrublands have been reported to be among the most fire prone vegetation types and the main carriers of fires in the Eastern Mediterranean, with a high contribution in the annually burned area [18,19]. Furthermore, Moreira et al. [74] reported a higher probability of an ignition to turn into a big fire when it occurs in shrubland forest formations.
Broadleaved forests face a relatively low hazard by wildfires. In the study area these forests are mainly confined across the streams, from high altitudes to the lowlands. These forests, despite the relatively high dead fuel load, which is generated by the annual drop of leaves and the relatively slow decomposition rate of the Mediterranean region, have humid conditions that locally prevail in these microhabitats, which increase the fuel moisture content and prevent an ignition from turning into a high-intensity fire.
Low-and high-density pine formations are the most vulnerable to fires and are those with the highest possibility to result in a high-intensity fire, and subsequently face the highest fire hazard. The largest and the most catastrophic wildfires that have occurred in the Eastern Mediterranean region over the last 15 years have primarily burned pine forests, dominated by P. halepensis and P. brutia [18,19]. All conifers and especially the Mediterranean pines P. halepensis and P. brutia have a high resin content, both in the woody parts and in the needles, which result in high flammability compared to other vegetation types. At the same time, the shrub understory, rich in flammable essential oils, promotes flammability and increases the fire hazard [75]. Furthermore, the canopy height of pine forests, which often exceeds 20 m, allows fire brands to disperse in longer distances promoting the advancement of fire front by spotting.
Another important observation in the current study is the significant role of vegetation structural characteristics in determining the fire hazard in the vegetation formations of pine forests most prone to high-intensity fires. CC is a major determinant of fire hazard with low tree densities being associated with low fire hazard. CBH is also important in determining the fire hazard in pine forests, with a CBH that exceeds 3.35 m resulting in a lower hazard on relatively gentler slopes. The trend of Mediterranean pines to retain dead branches in the crown and along the trunk has been thought of and reported to be a fire adaptive trait which promotes fire as a mechanism to create a favorable regeneration environment (niche construction) for the seeds retained in serotinous cones in the crown [76,77]. However, it is more likely that this trend is related to the slow decomposition rate of the Mediterranean pine forests rather than a niche construction mechanism. Nevertheless, the negative correlation between CBH and fire hazard provides an indication on the management measures that could prevent large high-intensity fires. Sheltered fuelbreaks on both sides of roads and, in some cases, in the cores of forest patches, where the CC would be kept low and the CBH high, have been proposed by Xanthopoulos et al. [78] as a method for the horizontal disruption of fuel continuity in fire-prone ecosystems. Such a practice would prevent the evolvement of an ignition into a high-intensity crown fire without exposing the soil to conditions of high soil erosion and degradation, often associated with clear cut firebreaks. This aspect was also found to have a significant role in determining the fire hazard with southern slopes scoring lower values. This was a little surprising since southern slopes are generally drier which should favor high-intensity fires. However, the dry conditions of southern slopes most likely prevent the accumulation of high fuel loads and, as a result, reduces the risk of a high-intensity fire. Estimating fire hazard and identifying the most vulnerable areas to high-intensity fires is probably the necessary first step towards developing a more effective fire management strategy. Traditional approaches of wildfire management which rely almost exclusively on fire suppression have reached their limits of success [79], despite the improvement in fire-fighting tactics and methods and the increased amount of resources allocated to it [80]. An alternative strategy is suggested by many experts, which should be based on a better balance between fire suppression, an increased level of preparedness involving fuel management and affordable cost in the burned area [81][82][83]. Fire-smart forest management has been proposed as a more effective approach of fire management, which employs forest management practices, including fuel reduction, conversion and isolation to decrease the probability of an ignition to turn into a high-intensity fire and increase the resistance of forests stands to fires [79,84]. Despite the uncertainties behind the effect of projected climate change on fire regimes, particularly with regards to the changes in the fire landscape ("firescape"), it seems that wildfires will become wilder in the future. Prescribed burning, tested primarily under fire simulation scenarios, seems to be reducing the intensity of future wildfires and, with the exception of the WUI, it deserves higher attention as a fuel reduction measure [84,85]. The same purpose of reducing fuel load could also be served by biomass extraction for bioenergy, especially when the intensity of extraction is high, and it is coupled with an intensive suppression effort [86].

Conclusions
The current study reveals the critical fire hazard faced by the most valuable protected area of Cyprus and one of the most valuable areas in the Mediterranean region. In total, 50% of the entire study area and the great majority of the designated FNP are in the high and extremely high class of fire hazard. Its protection relies on preventing megafires that could result in a significant loss of habitats, possible properties and infrastructure. The FDI employed in the study for the estimation of fire hazard is composed by components that can be easily interpreted, in relation to fire behavior and fire risk. It also integrates components related to local peculiarities in terms of ignition risk, so it also integrates local knowledge. The data required (satellite images and topographic data) are freely available while the software used are well documented and widely used. As a result, the proposed approach can easily be transferred to other fire-prone regions of Europe and elsewhere to support an effective wildfire management strategy. Increasing the level of preparedness through the appropriate allocation of resources and the adoption of management measures which could reduce the fire intensity and rate of spread, are steps in the right direction. When fuel management measures are adopted it is essential to ensure the ecological integrity of the ecosystems and avoid measures that could lead to significant degradation. Furthermore, the integration of the resulted mapping products in automatic fire detection systems could decrease significantly the time required for the initial attack. When wildfire management strategy is developed, one should always bear in mind that the main objective should not be the complete elimination of fires, which is both impossible and ecologically dubious, but rather the establishment of a fire regime with low intensity fires and, most importantly, a low socioeconomic and ecological cost [87].