Mapping forest fuels through vegetation phenology: the role of coarse-resolution satellite time-series.

Traditionally fuel maps are built in terms of 'fuel types', thus considering the structural characteristics of vegetation only. The aim of this work is to derive a phenological fuel map based on the functional attributes of coarse-scale vegetation phenology, such as seasonality and productivity. MODIS NDVI 250 m images of Sardinia (Italy), a large Mediterranean island with high frequency of fire incidence, were acquired for the period 2000-2012 to construct a mean annual NDVI profile of the vegetation at the pixel-level. Next, the following procedure was used to develop the phenological fuel map: (i) image segmentation on the Fourier components of the NDVI profiles to identify phenologically homogeneous landscape units, (ii) cluster analysis of the phenological units and post-hoc analysis of the fire-proneness of the phenological fuel classes (PFCs) obtained, (iii) environmental characterization (in terms of land cover and climate) of the PFCs. Our results showed the ability of coarse-resolution satellite time-series to characterize the fire-proneness of Sardinia with an adequate level of accuracy. The remotely sensed phenological framework presented may represent a suitable basis for the development of fire distribution prediction models, coarse-scale fuel maps and for various biogeographic studies.

Introduction

Phenological studies of vegetation can be carried out at the species level (bud break, flowering, leaf flush, etc.) using in-situ field techniques [1]. However, only remote sensing can offer information on ecosystem phenology and productivity over several temporal scales and continuous spatial coverage. Space-borne optical sensors such as Moderate Resolution Imaging Spectroradiometer (MODIS) provide daily measurements of a variety of biophysical parameters of the land surface [2]. MODIS is a key instrument aboard the Terra and Aqua NASA satellites viewing the entire Earth's surface every 1 to 2 days, acquiring data in 36 spectral bands (http://modis.gsfc.nasa.gov/about/). These earth observation systems examine coarse-scale phenomena that allow retrievals of whole-system phenological metrics, such as the timing and magnitudes of greening, peak activity, and drying phases of the growing season [1]. At this landscape scale, the spatio-temporal dynamics of the vegetated land surface is often referred to as ‘land surface phenology’, as the remotely sensed phenology deals with mixtures of land cover and vegetation [3]. Land surface phenology thus provides information on key aspects of vegetation functionality, such as seasonality, productivity and inter-annual variability. Also, as vegetation functioning responds faster to environmental change and variability than vegetation structure and composition, there is a high potential for characterizing, classifying and mapping vegetation based on remotely sensed phenology [4-6].

Although, a number of studies have used remotely-sensed indices, such as the Normalized Difference Vegetation Index (NDVI), for the monitoring of vegetation dynamics from regional to global scales [7-9], the capability of remote sensing indicators to capture ecosystem functional dynamics has been only partially faced so far. The NDVI has been usually considered a reliable indicator of plant biomass and vegetation primary productivity [10-12]. Therefore, NDVI time series have been routinely used to measure vegetation dynamics and ecosystem phenology over large geographic areas [13-15]. Remote sensing derived vegetation phenology may provide an integrated measure of ecosystem responses to climatic factors, such as temperature and rainfall, as well as a synthetic quantification of the determinant conditions for fires and other disturbances [6].

Phenological classification characterizes and stratifies the land surface based on similar phenological patterns. White et al. [16] developed a global pheno-region database as a geographic framework for studying global climate change. They used Advanced Very High Resolution Radiometer (AVHRR) NDVI time series (8 km of pixel resolution) from 1982 to 1999 in conjunction with an eight-element monthly global climatology to generate global pheno-regions with a minimized probability of non-climatic forcing. Hargrove et al. [17] derived 15 phenological ecoregions based on clustering the similarities in five years (2002–2006) of cumulative MODIS NDVI data (250 m of pixel resolution). Gu et al. [18] performed a phenological classification of the conterminous United States based on MODIS NDVI time-series and a digital elevation model. The resulting pheno-class map is composed of 40 pheno-classes, each having unique phenological and topographic characteristics. Clerici et al. [19] explored the use of MODIS NDVI-derived phenology metrics for the identification and classification of Forest General Habitat Categories in Europe. Finally, De Angelis et al. [20] also used MODIS NDVI-derived phenology metrics to classify the vegetation of Sardinia (Italy) into clusters showing uniform fire behavior.

In the framework of fire risk assessment, variations in satellite-based NDVI values proved to be indicative of variations in water and nutrient availability, plant disease, and other stress factors, which are in turn indicators of a marked vulnerability of the vegetation to fire [21]. For this reason, over the years, many authors have proposed the use of multitemporal NDVI profiles to assess the proneness of vegetation to fire (for a detailed review, see [22]). Accordingly, as vegetation phenological status represents the main driver affecting fuel availability and moisture content, any investigation on fire behavior over large areas requires the ability of capturing spatio-temporal changes in vegetation phenology [20, 23].

Fire managers have tried to summarize the physical parameters and spatial distribution of fuel in different fuel types or fuel models [24-25]. More specifically, a fuel model has been defined as “an identifiable association of fuel elements of distinctive species, form, size, arrangement, and continuity that will exhibit characteristic fire behavior under defined burning conditions” [26]. Two well-known fire behavior fuel type systems are the Northern Forest Fire Laboratory (NFFL) system and the Canadian Forest Fire Behavior Prediction (FBP) system [27]. Also european researchers developed a new system, in the framework of the Prometheus project, which is better adapted to fuels found in Mediterranean ecosystems. These fuel types were defined for surface fire modelling, taking into account fuel height and density. However, traditional fuel maps, being built in terms of fuel types, represent the static, structural properties of fuel [28], not involving its dynamic, functional characteristics, such as moisture content, or vegetation health status.

In the wake of the work of De Angelis et al. [20], the aim of this paper is to use MODIS NDVI time-series (2000–2012) to develop a phenological fuel map of Sardinia (Italy) that takes into account the functional (rather than the structural) aspects of fuels, and to demonstrate its potential in characterizing fire risk (i.e. the chance of a fire to start).

Study Area

Sardinia is an island located in the Mediterranean Sea between 38° 51' N and 41° 15' N latitude and 8° 8' E and 9° 50' E longitude, covering roughly 24000 km2. The island is characterized by a complex physical geography, with a prevalently hilly topography. The highest elevation is 1834 m, average elevation is 338 m. The island has a Mediterranean climate along the coast and a more continental climate on the interior mountain ranges. The average temperature is between 11 and 17° C, rainfalls have generally a Mediterranean distribution, with dry summers and two peaks of rainfall in spring and autumn. Mean annual precipitation is about 800 mm.

Land use along the coast and the main river valleys is dominated by agriculture that covers about 45% of the study area. Most urban areas of Sardinia are located in the coastal zone. In the interior areas, forest stands combined with pastures and shrublands prevail. The principal formations include Quercus ilex and Quercus suber forests. At higher elevations the sclerophyllous oak forests merge with broadleaved forests of Quercus pubescens and Castanea sativa.

The island is considered a wildfire hotspot, with a mean of more than 2500 events per year [29-30]. Fires are mainly of limited size, more than 67% of the events are smaller than 1 ha. Although most of the recorded fires are human-caused rather than ‘natural’ in origin, fire is a process largely correlated with climate. Like in most Mediterranean regions, in Sardinia fire is strongly seasonal, with the majority of fires occurring during the hottest and driest portion of the year. More than 98% of fires occur between May and October with a peak of fires (54%) in July and August.

Data

In this study, information on the remotely sensed phenology of the vegetation of Sardinia was extracted from the 16-day MODIS 250 m NDVI maximum value composite product (MOD13Q1; https://lpdaac.usgs.gov/products/modis_products_table/mod13q1). Twenty-three NDVI images per year from 2000 to 2012 were acquired, for a total of 299 MODIS images. Late fall and winter scenes revealed several low quality and noisy pixels regions due to bad weather events. Therefore, we focused only on late spring-early autumn scenes, i.e. eleven yearly MODIS NDVI images (143 images in total) from Julian day 113 (April 23rd) to Julian day 273 (October 1st), which correspond to the months of highest fire occurrence. By averaging the NDVI values of each 16-day composite over the period 2000–2012, for each image pixel we derived a mean NDVI temporal profile showing mean seasonal variations in the NDVI values Wildfire data were provided by the Forest Service of Sardinia and contain all fires recorded from 2000 to 2010, together with the geographic (UTM) coordinates of their ignition points, and is assumed to be complete and reliable down to the smallest fires. The database used in this study includes all the fires recorded from April 23rd to October 1st during the period 2000–2010, and contains 28744 events.

Methods

The following procedure was used to develop the phenological fuel map: (i) image segmentation on the Fourier components of the NDVI profiles to identify phenologically homogeneous landscape units, (ii) cluster analysis of the phenological units and post-hoc analysis of the fire-proneness of the phenological fuel classes (PFCs) obtained, (iii) environmental characterization (in terms of land cover and climate) of the PFCs.

Fourier analysis and image segmentation

The vegetation (i.e. fuel) seasonality can only be captured by repeated measurements over time. Serial correlation among successive observations taken over a period of time reduces the statistical utility of multitemporal imagery [31]. Therefore, data reduction (usually ordination methods, like Principal Component Analysis) is typically applied to remove image correlation. However, in the ordination process, explicit reference to vegetation seasonality is lost, such that data reduction is achieved at the expense of ecological relevance. Alternative methods for retaining information on vegetation seasonality include polynomial functions [32] and Temporal Fourier Analysis (TFA) [23, 33–34]. Therefore, to analyze remotely-sensed fuel seasonality, we applied Temporal Fourier Analysis to the average NDVI profiles of each pixel.

TFA allows a complex signal to be expressed as the sum of a series of sine and cosine waves (harmonics) and an additive term [33]. The constant additive term corresponds to the arithmetic mean of the time series, while the harmonics correspond to the amplitudes into which the temporal signal is decomposed. TFA transforms a series of N observations taken at regular time intervals into a spectrum of (uncorrelated) harmonics at frequencies of 1, 1/2, 1/3, …, 2/N times the duration over which the N observations were made [31]. The obtained spectrum is a true partition of the variance of the series among frequencies, therefore spectral analysis, like TFA, represents a type of variance decomposition analysis [35]. In TFA there are as many harmonic variables as the data points. However, in practice only a few harmonics usually contribute substantially to the overall variance, and only these need to be calculated for capturing the main features of the vegetation seasonality [23, 31].

In this paper, we used late spring-early autumn NDVI profiles as our basic unit for TFA, the amplitude at one cycle period A1 (i.e. the first Fourier component) thus corresponds to the maximum variability of the NDVI values over the whole fire season, while the second Fourier component A2 summarizes the NDVI variability within the half-period of the NDVI profile. Accordingly, significant amplitude values associated to the first and second periodic terms indicate major variability of the NDVI profiles in response to the seasonal patterns of bioclimatic conditions. Finally, the additive term A0 is related to the summation of the NDVI values over the whole temporal profile [36] and can therefore be considered as a good indicator of the coarse-scale net primary productivity of vegetation on a seasonal basis. For additional details on Fourier analysis, see e.g. [37].

We used the TFA module of IDRISI Selva (Clark Labs, Worchester, MA) to extract the constant additive term (A0) and the amplitude values for the first and second Fourier harmonics (A1 and A2) of the seasonal NDVI profiles of each pixel, with 11 as number of time series images and 16 days as time lag between images. The constant and the first two amplitude output images of the TFA were then used for segmenting the territory of Sardinia into phenologically homogeneous units in the multispectral Fourier feature space.

Image segmentation was conducted using the Fractal Net Evolution Approach [38] of the Definiens eCognition software (Definiens Imaging, Germany). In this bottom-up segmentation technique, individual pixels are perceived as the initial regions, which are sequentially merged pairwise into larger ones, with the intent of minimizing the heterogeneity of the resulting objects. The sequence of the merging objects, as well as the size and shape of the resulting objects, are empirically determined by the analyst [39]. A crucial unitless threshold, named scale parameter, is determined to specify the maximum allowed increase in the heterogeneity after a pairwise merge of the objects [40]. For the segmentation of the three TFA components two criteria were used: color (i.e. the value of each Fourier component for each pixel) and shape. The selection of the best value for each segmentation criterion was based on trial-and-error and on the visual interpretation of the results. Scale parameter is considered as the most effective parameter that is used to control the average image object size [41-42]. It depends on spatial resolution of the corresponding image and features of the study area. The higher value is selected for scale parameter, the larger the object is obtained. On these basis, the following values were selected for all three TFA components: scale parameter = 500, shape criterion = 0.2, and compactness = 0.8. The segmentation produced 60 phenological units (PUs) of at least 25 pixels in size, homogeneous in terms of mean seasonal NDVI timing (Fig. 1).

Fig 1

Spatial distribution of the 60 phenological units (PUs) and of the four phenological fuel classes (PFCs).

To quantify the fire-proneness of each PUs, the coordinates of the ignition points of all fires that occurred in the period 2000–2010 were overlaid on the PUs and the number of fires within each PU was computed. Next, we calculated a summary statistic of fire risk as the proportion of fires in each PU divided by the relative area of that PU in the study site [20]. This ratio, named ‘selection ratio’ (σ) [43, 29], takes values in the range [0, ∞]. If all PUs were equally fire prone, fire ignitions would occur randomly across the landscape, such that the number of ignitions within a given PU would be nearly proportional to the relative area of that PU and σ ≈ 1. Values larger than 1 denote PUs where ignitions occur more frequently than expected by chance alone, values lower than 1 denote PUs where ignition is less frequent than expected [29].

Cluster analysis

In order to identify groups of phenological units with similar seasonality, we performed a hierarchical cluster analysis on the mean values of the Fourier components of each PU. Clustering is an agglomerative classification method for determining the intrinsic grouping in a set of data and identifying subpopulations of related objects with similar behavior with respect to some user-defined ecological property [44]. The cluster analysis was performed within the statistical software package Past 3.04 (free downloadable from http://folk.uio.no/ohammer/past/) using the Euclidean distance as similarity measure and the UPGMA method (average linkage between groups) as linkage criterion (see [35]). This analysis produced a dendrogram in which three distinctive phenological fuel classes (PFCs) homogeneous in terms of phenological patterns were detected: PFC1, PFC2 and PFC3. Based on our territorial knowledge, we divided the largest cluster into two subgroups: PFC2a (mainly distributed in the coastal zones) and PFC2b (mainly located in the inner parts).

Next, to evaluate the pyrological uniqueness of the four PFCs, we performed a permutational analysis of variance on the selection ratios of all PUs. Like for traditional univariate ANOVA, we calculated an F-value by dividing the between-group variability of the selection ratios by the corresponding within-group variability. A p‐value is then computed by permutation of the PUs among the different phenological classes (9999 randomizations, one-tailed test, see [45]).

Finally, in order to test the fire-proneness of the four PFCs, we first calculated the selection ratio of each class as the ratio between the relative proportion of fires in each class and the relative area of that class in the study site. Next, to determine whether the number of fires in each class is significantly different from random, we constructed the following Monte Carlo simulation. The 28744 fires that occurred in Sardinia during 2000–2010 were randomly reassigned to the PFCs, such that the probability of assignment of each fire to a given PFC was kept equal to its relative area. The null hypothesis is that fires occur randomly across the landscape, such that there is no difference between the relative abundance of fires in each class and its relative area within the analyzed landscape. Based on this procedure, we compared the actual value of σ for each PFC with the results of 9999 randomizations using a Microsoft Excel macro developed by Pezzatti et al. [46]. For each class, p-values (two-tailed test) were computed as the proportion of random values that were as low or lower (as high or higher) than the actual values.

Environmental characterization of the phenological classes

To frame the PFCs within the environmental context of the study area, we analyzed the degree of association of each class with the land cover types and the climatic zonation of Sardinia. The land cover map was derived from the CORINE Land Cover (CLC) 2006 data [47]. The CLC data were aggregated into seven macro-classes: urban areas, arable lands, permanent crops, heterogeneous agricultural areas, forests, shrublands, and natural grasslands and pastures. Non-combustible classes, such as bare soils, wetlands and water bodies were excluded from further analysis. Because of the relative homogeneity of the CLC (macro-)classes in terms of fuel load and continuity, these classes were considered adequate for coarse-scale analysis of fire risk (see [29]).

The climatic zonation of the Island was obtained from the climatic map of Sardinia at scale 1:250,000 in Bajocco et al. [48]. According to this map, Sardinia is classified into three main climatic regions (CRs): Mediterranean, Transitional Mediterranean, and Transitional Temperate that are mainly characterized by a clear gradient of decreasing summer drought and mean annual temperatures.

To analyze the association of fuel phenology with land cover we first constructed a two-way contingency table in which the PFCs of all fire ignition points were associated to the corresponding CLC types (see [48]). Next, the actual values of all entries of the contingency table were compared to a distribution of 9999 random values (two-tailed test) obtained under the null hypothesis of no association between the phenological fuel classes and the land cover types. For all ignition points the association between the PFCs and the CLC types was randomized keeping constant the number of fires in each phenological and land use class. In addition, we also calculated a χ2 test statistic for the entire contingency table, the significance of which was tested with the same randomization procedure. The same approach was then used for analyzing the association between the phenological classes and the climatic regions of Sardinia.

Results

Fig. 2 shows the dendrogram obtained from the cluster analysis of the phenological units of Sardinia (see also Fig. 1). Four clusters, PFC1 = (σ = 1.98), PFC2a (σ = 1.40), PFC2b (σ = 0.67) and PFC3 (σ = 0.32), were identified with selection ratios denoting decreasing fire risk from PFC1 to PFC3. According to the randomization test, fire incidence in terms of ignition is selective for all PFCs with very high significance (p < 0.001). The classes PFC1 and PFC2a showed high to moderately high fire-proneness, whereas PFC2b and PFC3 were characterized by significant fire-avoidance, meaning that the total number of fires observed in both classes is lower than expected by chance alone.

Fig 2

Dendrogram of the hierarchical cluster analysis of the mean values of the Fourier components for each phenological unit.

The dashed line indicates the cutting level of the dendrogram. Circles indicate the grouping of the phenological units (PUs) into four clusters, named phenological fuel classes (PFCs), based on their phenological similarity in terms of vegetation seasonality and productivity.

Fig. 3 shows the mean values of each Fourier component A0, A1 and A2 for each phenological fuel class. Results highlight that the additive term follows an increasing gradient from the highest to the lowest fire risk class, to the contrary, the values of the first and second harmonics progressively decrease from PFC1 to PFC3.

Fig 3

Graph of the mean values of each Fourier component for each phenological fuel class (PFC).

The first (A1), and the second (A2) Fourier component represent the vegetation seasonality trend; the additive term (A0) expresses the vegetation productivity. The y-axis is in arbitrary units.

Overall, the permutational analysis of variance showed a significant pyrological difference between the PFCs (F = 30.56, p < 0.001). Likewise, the a-posteriori pairwise comparison between PFCs revealed that all classes are statistically unique in terms of fire-proneness at p < 0.05 (Table 1). The analysis of the degree of association between the phenological fuel classes and CLC types gave an observed value of χ2 = 9423.42 (p < 0.001). The same level of significance was found for the association between the PFCs and the CRs with an observed value of χ2 = 8097.30 (Table 2 and 3), denoting that the remotely sensed phenological fuel classes proved to be strongly associated both with the corresponding land cover types and the climatic regions of the study area.

Table 1

F-values of the a-posteriori pairwise comparisons between the phenological fuel classes (PFCs) identified by the cluster analysis.

PF classes	PFC1	PFC2a	PFC2b
PFC2a	7.841*	-	-
PFC2b	27.94	6.294*	-
PFC3	100.5	36.39	14.73

*significant at p = 0.05 (9999 randomizations, one-tailed test). The remaining values are all significant at p = 0.001. The value of the test statistic for the permutational ANOVA among all groups is F = 30.56 (p < 0.001).

Table 2

Contingency table showing the number of fires shared by the phenological fuel classes (PFCs) and the Corine Land Cover (CLC) types during the period 2000–2010.

		PF classes
		PFC1	PFC2a	PFC2b	PFC3
CLC types	Urban areas	852 NS	880	353 NS	110
	Arable lands	5921	3350	754	129
	Permanent crops	521	760	216	27
	Heterogeneous agricultural areas	2598	3324	1699	673
	Forests	176	242	334	872
	Natural grasslands and pastures	241	808	386	71
	Shrublands	381	977	1186	903

Italic characters indicate a positive association between PFC and land cover.

Normal characters indicate a negative association between PFC and land cover.

NS indicates not significant at the p = 0.05 level (9999 randomizations, two-tailed test).

The remaining values are all significant at p = 0.001.

The value of χ2 for the whole contingency table is 9423.42 (p < 0.001).

Table 3

Contingency table showing the number of fires shared by the phenological fuel classes (PFCs) and the Climatic Regions (CRs) of Sardinia during the period 2000–2010.

		PF classes
		PFC1	PFC2a	PFC2b	PFC3
CRs	Mediterranean	8313	6093	1290	561
	Transitional Mediterranean	0	317	633	939
	Transitional Temperate	2377	3931	3005	1285

Italic characters indicate a positive association between PFC and climate.

Normal characters indicate a negative association between PFC and climate.

All values are significant at p = 0.001 (9999 randomizations, two-tailed test).

The value of χ2 for the whole contingency table is 8097.30 (p < 0.001).

In particular, the class of highest fire risk, PFC1, shows a marked positive association with arable lands and a negative association with the natural and semi-natural land cover types. From a climatic viewpoint, PFC1 shows a significant positive association with the Mediterranean CR, and a negative association with both transitional CRs. The class of moderately high fire risk, PFC2a, is driven mainly by fuel size rather than by climate. Accordingly, PFC2a is positively associated to urban areas and to all land use types that are characterized by fine fuels, such as permanent crops and heterogeneous agricultural areas, which are usually located in the Mediterranean portion of the study area, and to natural grassland and pastures, which are mainly located in the Transitional Temperate region. At the same time, PFC2a is negatively associated to all coarse-fuel land use types, such as shrublands and forests.

To the contrary, the moderately low fire risk class, PFC2b, has a positive association with the Transitional Mediterranean and the Transitional Temperate CRs and a negative association with the Mediterranean CR. As concerns land use types, PFC2b is positively associated with all classes of low to very low human pressure, such as heterogeneous agricultural areas, pastures, shrublands and forests, which are mainly located at higher altitudes in the transitional CRs. Finally, PFC3 shows the same climatic preferences as PFC2b. However, unlike PFC2b, this latter class shows positive association only with shrublands and forests. The environmental characteristics of each PFC are summarized in Table 4.

Table 4

Summary of the environmental characteristics of the phenological fuel classes (PFCs) obtained from the analysis.

	Number of fires	Fire risk	Fuel type	Seasonal NDVI variability	Main Land Cover types	Main Climatic types
PFC1 (4210.25 km 2)	10690	High (σ = 1.98)	Fine fuel	Very high	Arable lands	Mediterranean
PFC2a (5759.63 km 2)	10341	Moderately high (σ = 1.40)	Fine fuel	High	Urban areas, Permanent crops, Heterogeneous agricultural areas and Natural grasslands and pastures	Mediterranean Transitional Temperate
PFC2b (5743.31 km 2)	4928	Moderately low (σ = 0.67)	Coarse fuel	Low	Permanent crops, Heterogeneous agricultural areas, Natural grasslands and pastures, Forest and Shrublands	Transitional Mediterranean Transitional Temperate
PFC3 (6691.31 km 2)	2785	Low (σ = 0.32)	Coarse fuel	Very low	Forest and Shrublands	Transitional Mediterranean Transitional Temperate

Discussion

This paper follows the wake of the work of De Angelis et al. [20] and shares with it the working hypothesis and the methodological flow; the general objective of both studies, in fact, is to classify the vegetation by means of MODIS NDVI time-series. However, unlike De Angelis et al. [20], the specific aim of this paper is not to investigate the phenological metrics that best characterize the fire selectivity behavior, but to introduce a coarse-resolution satellite-based framework that can be reliably used for characterizing fire risk (i.e. the chance of a fire to start) and easily developed regardless the availability of complete wildfire database or specific fuel data.

Fire risk (or ignition risk) is defined as “the chance of a fire starting as determined by the presence and activity of any causative agent” [49]. When dealing with fire fighting strategies, fire risk is the most important variable to consider. Knowing the fire ignitions patterns and their determinants may represent a decision support tool about where the intervention forces should be concentrated for time- and cost-efficient actions. Furthermore, unlike fire risk, the fire spreading probability is highly dependent on the presence of severe weather conditions, rather than on fuel characteristics alone [48, 50–52]. Therefore, focusing on fire ignition may allow to analyze the correlation between (remotely sensed) fuel conditions and fire occurrence without the confounding effects of extreme meteorological events [49].

Understanding the geographic distribution of vegetation phenology (i.e., the timing of fuel availability) may be key to planning effective fire fighting and prevention strategies, since regions with similar fuel phenology patterns, and hence with similar fire incidence, are regions where similar management policies should be applied [53]. According to our results, NDVI seasonal timing represents a main functional driver of fire risk, in terms of both fuel load and flammability. Hence, it can be adequately used to classify the territory into fuel classes with distinctive phenological and pyrological behaviour, which in turn reflects underlying differences in land cover and climate.

From a biological viewpoint, the discriminating factors among the four PFCs were represented by the vegetation seasonal productivity (i.e. the fuel amount cumulated during the observed season) and the seasonal variability of the NDVI values, which expresses the fuel flammability in terms of transition from moist to dry vegetation. The highest fire risk classes, PFC1 and PFC2a, were characterized by a low fuel amount during the fire season and a high variability in the NDVI values. To the contrary, the lowest fire risk categories, PFC2b and PFC3, were characterized by high seasonal fuel load and low variability in the NDVI timing.

Looking at land cover and climate, the classes with the highest fire risk are mainly associated to Mediterranean climates and to urban areas, agriculture and pastures. More specifically, while PFC1 is mainly associated with urban and agricultural areas, PFC2a is associated with a more heterogeneous vegetation cover, including grasslands and pastures in more temperate areas at higher altitudes. These land cover types are generally characterized by high human pressure and fine fuel that dries quickly and can therefore burn easily after short periods of dry weather [29]. In particular, in Mediterranean regions fire has been traditionally used as a land management tool for creating pastures or eliminating agricultural waste. On the other hand, the classes with the lowest fire risk are mainly located in areas with higher rainfall and mostly covered by woody vegetation. Such coarser forest fuels take longer to dry. Therefore, they can burn only after longer (and thus infrequent) periods of dry weather. While PFC2b is mainly covered by agricultural and more natural land cover types with different degrees of human impact, PFC3 is mainly associated with shrublands and forests, showing the highest fire-avoidance due to the coexistence of reduced human impact, coarser fuels and limited summer drought.

These findings are in agreement with previous studies on fire risk in the Mediterranean basin. For instance, Bajocco & Ricotta [29] and Guglietta et al. [54] showed that the highest probability of fire occurrence in Sardinia is associated with agricultural and urban areas, due to the heaviest human presence, which represents the main source of fire ignition. In contrast, the lowest fire ignition risk can be found in the inner and higher portions of the study region, where population density is much lower. In Portugal, Catry et al. [55] emphasized the role of agriculture in driving the spatial distribution of fire events, due to the fine fuel load and the land management practices (e.g. burning crop residues). Finally, Bajocco et al. [48] in Sardinia (Italy), Pausas & Fernández-Muñoz [51] in the province of Valencia (Spain) and Oliveira et al. [56] across Mediterranean Europe, all emphasized the combined effect of land cover (which determines the fuel amount) and climate (which determines the fuel-proneness to burn) on fire occurrence.

In conclusion, MODIS NDVI time-series may represent a suitable basis for the development of coarse-scale dynamic fuel maps and fire prediction models. In addition, remotely sensed phenological data may also have additional value for various wildfire-oriented applications related to land surface monitoring and assessment (see [57]), including drought monitoring, ecological modeling of the impact of global change on burning patterns [58-59], or the identification of main areas of pre-fire intervention [18]. Furthermore, the possibility of updating the remotely-sensed fuel maps year by year with low-cost, easy to retrieve, and homogeneous data allows to keep them current so that they will retain their value for users as landscapes change.