Assessing the predictability of fire occurrence and area burned across phytoclimatic regions in Spain

Most fire protection agencies throughout the world have developed forest fire risk forecast systems, usually building upon existing fire danger indices and meteorological forecast data. In this context, the daily predictability of wildfires is of utmost importance in order to allow the fire protection agencies to issue timely fire hazard alerts. In this study, we address the predictability of daily fire occurrence using the components of the Canadian Fire Weather Index (FWI) System and related variables calculated from the latest ECMWF (European Centre for Medium Range Weather Forecasts) reanalysis, ERA-Interim. We develop daily fire occurrence models in peninsular Spain for the period 1990–2008 and, considering different minimum burned area thresholds for fire definition, assess their ability to reproduce the inter-annual fire frequency variability. We based the analysis on a phytoclimatic classification aiming the stratification of the territory into homogeneous units in terms of climatic and fuel type characteristics, allowing to test model performance under different climate/fuel conditions. We then extend the analysis in order to assess the predictability of monthly burned areas. The sensitivity of the models to the level of spatial aggregation of the data is also evaluated. Additionally, we investigate the gain in model performance with the inclusion of socioeconomic and land use/land cover (LULC) covariates in model formulation. Fire occurrence models have attained good performance in most of the phytoclimatic zones considered, being able to faithfully reproduce the inter-annual variability of fire frequency. Total area burned has exhibited some dependence on the meteorological drivers, although model performance was poor in most cases. We identified temperature and some FWI system components as the most important explanatory variables, highlighting the adequacy of the FWI system for fire occurrence prediction in the study area. The results were improved when using aggregated data across regions compared to when data were sampled at the grid-box level. The inclusion of socioeconomic and LULC covariates contributed marginally to the improvement of the models, and in most cases attained no relevant contribution to total explained variance – excepting northern Spain, where anthropogenic factors are known to be the major driver of fires. Models of monthly fire counts performed better in the case of fires larger than 0.1 ha, and for the rest of the thresholds (1, 10 and 100 ha) the daily occurrence models improved the predicted inter-annual variability, indicating the added value of daily models. Fire frequency predictions may provide a preferable basis for past fire history reconstruction, long-term monitoring and the assessment of future climate impacts on fire regimes across regions, posing several advantages over burned area as a response variable. Our results leave the door open to the development a more complex modelling framework based on daily data from numerical climate model outputs based on the FWI system.


Introduction
Fire is a global phenomenon that has a decisive influence on the ecosystems throughout the world (Bond et al., 2005;Beerling and Osborne, 2006) and as such, it must be regarded as an integral earth system process (Bowman et al., 2009).At Published by Copernicus Publications on behalf of the European Geosciences Union.
the same time, wildfires are also the cause of important damages and economic losses in many fire-prone regions of the world, arising public concerns and requiring important economic efforts towards fire prevention, protection, suppression and restoration (Hardy, 2005;Barbati et al., 2010).
It is widely accepted that fire activity is strongly influenced by climate, although human activities and land uses are also important factors in determining fire regimes (Marlon et al., 2008;de Torres Curth et al., 2012), particularly in the Mediterranean areas (Catry et al., 2010).Moreover, this relationship is not constant across landscapes, and tends to weaken in resource-limited ecosystems, in which fuel availability becomes the major constraint to fire activity (Krawchuk and Moritz, 2010;Pausas and Paula, 2012).
Even though predicting daily fire occurrence is more challenging than considering other time aggregations (e.g.weekly, monthly) that tend to reduce data variability improving model fits, from an operational point of view the daily predictability of wildfires is crucial in order to issue timely forecasts to the fire protection agencies.For instance, the European Forest Fire Information System (EFFIS) produces and publishes on the web daily maps of fire danger in Europe based on the FWI system (Camia et al., 2006;Camia and Amatulli, 2009) as the method to assess the fire danger level in a harmonized way throughout Europe.These fire danger forecasts are fed with meteorological forecasted data of similar characteristics to reanalysis outputs, although the evaluation of the FWI system for the prediction of actual fire occurrence using numerical model outputs at an operational timescale remains unexplored in many fire-prone regions of the Mediterranean (see, e.g.Viegas et al., 1999, based on meteorological records).
On the other hand, model validation is sensitive to the spatial extent on which the modelling is undertaken.The possibilities range from different levels of administrative boundaries like countries (e.g.Amatulli et al., 2013) and the EU-ROSTAT Nomenclature of Territorial Units for Statistics (NUTS3) (e.g.Carvalho et al., 2010) to other features such as ecozones (e.g.Littell et al., 2009), potential vegetation areas (Vázquez et al., 2002), hydrological basins (Pausas and Paula, 2012), etc., and largely depend on the aims and scope of each particular study.In this context, phytoclimatology is a scientific discipline focused on the establishment of links between natural vegetation and climatic types.As a result, the phytoclimatic classification provides a potentially useful spatial framework for the characterization and analysis of wildfires, identifying regions with homogeneous characteristics in terms of fuel types and climatic conditions.
In this study, we investigate the adequacy of the FWI system for the prediction of fire occurrence at a daily timescale and monthly burned area, as well as the practical applicability of numerical model output -ERA Interim (European reanalysis) -to this aim, following a previous study showing that this reanalysis product is adequate for the characterization of fire danger conditions in Iberia as compared to weather observations (Bedia et al., 2012).We also analyse the contribution of socioeconomic and land use / land cover variables in order to assess model performance and assist in the interpretation of model results.
In addition, we also analyse the suitability of monthly burned area models, potentially useful for the reconstruction of past fire history from climate records and often used for the assessment of future fire danger impacts (e.g.Amatulli et al., 2013;Balshi et al., 2009;Flannigan et al., 2005;Carvalho et al., 2010).
Furthermore, we analyse whether the phytoclimatic regions can be used as convenient generalization units for the development of accurate models of fire occurrence and burned area, by means of a spatial aggregation experiment.We determine the predictive ability of the models by applying cross-validation procedures, and calculate the contribution of the different explanatory variables to the total explained variance in order to ascertain which are the most important climatic controls of fires across phytoclimatic regions.

Fire data
We extracted fire data from the National Wildfire Database of the Spanish Environmental Agency (Mérida et al., 2007), stored at a 10 km resolution grid-cell size.We selected all fires since 1990, the year when a rigorous fire reporting protocol with a normalized form started to be applied at a national level, thus ensuring the maximum homogeneity of the fire database across Spain.

Phytoclimatic regions
The phytoclimatic regions used in this study were delimited according to the classification performed by Allué (1990) in Spain (Table 1), built upon meteorological data from the Spanish Meteorological Agency and the potential vegetation series elaborated by Rivas Martínez (1987).These phytoclimatic regions encompass a long gradient of bioclimatic conditions, ranging from the Atlantic area of influence, characterized by high precipitations and mild temperatures throughout the year, where potential vegetation is represented by broad-leaved deciduous forests, to the most arid areas of the south-east and the Ebro depression, where the natural vegetation potential corresponds to sparse formations of spiny shrubs.The resulting classification consists of 19 different subtypes of vegetation, each of them linked to characteristic climatic conditions, which are grouped in four general phytoclimatic types, then subdivided into more specific types.Due to the small area encompassed by some of the phytoclimatic zones, we made an aggregation of some of them with the neighbouring units, based on the spatial proximity and ecological affinity in terms of vegetation types.As a result, phytoclimatic types 2-3, 7-8 and 10-11-12 were merged together for the analyses.In addition, due to the very low fire numbers in the oro-boreal phytoclimatic region, the corresponding types were merged together with class 13-14-15 (Fig. 1).Similarly, type 1 was restricted to a very small coastal Mediterranean area, falling out from the 25 kmresolution land mask applied, and therefore it was discarded from the analysis.

Climate data
We obtained climate data from the ERA-Interim reanalysis, produced by the European Centre for Medium-Range Weather Forecasts in collaboration with many institutions (Dee et al., 2011).
Reanalysis products are generated from the assimilation of observational data from various sources (satellite, ships, buoys, radiosonde, surface meteorological records) through numerical simulation models in order to reproduce the state of the atmosphere with variable vertical and horizontal spatial resolution and spanning an extended historical period that covers several decades or more.The performance of different reanalysis products for FWI representation in the Iberian Peninsula and further details on their use for fire danger estimation are described in Bedia et al. (2012).The main advantages of reanalysis data are the wide geographical coverage and the homogeneity of the time series provided.In the con-text of this study, reanalysis data are an adequate surrogate for weather forecast models as they yield similar variables at a sub-daily time resolution, with the additional advantage of providing homogeneous long time series for the development and validation of daily models for an extended period for which fire records are available (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008).

Socioeconomic and LULC data
Socioeconomic and LULC data were retrieved from the "EUMED FireDrivers" database (Vilar et al., 2014), developed within the EU-funded FP7 Project FUME (Forest fire under climate, social and economic changes).It contains socioeconomic data from EUROSTAT (http://epp.eurostat.ec.europa.eu/portal/page/portal/eurostat/home/)and LULC drivers from the CORINE (coordination of information on the environment) Land Cover programme (http://www.eea.europa.eu/publications/COR0-landcover),as well as fire occurrence derived products and LULC changes data, all referred to a 10 km grid layer.The database covers the following countries of the EU-Mediterranean region: Portugal, Spain, (southern) France, Italy and Greece.The selection of the variables has been performed on the basis of their relevance for forest fire risk.An analysis of fire occurrence in 6066 municipalities in Spain of the Spanish peninsular territory and Balearic Islands revealed the strong links of agricultural landscape fragmentation, agricultural abandonment and development processes as the main anthropogenic factors (Martínez et al., 2009).Details on the variables considered in this study are provided in Table 2.

The Canadian Forest Fire Weather Index System
The Canadian Forest Fire Weather Index (FWI) System consists of six components rating the effects of fuel moisture content and wind on a daily basis, based on various factors related to potential fire behaviour (van Wagner, 1987;Stocks et al., 1989).The first three components, referred to as the Fine Fuel Moisture Code (FFMC), the Duff Moisture Code (DMC) and the Drought Code (DC), rate the average moisture content of different soil layers, respectively fine surface litter, decomposing litter, and organic layers.Wind effects are then added to FFMC to form the Initial Spread Index (ISI), which is an indicator of the rate of fire spread.The remaining two fuel moisture codes (DMC and DC) are combined to produce the Buildup Index (BUI), which rates the total amount of fuel available for combustion.BUI is finally combined with ISI to produce the Fire Weather Index (FWI), a dimensionless index rating the potential fire line intensity given the meteorological conditions at noon local standard time in a reference fuel type (mature pine stands).The FWI system uses as input four meteorological variables: daily accumulated precipitation, instantaneous wind speed, instantaneous humidity and instantaneous temperature.Further details on the procedure for FWI system calculation from reanalysis outputs are described in Bedia et al. (2012).A schematic diagram is also presented in the Supplement.

Drought indices
Drought is an important factor related to wildfire occurrence and magnitude (see, e.g.Pereira et al., 2005;Littell et al., 2009;Meyn et al., 2010).Unlike the previous daily fire danger indices, droughts are climatic phenomena difficult to quantify in terms of intensity, magnitude, duration and spatial extent, partly because there is no a straightforward manner to identify their onset, duration and end (Vicente-Serrano et al., 2010).A number of specific indices have been developed in order to quantify and properly describe drought episodes, and in this study we have included two of them: The standardized precipitation index (SPI) and the standardized precipitation evapotranspiration index (SPEI).Unlike other popular drought indices, both SPI and SPEI account for the widely accepted multi-scalar nature of droughts (see, e.g.McKee et al., 1993).Both indices are calculated on a monthly basis, and therefore they have been only tested as predictors in the burned area models.SPI is an index based on the probability of recording a given amount of precipitation at a specific point and can represent precipitation dynamics over user-selected time frames (McKee et al., 1993).However, SPI does not take into account temperature, and thus may present important devi-ations from the true water deficits derived from evapotranspiration.As a result, the recently developed SPEI (Vicente-Serrano et al., 2010) was also introduced.In practice, both indices were highly correlated and therefore only SPEI was used for burned area model building (Table 2).

Data analysis
Climate and fire data were interpolated to a regular grid of 25 km resolution, representing a compromise between the 10 km resolution of the fire and socioeconomic/LULC information and the 70 km horizontal resolution of ERA-Interim data.For each pixel, we computed daily time series of climate predictors (except SPEI and SPI, which are monthly, Table 2), fire occurrence and burned area.For the burned area models, the total monthly burned areas were calculated for each pixel and the climatic predictors were monthly averaged.
We tested two different algorithms for model development: generalized linear models (GLM, McCullagh and Nelder, 1989) and multivariate adaptive regression splines (MARS, Friedman, 1991).On the one hand, GLMs constitute a parametric method widely used, thus constituting an adequate benchmarking method.On the other hand, MARS is a non-parametric method for regression.Unlike GLMs, MARS is able to model non-linearities in the data by approximating the underlying function through a set of adaptive piecewise where the slope of each piecewise can change in a set of points called knots.The popularity of this technique is due to the efficient optimization procedure used for the iterative search for basis functions and knots.
A comparative study of both MARS and GLM in the context of binary response predictions is done in Bedia et al. (2011Bedia et al. ( , 2013)).In both types of models, the presence of redundant (highly correlated) predictors may introduce inconsistencies in variable importance estimates (see Sects.2.7.1 and 2.7.2). Thus, we first computed the pairwise-correlation matrix with all candidate explanatory variables averaged at the country level, and eliminated one of each pair attaining correlation coefficients (Spearman's ρ) greater than 0.7.We decided to preserve variable pairs below this threshold in order to avoid the loss of useful information.The resulting subset of explanatory variables is indicated by the asterisks in Table 2.

Fire occurrence models
We used generalized linear models (GLM) for fire occurrence model development, considering the logit link function for the binary response variable (fire/no fire), at a daily resolution.We selected GLMs after finding that model performance was similar than with the use of the more sophisticated MARS approach.
In order to analyse the effect of spatial aggregation of data in the models, we tested two different approaches: -Grid-box models: a full matrix of occurrence/absence of fires was constructed for each phytoclimatic zone, considering all the grid-boxes encompassed within the zone and the full daily time series.This will be referred to as the grid-box approach hereafter.
-Areal models: occurrence/absence data were aggregated at the phytoclimatic zone level, considering as occurrences all days in which at least one fire at one grid box took place, and absences those days in which no fire took place at any of the grid boxes.This will be referred to as the areal approach hereafter.
In order to test the sensitivity of the models to fire size, we set different burned area thresholds for occurrence definition: 0.1, 1, 10 and 100 ha.As a result, only fires above the corresponding area thresholds were computed as occurrences.
Fire occurrence models were trained using all occurrence samples and fire absences randomly chosen in an equal number to the fire occurrences, thus using balanced data sets for model training to avoid an artificial inflation of model skill (see, e.g.: Manel et al., 2001;McPherson et al., 2004).For the grid-box model training, fire absences were sampled only from those days in which no fires occurred in any of the grid-boxes of the phytoclimatic zone and this process was repeated 100 times in order to get a confidence interval of the sampling error.We then tested the resulting models using a random sample containing all possible cases.We undertook a one-year out cross-validation procedure, using 18 yr for training and the remaining one for testing, repeating this process 19 times, exactly one per year.In the case of socioeconomic/LULC variables, we performed a 3-fold cross-validation, each fold corresponding to the periods at which these statistics are given (Table 2), in order to avoid model overfits due to the repeated values.
From the resulting probabilistic predictions we computed the area under the receiver-operating characteristic curve (ROC score area, RSA hereafter), which provides a quantitative measure of model performance (Swets, 1988).The possible range of RSA is (0,1).A null performance is indicated by RSA = 0.5, when the ROC lies along the positive diagonal, whereas RSA = 1.0 corresponds to a perfect performance.A RSA value below 0.5 corresponds to a ROC curve below the diagonal, indicating the same level of discrimination capacity as if it were reflected about the diagonal, but wrongly calibrated (Jolliffe and Stephenson, 2003).In order to compute the predicted fire frequencies, the probabilistic model predictions were converted into a binary prediction using two different approaches for decision threshold determination.Further details on the modelling approach are provided in the Supplement.
-A global fixed probability threshold was determined by calculating the likelihood ratio of fire occurrence as given by the observations, and then applied to the full vector of predictions.This approach will be termed as global threshold hereafter.
-A monthly-varying threshold, corresponding to the likelihood ratio of fire occurrence as given by the observation, conditioned to the month.As a result, 12 different decision thresholds were obtained for each month of the year, which were applied to the predictions of the corresponding months to obtain the binary prediction.This threshold is referred to as monthly threshold in the following.
In order to estimate variable importance in the context of GLMs, we applied the method of hierarchical partitioning, by which the independent effect of each variable is calculated by comparing the fit of all models containing a particular variable to the fit of all nested models lacking that variable (Chevan and Sutherland, 1991).For instance, for variable X 1 , its importance I would be calculated as follows: where X h is any subset of i predictors from which X 1 is excluded.As a result, the variance shared by two or more correlated predictors can be partitioned into the variance attributable to each predictor.This method provides a robust assessment of variable importance and has been shown to outperform other methods used for variable importance estimation in the context of regression analysis (Murray and Conner, 2009).Finally, we also fitted multiple linear regression models to the monthly fire counts in order to test whether the daily occurrence models do provide or not an added value to the predictability of the inter-annual fire frequency variability with respect to the monthly ones.
In order to obtain robust estimates of model performance, we carried out a leave-one-out cross-validation procedure (LOOCV) to compute the error (Michaelsen, 1987).LOOCV is a resampling technique in which n−1 instances out of the total of n are used as the training data set and the remaining one is used for testing.The procedure is repeated n times, one per observed instance, producing a more precise estimation of the classification accuracy.The method assumes that each sample is independent, so prior to its application we constructed autocorrelation plots of the monthly burned areas.We found a slight autocorrelation (maximum of 0.26) at some phytoclimatic types that were significant at the α = 0.05 level.However, time series with autocorrelation of 0.25 or less will have an effective sample size at least of 90 % of the original sample size (Michaelsen, 1987), and thus it can be considered that this does not produce any measurable effect on the LOOCV estimates.
For variable importance estimation in the context of MARS, we looked at the reductions in the generalized crossvalidation estimate of error (GCV) in the selection routine performed by the MARS algorithm (Milborrow, 2013b).The GCV is reduced each time a new variable is entered into the model.The accumulated reductions in GCV can be used as an estimate of variable importance, a value that is scaled to have a maximum of 100 and a minimum of zero (this mini-mum is reached when the variable is not used at all, or produces positive changes in GCV).
All the analyses were conducted in the R language and environment for statistical computing (R Core Team, 2013).The hierarchical partitioning was undertaken using the R package hier.part(Walsh and Mac Nally, 2013).For the MARS models, we used the implementation of the algorithm included in the R earth package (Milborrow, 2013a).The drought indices SPI and SPEI were computed using the SPEI R package (Beguería and Vicente-Serrano, 2013).

Fire regimes of the different phytoclimatic regions
We found two contrasting fire regimes in terms of area burned and number of fires across phytoclimatic zones: one characterized by a bimodal annual pattern, and another one exhibiting an unimodal annual cycle, with the fire season concentrated in the summer months (Fig. 3).The first case corresponds to the phytoclimatic types under the Atlantic influence (10-11-12 and 13-14-15), with two marked peaks of fire activity in March and August.These regions are characterized by temperate and wet conditions during most of the year, and also by relatively low fire danger conditions.In spite of the less suitable conditions for fire activity of these regions, there is a high number of fire records and also large burned areas.In this regard, previous studies have highlighted the strong influence that humans exert on fire regimes, and how fire incidence can be greatly enhanced for this reason even when climate conditions are not the most favourable (Vázquez et al., 2002).On the other hand, in the remaining phytoclimatic zones, which belong to the Mediterranean climate, fire activity is concentrated in the summer months and exhibits a marked unimodal annual cycle, coincident with the most favourable climatological conditions for fire.
Thus, the bimodalities in fire occurrence can be attributed essentially to anthropogenic factors, as depicted in Fig. 3, where the largest contributions to total explained variance are in most cases done by climatic variables, and only in the bimodal ones (and especially in zone 13-14-15), the contribution of some socioeconomic/LULC variables is more important.This is more accentuated in the case of the grid-box approach, which operates at a lower spatial scale, whereas in the case of the areal models their contribution is always overrun by climatic variables, highlighting the finer spatial scale at which these variables become relevant.The strong anthropogenic influence of fire regimes in this area of Spain (and especially in the north-western corner) is a phenomenon well described by previous authors (see, e.g.Martínez et al., 2009).For the rest of zones, variable importance assessment revealed that FWI and temperature are the chief climatic controls of fire occurrence, followed by other drought-related components of the FWI system (DC, FFMC).For the sake of conciseness, the results for the occurrence models are illustrated for the 10 ha burned area threshold models, although very similar results were obtained for the remaining burned area thresholds.Due to the marginal contribution of socioeconomic/LULC variables to model performance, in the following the results presented correspond to the climate-only models.

Fire occurrence models
Most phytoclimatic zones attained moderate to good model performance, and only zones 2-3, 10-11-12 and 13-14-15 yielded RSA values below 0.70 in the case of the grid-box models, whereas the areal models yielded higher RSA values in most cases for the larger area thresholds (Table 3).The inclusion of socioeconomic/LULC variables in the models did not contribute at all to model performance in the case of areal models, and only marginally in the case of the grid-box models, even in the case of phytoclimatic zone 13-14-15, in which these variables attained higher importance than the climatic ones.In consequence, for brevity the corresponding RSA values using socioeconomic/LULC variables are not shown.In general, all models attained higher skills with increasing burned area thresholds, showing that the fire weather predictors used are more sensitive to the detection of larger fires than to smaller ones, being the latter not so closely dependent on favourable climate conditions for their occurrence.For the climate-only models, the RSA scores attained were similar when considering the 3-fold and the leave-oneyear-out procedures.
The sampling error related to the random selection of days without fires was very low, and increased slightly for higher area thresholds, although they were always in the range of hundredths of RSA (not shown).
In order to obtain the deterministic binary predictions of occurrence/absence of fires from the probabilistic output of the models, we applied both the global and monthly probability thresholds as cutoff value.We did not find a clear advantage of using the monthly probability threshold as compared to the global one in terms of correlations between observed and predicted series, and therefore in the following the results from the global threshold are presented (see the Supplement for extended information).
We found that the good model skills are not directly linked to a good reproducibility of observed inter-annual fire frequencies when working at the grid-box scale.In this regard, all models tended to a large overestimation of fire occurrence at this spatial scale because all events of high danger potential are given a high probability of occurrence, in close relationship with the annual cycle of fire weather danger.However, this effect is overridden when considering a larger spatial aggregation unit because the probability of having at least one fire in a larger region when the conditions are favourable is much higher (this is illustrated in more detail in the Sup-plement).Thus, the use of the grid-box approach provides an adequate basis for model assessment, as the good performance of the models at the grid-box scale translates into a good reproducibility inter-annual of fire frequencies across most phytoclimatic zones within the areal context (Table 4).The results of the areal models will be presented in the following (Fig. 4).
In spite of the higher RSA attained by the models for the 100 ha area threshold, their ability to reproduce the interannual fire frequencies was poor, contrasting with the good reproducibility of fire frequencies for lower area thresholds.In this sense, there are two different issues jointly affecting the reproducibility of inter-annual fire frequencies: on the one hand, this is affected by model performance (RSA).Obviously, those models with low RSA have less ability to adequately reproduce the inter-annual fire frequencies.This can be seen in the examples of zones attaining low RSA values (e.g.13-14-15), especially with the area threshold of 0.1 ha (Table 3).However, it is also important to take into account the inter-annual variability of fire occurrence, which in turn is related with the prevalence of the phenomenon.For instance, most models attain relatively high RSA values for the 100 ha area threshold, but the inter-annual variability of these large fires is lower because there are few fires larger than 100 ha every year (even zero some years at some regions; see Fig. 4 and also the Supplement for fire count data).As a result, the estimated inter-annual frequencies are more drastically affected by false negatives/positives in the case of large fires than in the case of small ones, even though in both cases RSA values are good.In Fig. 4 it can be seen how the best inter-annual correlations correspond to those zones where the inter-annual variability is higher.This variability tends to reduce as the fire area threshold increases, leading to worse results.
Besides, in general the monthly fire count models captured better the inter-annual variability of fire frequency when considering the 0.1 ha burned area threshold (with the notable exception of the phytoclimatic zone 13-14-15, Table 4), although daily models improved this correlation for higher area thresholds, suggesting that the loss of information when considering monthly aggregated statistics may be harmful for an optimal reproducibility of this events due to their dependence on very particular weather conditions that require a finer temporal scale to be represented.

Burned area models
The predictability of burned area was variable across zones, although in general model performance was poor and in most zones we obtained non-significant inter-annual burned area cross-correlations between observed and predicted time series.The explanatory variables provided an added explained variation when compared against a simple model incorporating only the annual cycle (not shown), and therefore the skill of the models cannot be solely attributed to the seasonal Gridbox 13-14-15 Fig. 3. Variable importance (% of total explained variance) in the fire occurrence models (10 ha area threshold) for each phytoclimatic zone and considering both the grid-box and areal models and the inclusion of socioeconomic/LULC as co-variables (pink colour).The results presented correspond to the 3-fold-cross validation, each fold corresponding to a different period of the socioeconomic LULC statistics (Table 2, Sect.2.7.1). cycle of the fire danger predictors.Nevertheless, the inclusion of the socioeconomic/LULC variables provided none or little improvement to model performance (Fig. 5), and therefore model performance, when any, can be attributed solely to the climate variables.
Only in phytoclimatic zones 6 and 10-11-12 a fair reproducibility of inter-annual burned area series was attained.In both cases, temperature, FWI, DC and FFMC were important explanatory variables (Fig. 6).In addition, the drought indicator SPEI was also a relevant predictor in the case of phytoclimatic zone 10-11-12.Several studies in the Mediterranean have highlighted the link of burned area with an-tecedent conditions (see, e.g.Pausas, 2004;Turco et al., 2013;Koutsias et al., 2013), particularly lagged precipitation, which determines biomass production and accumulation during the growing season, and hence fuel availability during the summer.Models presented in this study do not incorporate lagged variables, and only some of them (DC, SPEI) are more dependent on antecedent conditions.The lack of these variables may partly explain the poor predictive performance of the burned area models.In this regard, Carvalho et al. (2011) indicate that DC estimated in summer still reflects springtime atmospheric conditions, justified by the slow reacting character of this indicator that models variations on Fig. 4. Observed and 1 yr-out cross-validated predicted annual fire frequencies at each phytoclimatic zone, considering the climatic predictor data set and the four different burned area thresholds tested according to the areal models.Cross-correlation values between observed and predicted series are presented in Table 4.
moisture of deep organic soil layers.Other studies conducted in Portugal show that the DC gives a good indication of wildfire behaviour and propagation and also of the relative hazardousness of a fire season due to its long-term response to daily weather variations (Viegas et al., 2004).Inter-annual correlations between observed and modelled fire frequencies (Fig. 4, Table 4) are much higher than those obtained with burned area (Fig. 5), showing that fire occurrence is better modelled from climate data alone than burned area.Within the framework of future climate impact assessment, the projections of future fire danger scenarios are most often based on the simulation output of GCMs (either downscaled or not) run in transient mode (e.g.Bedia et al., 2013).This implies that model predictions do not have a day-today correspondence with real climate, and their value lies in the ability of the models to represent multi-decadal climatic features (mean state, variability, trend, etc.).This suggests that the estimation of inter-annual frequencies from simulated model outputs for sufficiently long time slices (typically 30 yr periods) is able to provide a more robust estimation of future fire impacts than area burned, the latter often yielding too inflated, unrealistic future estimations, as shown in several studies (see, e.g.Amatulli et al., 2013;Balshi et al., 2009;Flannigan et al., 2005;Carvalho et al., 2010).

Conclusions
Our results support the use of ERA-Interim reanalysis and the FWI system for fire modelling applications in Spain, provided an adequate temporal and spatial scale of data analysis.The added value of using daily occurrence data vs. monthly fire counts for modelling inter-annual fire frequency variability was demonstrated for fires of medium to large sizes, and even for all fire area thresholds in some particular zones.Grid-box models of fire occurrence yielded in general good model skill in terms of RSA, although this fact does not translate directly into a good reproducibility of fire frequencies due to the inherent tendency of the method to overestimate fire occurrence.Nonetheless, the annual cycle was adequately modelled regardless of the distributional characteristics of the different annual fire regimes of each phytoclimatic zone.Temperature and some components of the FWI system (FWI itself and DC and FFMC in particular) were the most important predictors of fire occurrence.
The inclusion of socioeconomic/LULC covariates did not significantly contribute to the improvement of model performance, neither for the fire occurrence models nor for the burned area ones.
Areal models yielded accurate predictions of the interannual fire frequency series, showing that the aggregation of the data into larger spatial units is needed for an adequate analysis the climatic drivers of fires, and that the phytoclimatic zones used in this study constitute representative units of the climate-fire relationship, although other convenient aggregation units may be also adequate.In addition, the good reproducibility of inter-annual fire frequencies using the fire occurrence models developed may prove useful for fire history reconstruction in fire-prone areas, for instance in order to obtain long fire frequency annual series extending before satellite data -often used to this aim.Most of burned area models failed to adequately reproduce the inter-annual burned area series, with some exceptions.In these cases, burned areas were mostly explained by temperature and drought-related indices bearing some sort of "memory" on the antecedent conditions, such as the FWI system component DC or the recently developed drought index SPEI.As a result, the practical application of burned area models, building solely on the FWI system for the prediction of future burned area scenarios, poses important limitations, being the fire frequency models are more robust for this aim.

Fig. 2 .
Fig. 2. Mean monthly area burned (vertical bars, main y axis) and mean number of fires (lines and dots, secondary y axis) recorded for the period 1990-2008 in Spain, at the different phytoclimatic regions.

Fig. 6 .
Fig.6.Variable importance of climatic and socioeconomic/LULC variables of the burned area models for phytoclimatic zones 6 and 10-11-12 (the only ones attaining significant inter-annual correlations of observed and modelled fire frequencies).Values for each variable correspond to the scores obtained by each of the 3-fold cross-validated models according to the three periods of socioeconomic/LULC data (see Sect. 2.7.1 for details).

Table 1 .
Summary of the phytoclimatic regions considered for spatial aggregation in this study.The spatial distribution of the phytoclimatic types is displayed in Fig.1.The area occupied by each phytoclimatic type is also indicated.

Table 2 .
Summary of predictors tested in this study for fire model building.After checking for redundancy, a final subset of weakly-correlated variables was used for model building, marked with an asterisk (Note that SPEI is a monthly indicator, and therefore it was only used as a predictor for the burned area models).In the case of the socioeconomic and LULC variables, data correspond to three periods:1990-1999,  2000-2005 and 2006-2012.Note that the relative areas occupied by communication networks have been rejected because of the coarse spatial resolution of the original data (NUTS2 level).ET0 = potential evapotranspiration.www.nat-hazards-earth-syst-sci.net/14/53/2014/Nat.Hazards Earth Syst.Sci., 14, 53-66, 2014 linear regressions known as basis functions of the form:

Table 3 .
ROC skill area (RSA) attained by the fire occurrence models for each phytoclimatic zone.Results are presented for the different burned area thresholds used for fire occurrence definition, and considering the leave-one-out cross-validation approach (see Sect. 2.7.1 for details).

Table 4 .
Observed and cross-validated predicted annual burned areas at each phytoclimatic zone, considering the climatic predictor data set.Cross-correlation values (Spearman's rho) between observed and predicted series are indicated by the figures within each panel.Spearman's rho cross-correlation coefficients between the observed and predicted annual fire frequencies determined by the daily (D) and monthly (M) occurrence models, determined using the global probability threshold for the areal models.The number of grid boxes comprising each phytoclimatic zone are indicated by the N grid-boxes column.The results are presented for the different burned area thresholds used to define fire occurrence/absence.