Wildfires and the role of their drivers are changing over time in a large rural area of west-central Spain

During the last decades, wildfires have been changing in many areas across the world, due to changes in climate, landscapes and socioeconomic drivers. However, how the role of these drivers changed over time has been little explored. Here, we assessed, in a spatially and temporally explicit way, the changing role of biophysical and human-related factors on wildfires in a rural area in west-central Spain from 1979 to 2008. Longitudinal Negative Binomial (NB) and Zero-Inflated Negative Binomial (ZINB) mixed models, with time as interacting factor (continuous and categorical), were used to model the number of fires of increasing size (≥1–10 ha, >10–100 ha, >100 ha) per 10 × 10 km cell per year, based on fire statistics. We found that the landscape was rather dynamic, and generally became more hazardous over time. Small fires increased and spread over the landscape with time, with medium and large fires being stable or decreasing. NB models were best for modelling small fires, while ZINB for medium and large; models including time as a categorical factor performed the best. Best models were associated to topography, land-use/land cover (LULC) types and the changes they underwent, as well as agrarian characteristics. Climate variables, forest interfaces, and other socioeconomic variables played a minor role. Wildfires were initially more frequent in rugged topography, conifer forests, shrublands and cells undergoing changes in LULC types of hazardous nature, for all fire sizes. As time went by, wildfires lost the links with the initial fire-prone areas, and as they spread, became more associated to lower elevation areas, with higher solar radiation, herbaceous crops, and large size farms. Thus, the role of the fire drivers changed over time; some decreased their explaining power, while others increased. These changes with time in the total number of fires, in their spatial pattern and in the controlling drivers limit the ability to predict future fires.

have attempted to assess the ongoing evolution of the roles that biophysical conditions and human factors play in fire occurrence (presence/absence) and frequency (number of fires per unit area and time) [16][17][18][19][20] from a stationary temporal point of view. In addition, other studies at EUMed have applied different models for various time periods and later comparing them over the longer term [21][22][23][24][25] . However, models that do not assume stationarity in the controlling variables (e.g., panel or longitudinal models) can be most useful [26][27][28][29][30] . On the other hand, to account for spatial heterogeneity (spatial non-stationarity), local Geographically Weighted Regression (GWR) has been commonly applied 18,20,23 . Nevertheless, to explicitly model spatial heterogeneity, hierarchical, multilevel modeling or the General Linear Mixed Models (GLMM) are also appropriate. These models have been successfully applied for predicting fire occurrence [31][32][33] . Finally, for dealing with high number of zero observations, zero-inflated models have been used for predicting fire spatial patterns 4,[34][35][36][37][38] .
In this work, we modeled the variation of wildfires in space and over time by applying longitudinal Negative Binomial (NB) and Zero-Inflated Negative Binomial (ZINB) mixed models over a thirty year period, using wildfires of three size categories (≥1-10 ha, >10-100 ha, >100 ha) that occurred in a large rural area in west-central Spain. The objective was to determine the relationship between occurrence and frequency of wildfires as a function of various factors, including topography, LULC types and their changes, socio-economic variables, forest interfaces, linear features and climate in a dynamical way. We expected that wildfires would have increased in the area and during the study period due to changes in some of the main fire drivers: i) an increase in landscape hazardousness as result of LULC changes (i.e., land abandonment, afforestation, etc.), ii) socio-economic restructuring in the area due to rural exodus (i.e., population decline, shifts in economic activities), and iii) changes in climatic conditions (i.e., warming).

Methodology
The Study area. The study area covers 56,000 km 2 in west-central Spain (UTM coordinates 4369-4551 and 201-394 30 N) (Fig. 1A). The area is characterized by the mountainous landscapes of Sierra de Gredos, and the gentler mountains at the SE, both flanked by relatively flat areas (Fig. 1B). The climate is Mediterranean, with colder and wetter conditions up in the mountainous areas, and warmer and drier at the low-lands (Fig. 1C). Soils in the mountain areas are shallow, with high stoniness and coarse texture (cambisol, regosol and lithosol), whereas, in the low-lands, they are deep and fine textured (luvisol and fluvisol) (Fig. 1D).
Fire data and explanatory variables. We used the yearly number of fires recorded at 10 × 10 km grid cells (n = 443) of the Spanish Fire Statistics database ( Table 1). These variables proved to be important in explaining wildfires in other areas 12,13,18,19 . An exploratory analysis based on Principal Components (PCA) was carried out to assess collinearity among explanatory variables (i.e., co-variates). Moreover, trends over time of the number of fires per cell for each fire size category and of each explanatory variable were assessed by means of the Mann-Kendall test (Kendall package 39 in R software 40 ). Additionally, we assessed the significance of the proportion of cells that changed over the years against the null hypothesis of an even proportion of positive/negative or no change, by a G test of goodness of fit 41 . The datasets generated during the current study are available in the figshare repository https://doi.org/10.6084/m9.figshare.7159727.

Statistical approach.
Longitudinal NB and ZINB mixed models. In gridded fire data, like the ones used here, fire counts are often characterized by a high number of zero observations, as no fires may occur in some years in a high number of spatial units. The excess of zeros yields a variance greater than the mean (over-dispersion) 42 . To handle over-dispersion in count data, a Negative Binomial (NB) distribution may be appropriate. Moreover, when in addition to over-dispersion an "excess of zeros" may be present, Zero-Inflated Negative Binomial (ZINB) models are theoretically most suitable 43 . ZINB models are composed of two parts: 1) a logistic regression for the "excess zeros" over and above what would be predicted by the count distribution, and 2) a regression (either Poisson or Negative Binomial) for the counts part 42,44 . Likewise, the number of fires per year at a given cell (i.e., fire frequency) is a type of longitudinal data with repeated measures, in which the observations are taken over time on the same subject (i.e., 10 × 10 km grid cells), thus observation are not independent. In these cases, longitudinal NB or ZINB mixed models handle dependency by incorporating an appropriate random component in the model structure [45][46][47] . To evaluate suitability of one modelling approach over the other the Vuong statistic 48 was used. However, we did not find strong evidence to conclude that one model fitted the data better than the other 49 , and thus we decided to test both NB and ZINB models.
Modelling strategy. We applied longitudinal NB and ZINB mixed models to three fire size categories (small [≥1-10 ha], medium [>10-100 ha], and large [>100 ha]) in order to explain the number of fires per cell per year as a function of several co-variates fitting univariate models (Fig. 2). To account for the spatial and temporal dependencies of our data, two random components were included: (i) time (T, years) as level 1 (repeated measures), which represents the "within-subject variance", and (ii) an identifier for each grid cell as level 2 (subjects), which represents "between-subject variance" 45  and co-variates, time (T) was additionally included as interacting factor. This was done in three different ways, either it was not considered (T = 0), or considered as a continuous (T = 1, 2, …, 30) or categorical (T = 1; T = 2; …, T = 30) factor. This permitted obtaining three types of models, respectively: (i) "Average" models, which provided the average response and the average effects of co-variates over the whole period; (ii) "Average change" models, which provided the initial response and the initial effect of each co-variate, as well as the average rate of change (positive or negative) over the years; and (iii) "Annual change" models, similar to the previous models except that they provided coefficients of change between the initial year and any other given year separately. "Annual change" models permitted a finer analysis of the trends, which were later explored by means of a Mann-Kendall test (Fig. 2). Null models (without covariates, thus being pure random effects models) provided the baseline to assess the performance of univariate models for each co-variate. Model outputs were: intercept, slope of the time factor, and logit coefficient of the co-variates. The logit-coefficients, which were standardized by z-score and were conditional on the random effects, represent the following: in the "Average" models, the average effect of covariate (X) over the response variable (Y) when X changes between subjects and across time by one standard deviation (SD) unit 47 . In "Average change" and "Annual change" models, the regression coefficients of X were also conditional terms and must be interpreted along with the interaction term (X*Z [time interaction term]) 46 . Finally, the logit coefficients were exponentiated to obtain the Odds Ratio (OR) and the Rate Ratio (RR). The OR represents the increasing or decreasing odds of the excess of zeros in the outcome (Y) for the zeros part in ZINB mixed models, and the RR represents the percentage increase or decrease in the outcome (Y) in NB mixed models and in the count part of ZINB mixed models. Fire occurrence (presence/absence) was assessed by the intercepts and coefficients of the zeros part in ZINB mixed models and fire frequency (number of fires) by the intercepts and coefficients of NB mixed models and those of the count part in ZINB mixed models.
Model evaluation and goodness of fit. Longitudinal NB and ZINB mixed models were calculated based on a Bayesian approach, using strong normal priors (V = diag (2) and nu = 4) after a sensitivity analysis. The posterior distribution of the model parameters were obtained using the Markov Chain Monte Carlo (MCMC) method, whose accuracy was assessed by density plots, autocorrelation analysis and the Geweke diagnostic 50 . Furthermore, we assessed the accuracy of the longitudinal NB and ZINB mixed models by linear regressions between the observed versus predicted values, using annual data (n = 13,290), proportion of zeros and counts correctly estimated, and by spatial overlays between predicted and observed values using aggregated (whole period) data. All models were calculated using the MCMCglmm package 50 in R software 40 .
Given the large number of models, we focused on those univariate NB and ZINB mixed models that had a reduction [<−10] of DIC (Deviance Information Criterion), showed significant regression coefficients of the co-variates (p < 0.05), and explained >20% of random variance relative to that of null models (pseudo-R 2 ). For each model, we checked that Markov chains did not show autocorrelation or any patterning in parameter estimations.

Results
Changes in wildfires and drivers. During the 30-yr period, there were 23,435 fires, which burned a total of 456,253 ha. Number of fires per cell over the study period varied between 0 and 326. Annually, several cells contained a high proportion of zeros, thus the distribution was heavily left skewed (74%, 88% and 96% of total annual counts were zeros for fires ≥1-10 ha, 10-100 ha and >100 ha, respectively). The three fire size categories accounted for 34%, 11%, 3% of all fires, and 7%, 18%, and 74% of the burned area, respectively (Fig. 3A). In the whole area, temporal trends of small fires showed a significant increase (Mann-Kendall Tau: +0.54, p < 0.001), no significant trend was obtained for medium fires, and a significant decreasing trend for large fires (Tau: −0.32, p < 0.01) (Fig. 3B). Temporal trends at cell level indicated that 26%, 4% and <1% of cells showed a significant positive trend, whereas 3.8%, 8% and 7% of cells showed a significant negative trend for each fire size, respectively (Fig. 3C). PCA of the explanatory variables showed that these could be merged in 3 main axes (SM Fig. 1). Axis 1 (38% of total variance) was positively related to areas with high solar radiation and high density of livestock, and negatively related to high elevation and slope ranges, with hazardous vegetation (shrublands and dense forests), which experienced hazardous stability and hazardous LULC changes (i.e., densification and afforestation). Axis 2 (24%) was positively related to areas with herbaceous crops, crop-trees and WAIs, and negatively related to pastures, agroforestry, large farms, agriculture abandonment and WGIs. Axis 3 (19%) was mainly related to artificial uses, development changes (i.e., artificialization), population density, unemployment and WUIs.
Moreover, we found that landscapes in the study area were rather dynamic (SM Fig. 2 and SM Table 2). In general, landscapes became more hazardous over time due to the increase of hazardous stability (i.e., hazardous LULC types that remained so over time) thus accumulating fuel, the spread of open and broadleaved forests as well as pastures, which resulted in important spatial changes in forest interfaces. Other important changes were the appearance of development in former more natural areas, hence increasing WUIs; the spread of non-hazardous (i.e., herbaceous crops and agroforestry) over areas in which they were not so abundant, thus raising the WAIs; and the loss of abundance of some hazardous LULC types (i.e., coniferous forests, shrublands) in some cells due to wildfires 10 , which was compensated by increases in other cells not traditionally occupied by such LULC types at the beginning of the study. Regarding socioeconomic factors, we found: population increases in peri-urban areas around the main towns and cities; a reduction of employees in the primary sector, with an increase employment in the building sector; a significant reduction of farm density, together with a decrease of small farms and an increase of larger ones; a significant increase of livestock density and a reduction of agricultural machinery. Finally, in relation to climate, there was not any clear trend for any of the variables analyzed (SM Fig. 2).
Wildfires modelling. The number of fires per cell per year was small at the beginning (1979) for all fire size categories, particularly for large fires, as is reflected in the low rate ratio of the intercepts of "Average change" null NB and ZINB mixed models. According to NB mixed models, small fires tended to increase over time, medium fires did not show any significant trend, and large fires decreased, as shown by the effect of time (Table 1). Similarly, time effects in the "Average change" null ZINB mixed models indicated that the probability of small and medium fire occurrence (presence/absence) increased over time (lower excess of zeros), whereas the frequency (number of fires per cell per year) of medium and large fires (counts part of the model) slightly decreased (Table 1, Fig. 4). Based on reduction of DIC, the best univariate models were the NB mixed models for small fires, and the ZINB mixed models for medium and, particularly, for large fires (Fig. 5A, Table 2). In relation to the way "time" was included in the models, "Annual change" models were always better than the other two models for all fire sizes, although for large fires, "Average" and "Average change" models were also appropriate (Fig. 5B). Linear regression between observed versus predicted values indicated that, on average, for all co-variates and significant models, small fires were better explained (R 2 = 62%), followed by medium (R 2 = 35%), and, to a lesser extent, by large ones (R 2 = 18%) (Fig. 6A). Moreover, from ZINB mixed models we ascertained that the number of zeros was better predicted for both medium and large fires than for small ones (91% 103% and 69%, respectively), whereas the number of fires (counts) was better predicted for small fires than for medium and large fires (89%, 57%, 19%, respectively). Spatially, there was an error of ±15 fires per cell for small fires, −9/+2 for medium fires and −10/+5 for large fires (Fig. 6B-D). The role of fire drivers on wildfires. Some of the co-variates were able to largely reduce DIC and explain a large part of the random variance (pseudo-R 2 ) (SM Tables 3-8). In general, topography, LULC types, LULC changes and agrarian characteristics were the group of variables that were more related to fires of all sizes, although with differences depending on fire size ( Table 2, SM Tables 3-8). We found that topographic features were strongly linked to fire occurrence and frequency for all fire sizes, particularly for large fires. High fire occurrence (low excess of zeros) and frequency were linked to high elevation and slope ranges, whereas in low elevation areas with high radiation occurred the opposite. Regarding LULC types, fires of all sizes were positively related to conifer forests, which mainly explained fire occurrence of small fires, and fire frequency of medium and large ones. Similarly, shrublands were positively related to the occurrence of all fire sizes, particularly for small and medium fires. On the contrary, herbaceous crops were negatively linked to fires of all sizes, mainly large fires. Regarding to LULC changes, hazardous stability and afforestation explained positively fires of all sizes, whereas agriculture conversion showed the opposite. Concerning socio-economic variables, agrarian characteristics showed that large farms, usually under leasing, were negatively related to fires of all sizes. Population and employment variables were weakly related to fires. Finally, forest interfaces, linear structures and climate variables were marginal in explaining fires ( Table 2, SM Tables 3-8).
Changes through time in the role of fire drivers. The changing role through time of fire drivers on wildfires was ascertained using Mann Kendall test on the logit coefficients of "Annual change" NB and ZINB mixed models (Table 2). Regarding topography (SM Fig. 7), areas with high elevation and slope ranges showed lesser probability of having fires (greater number of zeros) over time. In contrast, areas with high radiation (i.e., flatter areas) showed the opposite trend. Regarding LULC types (SM Fig. 8), coniferous forests tended to loose statistical weight (coefficients in zeros part tended to zero) for explaining fire occurrence for all fire sizes, and showed negative trends for fire frequency regardless of size. Additionally, shrublands tended to have less fire occurrence of small fires over time, not showing significant trend for medium and large ones. Conversely, herbaceous crops showed a clear trend to increase the probability of occurrence of small fires, and reduced their statistical weight for explaining fire frequency of fires of all sizes (coefficients tended to zero in the count part). More to, areas with hazardous stability and afforestation showed less probability of fire occurrence and less fire frequency, respectively, over time for small and medium fires. On the contrary, agriculture conversion showed a positive trend for the frequency of small and medium fires. Similarly, large farms reversed its sign, increasing fire frequency with time (SM Fig. 9).

Discussion
This study documents that the studied landscape of west-central Spain was rather dynamic. Wildfires were also very dynamic, as they changed in time and space, changes which were size dependent: small fires increased and spread over the landscape, while medium and large were stable or decreased. Models used to explain wildfires were also size dependent, with small fires being better explained by NB and larger fires by ZINB. Time was an important factor, and including it in the models as a categorical variable improved model performance. Not all sets of variables were able to equally explain wildfires, with topography, LULC types and their changes, and agrarian characteristics producing the best models. Within these sets, several variables had similar model performance, suggesting high collinearity between them. Additionally, the role of variables in explaining fires changed along the time, so that variables that were important at the beginning of the study period were not so at the end. Following, we discuss these findings in detail. We found that the landscape in west-central Spain was rather dynamic. The main changes were characterized by the maintenance of hazardous LULC types, an increment in open forest and grasslands, as well as spatial changes in forest interfaces. The observed LULC changes were compatible with an increase in landscape fire hazard. For example, the increment in open forests and grasslands likely contributed to enhance fuel continuity. Furthermore, the increase in forest interfaces, notably WAIs, elevated the accessibility of ignitions to forest and other flammable areas. These processes were similar to those described in other Mediterranean countries 7,10,51 . Such changes have often been related to increases in wildfires 52,53 .
Wildfires were also dynamic, but this was fire size dependent. Whereas small fires increased over time and became widespread much over the whole landscape, medium fires remained stable and large fires decreased. Medium and large fires exhibited a spatial concentration, since they tended to disappear from areas in which they were present at the beginning of the study period, without spreading to new areas. The fact that the number small fires tended to increase suggests that ignitions were facilitated by a landscape with greater number of interfaces. Despite of that, and the more hazardous landscape that has evolved with time, large fires tended to decrease, which supports that firefighting services are rather efficient at deterring fires 9 . Notwithstanding, when conditions are severe, very large, even megafires, are erupting in the region 54 , supporting that there is limit to what the firefighting services can accomplish when conditions pass certain threshold, during which extant landscape can't but facilitate fire spread. Our findings are compatible with wildfire trends across southern European Mediterranean-type countries. Indeed, wildfires have been decreasing in these countries 3,6,55 . This tendency, however, varies among countries and within countries. For example, while the reduction of fires and area burned in Spain, France and Italy are clear, that is not so much the case with Portugal 54,55 . Moreover, within Spain, while in Figure 5. Number of co-variates with the highest DIC (Deviance Information Criterion) reduction in: longitudinal Negative Binomial (NB) and Zero-Inflated Negative Binomial (ZINB) mixed models (A) and in models in which time was not considered ("Average") or considered as a continuous ("Average change") or categorical ("Annual change") factor (independently of whether they were NB or ZINB mixed models) (B). a majority of provinces wildfires have been decreasing, in other provinces, including the ones comprehended in our study area, were increasing, mostly due to the contribution of small fires 5,56 . ZINB mixed models explained better the occurrence (presence/absence) of medium and large fires and the frequency of small fires. Moreover, ZINB mixed models substantially reduced the variance in the random components, after incorporating the zero-inflated part 44 , particularly for small and large fires. ZINB models have been previously used to model wildfires 33,[35][36][37] ; although other works that used ZINB and NB models found no evidence to favor the first over NB models 34 . Additionally, including time as categorical interacting factor improved model performance for all fire sizes. Notwithstanding, average models were also appropriate for large fires. Applying zero-inflated models to longitudinal data, including random effects, should be further considered to study wildfires in dynamic landscapes like the Mediterranean ones.
Of the seven sets of variables used for modelling wildfires, four of them were the ones producing best model results, which is supported by the literature: topography 57-60 , land use-land cover 20,25,61 and its changes 7,18,51,62 , and agrarian characteristics 18,33,53,63 . Other fire drivers as socio-economy, forest interfaces and linear features showed lower explanatory power for all fire sizes, which confirms other works 64 , even if these findings are not general 12,19,22,24,65 . Additionally, we could not ascertain an effect or trend of climate variables, despite changes underwent by some of them in southern Europe in recent decades 66 . On the other hand, there were large differences in model performance of these sets of variables, depending on fire size and modelling approach. While topographic variables had the best model performance for large fires (largest DIC reduction), other sets of variables produced similar results for small and medium fires (e.g., LULC changes types and their changes). Moreover, within a given set of variables, not any one variable produced the best results (highest DIC reduction) across fire sizes and either for NB or ZINB models. This permits arguing that this area is very complex from a perspective of wildfires 61,67 . Fires are ignited by people, with varying sensibilities and needs, and with changing capacities to fight fire over time. In addition, the role of the various variables changed over time. Initially, wildfires were concentrated in mountainous terrain, where farming was limited and the vegetation was dominated by conifers and other hazardous types, which underwent important changes (e.g., afforestation) over time. Wildfires were abundant in landscapes with smaller farms, and rare in areas with larger ones. By contrary, fires were limited where agricultural activities dominated, proper of the low elevation and flatter areas. With time, some variables decreased their statistical weight such as in conifer forest, or elevation, the latter being so mainly for medium and large fires. Other variables reversed their sign, as shrublands, hazardous stability and large farms. The net result was that wildfires moved to areas at lower elevation, with herbaceous crops and large size farms. Our findings support studies that found that topographic variables tended to have lower contribution to explain wildfires in the last decades 16 , as well as those studies that found that wildfires have decreased on forests 22,65 as well as on shrublands 25,68 , while increased on less flammable areas 22 .
To conclude, the initial link of wildfires to areas of certain hazardous characteristics, in which they abounded at the beginning of the study period, was lost. This process has also been observed in other Mediterranean areas of Greece 52 and Portugal 68 , in which recent fires have moved partially to non-fire-prone areas, indicating a departure from historical burning patterns. This again supports arguing that the complex nature of this human dominated landscapes, in interaction with the ignitions generated by people, complicates understanding of how fires will evolve over time and in space. In this regard, deterministic modelling approaches to infer future fires due to climate change or other global change drivers will be plagued with uncertainties. Consequently, if projections of future fires are difficult to make in general 69 , in Mediterranean landscapes, like the one studied, will be even more challenging.