Introduction
Forest fires are a major disturbance in many parts of the world, with considerable impacts on vegetation and on human societies (Pyne, Andrews and Laven, 1996; Bowman et al., 2009). In sites with little direct human impact, variables associated with natural fire regimes, such as frequency, size, seasonality and intensity of fires, are governed by natural processes and biophysical conditions, such as climate, vegetation, and topography. In human dominated landscapes, anthropogenic activities modify fire regimes through actions such as fire suppression or an increase in the incidence of ignitions through activities associated with land use (Heyerdahl et al., 2003; Krebs, Pezzatti, Mazzoleni, Talbot and Conedera, 2010). Clearly, human pressure through land use and fire management activities have altered fire regimes both at local and global scales (Chuvieco, Giglio and Justice, 2008; Marlon et al., 2008; Martínez, Vega-Garcia and Chuvieco, 2009). Consequently, the current patterns of fire regimes have dramatically changed during the last decades (Moreno, Vázquez and Vélez, 1988; Westerling, Hidalgo, Cayan and Swetnam, 2006; Martínez-Torres, Cantú-Fernández, Ramírez-Ramírez y Pérez-Salicrup, 2015), becoming one of the main drivers of forest degradation.
Although human activities have clearly altered fire regimes, this process is not globally homogeneous and is dependent on the types of ecosystem and associated human activity. In fact, regional studies hint at a complex relationship between fire regime characteristics and human activities. For instance, Syphard et al. (2007), found in California that fire density was strongly associated with distance to the wildland-urban interface, and by population density and vegetation type, but fire frequency was only associated by vegetation type. Di Bella, Jobbágy, Paruelo and Pinnoc (2006) found that in South America, agriculture suppresses fire density in arid regions but enhances it in humid environments, suggesting that agriculture prevents biomass burning in semiarid areas but enhances it in humid environments, where biomass accumulates at faster rates. In general terms, Guyette, Muzika and Dey (2002), identified four ways in which humans influence the amount of burnt land (or the burned area fraction): anthropogenic ignitions, fuel production, fuel fragmentation, and cultural behavior. All these factors are linked to population density.
Human caused changes in fire regimes are particularly critical near or in natural protected areas. More important, understanding the effects of human activities on several aspects of a fire regime is essential and would allow for the generation of more realistic fire management plans. Protecting nature reserves from catastrophic disturbances associated with wildfires requires a management strategy that considers the driving factors of forest fire occurrence, and adjustment of local strategies. Under this approach there are several studies that have evaluated the relationships between different types of driving factors like landscape elements such as distance to roads with effects in the frequency of human-ignited fires (Alencar, Solorzano and Nepstad, 2004; Syphard et al., 2007). Also, fire occurrence appears to be positively associated with population density (Keeley and Fotheringham 2001; Guyette et al., 2002). Furthermore, the agricultural frontier provides abundant ignition sources, where fires frequently escape beyond their intended boundaries and ignite adjacent forests (Nepstad et al., 1999; Alencar et al., 2004).
There have been only few studies in Mexico that address the pattern of fire occurrence in coniferous forests based on landscape elements (Roman and Martinez, 2006; Ávila-Flores et al., 2010). Furthermore, Mexico is a country where a great proportion of forest resources is distributed in special land tenure system. In "ejidos" and "indigenous communities" land is owned by farmers and indigenous communities, decisions about forest management are taken collectively, and human presence in forested ecosystems is extremely conspicuous, particularly in the central and southern portions of the country (Martínez-Torres, Castillo, Ramírez and Pérez-Salicrup, 2016).
In this study, we evaluated which anthropic landscape factors are associated with the occurrence of fire in the Monarch Butterfly Biosphere Reserve (MBBR), where the implementation of fire management actions must overcome the daunting problems related to land use change (i.e., conversion to pastures and agriculture), illegal logging, and an inadequate management of forest resources (Ramírez, Azcárate and Luna, 2003; Merino and Hernández, 2004; Honey-Rosés 2009; Brower et al., 2016). We first assessed the driving factors from the landscape using historical fire records (2009-2013), as well as environmental and anthropogenic variables associated to fire ignition. Then, we developed a spatially explicit model of forest fire probability within the MBBR.
Objectives
The main objective of the current study is to assess the influence of several driving factors like anthropogenic and biophysical landscape elements on the performance of static, spatially explicit fire occurrence probability.
Materials and methods
Study area
Location. The MBBR is located in central Mexico, between Mexico and Michoacan States (Fig. 1). The MBBR constitutes the largest hibernation habitat of the monarch butterfly (Danaus plexippus L.), which migrates from the eastern portions of the United States and Canada to these forests (Slayback and Brower, 2007). Its territory encompasses a total area of 56 259 ha, divided in two buffer (42 707 ha) and three (13 551 ha) core zones (Comisión Nacional de Áreas Naturales Protegidas [Conanp], 2001). The core zones are the Altamirano Mountain in the north, the Chincua, Campanario, Chivati-Huacal, and Lomas de Aparicio Hills in the center, and the Pelón Mountain in the south. Cerro Altamirano constitutes a single mountain detached from the larger part of the reserve. In this study, we excluded it from our analyses.
Relief, Vegetation Types and Land Uses. The MBBR is part of the Trans-Mexican Volcanic Belt physiographic province (Ferrusquía-Villafranca 1990), which includes basically mountainous areas in central Mexico. It is characterized by a discontinuous mountain system, intensely dissected by strong tectonic processes, consisting of a set of hills and low hills, with a maximum elevation in the northern part of 3640 meters above sea level and with abrupt slopes. Vegetation types include dense conifer forests dominated or codominated by Abies religiosa (locally known as oyamel), Pinus spp., and Cupressus lusitanica; and mixed conifer and broad-leaved forests which include species of the genera Quercus, Alnus, Arbutus, and Cleyera, among others (Espejo, Brunhuber, Segura and Ibarra, 1992; Madrigal, 1994). Agriculture is the most important economic activity in the study area, where the dominant land uses are irrigated and rainfed agriculture with rotation cycles, followed by cattle pastures and forest management activities associated with timber extraction (Ramírez et al., 2003).
Forest Fire Regime. Knowledge regarding the forest fire regimes at MBBR is still limited. According to Martínez-Torres et al. (2016), the monitoring institutions of fire occurrence at MBBR are the National Forestry Commission (Conafor) and the Natural Protected Areas National Commission (Conanp) brigades' reports. Thereafter, Cantú (2013) assessed the reliability of the fire database of 2012. Most of the fires recorded were described as having low or moderate intensities in accordance with Agee (1993) and Sugihara, van Wagtendonk and Fites-Kaufman (2006).
Fire Database
We used the fire database from various federal institutions, such as Conafor and Conanp. The former government-based agency coordinates operations for fire prevention and control and keeps a registry of all fire events in the country along with latitude and longitude coordinates. We used information on fires from 2009 to 2013. The latter is the government agency in charge of natural protected areas. Only the fires of 2012 and 2013 were verified on the field. The heterogeneity of this fire data base makes it difficult to assess data quality in the extension of the surface area burned. Therefore, we used only records of geographic coordinates as occurrences fire events, with a total of 41 records of fire occurrence. We did not use the hotspot of MODIS sensor that are accessible through the Comisión Nacional para el Conocimiento y Uso de la Biodiversidad[Conabio] (2015).
Empirical Variables
Since we hypothesized that anthropogenic factors shape fire occurrence, we gathered data related to ignition sources and human presence (Table 1). To express these variables, we compiled spatial data of land use types (Ramírez-Ramírez, 2001), roads (Instituto Nacional de Estadística y Geografía [Inegi], 2011), human population per locality (Inegi, 2010), total tree biomass, and dead tree and stump density (Inventario Nacional Forestal y de Suelos 2009-2012; Conafor, 2012). We employed the digital elevation model from the Shuttle Radar Topography Mission (CGIAR Consortium for Spatial Information [CGIAR-CSI] undated) to estimate the slope and aspect factors. The friction map is based on the notion that movement in the terrain usually requires certain amount of effort and/or energy. Due to this "friction", spatial interactions tend to occur more often over short distances and low-grade slopes. We estimated this variable from road maps (scale 1:250 000) and 90 m. digital resolution elevation models (DEM) (Farfán, Mas and Osorio, 2012). We generated four classes for the friction map, 1- zero friction, 2 -moderate friction, 3- high friction, and 4 - very high friction. Finally, we defined edge density using a classification between forest and non-forest cover from the vegetation and land use map updated from the first version by Ramírez-Ramírez (2001) every three years. Therefore, we use the updated version for 2012 and then we estimated the mean forest cover binary raster (forest = 1) using ARCGIS (version 10.2) "focal statistics" tool with a 20 × 20 pixel window. The resulting raster contained values from zero to 1, where zero represents pixels in non-forest areas away from the edge of the forest, and one-valued pixels are away from forest edges but within the forest. As we wanted to assess the impact of edge, we set a value of 0.5 as the maximum value, where the influence of forest and non-forest is the same. Pixels with values > 0.5 were reclassified by their complement (1 - pixel value). The spatial resolution of the cells of each aforementioned variables was 90 m.
Types of variables | Variable | Description | Data Source | CODE |
---|---|---|---|---|
Environmental | Total biomass | Interpolated using IDW method by ArcGis 10.2 spatial analyst tool | Inventario Nacional Forestal y de Suelos (INFyS) | B |
Dead trees density | Number of dead trees per hectare interpolated using IDW method by ArcGis 10.2 spatial analyst tool | INFyS | DTD | |
Slope | Slope in degrees calculated by ArcGis 10.2 surface analysis tool | Shuttle Radar Topography Mission | S | |
Stump density | Number of stumps per hectare interpolated using IDW method by ArcGis 10.2 spatial analyst tool | INFyS | SD | |
Anthropic | Total population per locality | Interpolated using IDW method by ArcGis 10.2 spatial analyst tool | Instituto Nacional de Estadística y Geografía (Inegi 2010) | TP |
Friction | Four categories: 1- zero friction, 2 -moderate friction, 3- high friction, and 4 - very high friction | Farfán et al. (2012) | FR | |
Road density | Number of roads per hectare calculated using DINAMICA EGO | Inegi | RD | |
Distance to road | Euclidian distance in meters calculated by ArcGis 10.2 spatial analyst tool | Inegi (2011) | DR | |
Distance to core zone | Euclidian distance in meters | Comisión Nacional de Áreas Naturales Protegidas (Conanp) | DCZ | |
Distance to cattle pasture | Euclidian distance in meters | Ramírez-Ramírez (2012) | DP | |
Distance to agriculture | DA | |||
Distance to urban zone | DUZ | |||
Edge density | Average of the number of forest pixels in a 20 × 20 window calculated by ArcGis 10.2 with focal statistics tool | ED |
Data Analyses
Generalized Linear Mixed Models and Model Averaging. We analyzed the correlation degree between variables by Spearman coefficient analyses and assessed which variables were correlated. When two variables had a correlation coefficient ≥ 0.6 (r ≥ 0.6) we retained only one of the variables, and discarded the other to avoid variance inflation for correlated variables (i.e. redundant variables were excluded; Wintle, Elith and Potts 2005; Zhang, Lim and Sharples et al., 2017). We evaluated the remaining response variables by means of generalized linear mixed models (GLMMs). We fitted a binomial distribution with a logit function to build a set of models representing all possible combinations of one (N = 14), two (N = 77) and three (N = 285) response variables, and used the southern aspect as a random variable to control for variability among the explanatory variables most closely associated with fire occurrence (Demidenko, 1987). We corrected the models for overdispersion (Crawley, 2002), and dis-carded resulting models with non-significant intercept values and a variance inflation factor (VIF) ≥ 2 (Fox, 2008). We assessed the level of support of the remaining models (N = 44) through the small-sized bias-correction version of Akaike Index Criterion (AICc), and selected those with a ∆AICc < 2 (Burnham and Anderson, 2002). Then we estimated the conditional (GLMM R2c) and marginal (GLMM R2m) R2 of the best-supported GLMMs models, following Nakagawa and Schielzeth (2013) procedures. The first represents the amount of variability explained by the whole model (fixed and random effects), while R2m is the variability explained only by the fixed effects.
Finally, to identify the influential variables of fire occurrence, we built a subset of models representing all possible combinations of the response variables found in the best-supported models from above. We averaged all the models to obtain the mean and 95% confidence interval (CI) of the variables coefficients, and considered those variables whose values did not include zero as influential. We performed all statistical analyses in the program R (R Development Core Team, http://www.R-project.org/). We employed the coefficients derived from the best model to generate a fire probability map at one square kilometer grid cell resolution with a total of 836 cells.
Model Validation. The relative operating characteristic (ROC) analysis is a widely used quantitative method to assess the performance of spatial models that produce a probability map, which presents the sequence in which the model selects grid cells to determine the occurrence of a certain event (Pontius Jr. and Schneider, 2001; Mas, Soares Filho, Pontius, Farfán Gutiérrez and Rodrigues, 2013). In this study, ROC was employed to assess the effectiveness of the model in predicting the probability of occurrence of fire. ROC measures the accuracy of the model by comparing the map of predicted fire probabilities with the binary burned vs. unburned map of subsequent observed fires of 2014 and 2015. If the true events coincide perfectly with the higher ranked probabilities, then the area under the curve (AUC) is equal to one, since the ROC curve begins at the point (0,0), goes up the horizontal axis to the point (0,1), and shifts to the right reaching the point (1,1). A random probability map produces a diagonal ROC curve in which the true positive rate equals the false positive rate at all threshold points. Any probability map that has a ROC curve below the diagonal has less predictive power than a random map (Mas et al., 2013).
Results
Spearman Correlation Analysis
Two pairs of variables were correlated (r ≥ 0.6): road density (RD) with distance to roads (DR); and biomass (B) with density of dead trees (DTD) (Fig. 2). We selected distance to roads and biomass for inclusion in the GLMM regression analysis.
Generalized Linear Mixed Models
The GLMM regression analysis produced a prediction of fire occurrence with six best-fitted models (Table 2). The first best model suggested that this occurrence was significantly correlated with the variables ‘total population’ and ‘edge density’. Although models that are more complex explained greater variability, only population density and edge density were influential to explain fire occurrence, according to the averaged models (Fig. 3).
Modela | Intercept (SE) | AICcb | ΔAICc | - LLc | R2md | R2ce |
---|---|---|---|---|---|---|
TP + ED | ˗2.7 (0.96) | 103 | 0 | ˗ 48.32 | 0.28 | 0.28 |
S + TP + ED | ˗2.28 (1.03) | 103.8 | 0.84 | ˗ 47.63 | 0.30 | 0.30 |
DR + TP + ED | ˗3.31 (1.16) | 104.1 | 1.17 | ˗ 47.8 | 0.31 | 0.31 |
DUZ + TP + ED | ˗3.24 (1.17) | 104.6 | 1.61 | ˗ 48.02 | 0.28 | 0.3 |
DCZ + TP + ED | ˗2.79 (0.98) | 104.9 | 1.92 | ˗ 48.17 | 0.28 | 0.28 |
SD + TP + ED | ˗3.00 (1.09) | 104.9 | 1.95 | ˗ 48.19 | 0.28 | 0.28 |
aTP = Total population, ED = Edge density, S = Slope, DR = Distance to roads DUZ = Distance to urban zone, DCZ = Distance to core zone, SD = Stump density.
bAICc = Small-sized bias-correction version of Akaike Index Criterion
c-LL = Negative log-likelihood function.
dR2m = Marginal R2
eR2c = Conditional R2
(A and B) The points and the vertical lines represent the mean and 95 % confidence intervals (CI). *Influential variables, whose CI did not include zero. SD: stump density, DUZ: Distance to urban zone, S: slope, TP: total population per locality, DR: distance to roads, DCZ: distance to core zone, ED: edge density.
The best-fit model indicated that fire occurrence was favored where human population and forest fragmentation scored high values (Fig. 4). In Figure 5 it is possible to observe the spatial distribution of the values taken by these two variables.
The spatial model of fire occurrence from the best-fitted model (Table 2, Fig. 6) showed that the highest probabilities of fire occurrence are located at the surroundings of the MBBR core areas, where high total population and forest edge scores are present. We interpret the significance and direction of influence of these variables to indicate that increased human activity increased the probability of fire occurrence in the region.
The AUC of the ROC analysis of this model was 0.71 (Fig. 6), indicating a 71% of concordance between predicted probabilities and observed outcomes. The curve represents the cumulative tallies of false and true positives for all the cells in the modeled area. True positives are cells where the model and observed forest condition “burned” coincide.
Discussion
Identifying the driving factors of forest fires is essential for their management (Ávila-Flores et al., 2010). There are few studies that have focused on the causes of forest fire occurrence from the geographical perspective for Michoacan and Estado de Mexico states. Our results suggest that the occurrence of fire at the MBBR is influenced by a combination of factors that vary across the region and involve both physical and anthropogenic characteristics of the landscape. The significant positive correlation of fire occurrence with edge density is consistent with numerous studies where fires are more common along forest edges (Cochrane, 2001; Cochrane and Laurance, 2002). This association might be caused by the fact that forest edges are drier (Kapos, Ganade, Matsui and Victoria, 1993) and have higher tree mortality (Ferreira and Laurance, 1997; Laurance et al., 1997), than the core of a forest fragment, which might in turn increase fuel loads. In addition, forest edges in our study area are the boundary of agriculture and areas with cattle pastures, where land burning for crop and pasture regrowth is commonly employed as a management tool (Martínez-Torres et al., 2016). Fine litter, dead woody debris and canopy openness at forest edges might increase with the frequency and intensity of the adjoining land burning, thus creating a positive feedback where fire risk and danger increase. The spatial scale of fire occurrence is also driven by the shape of the forest fragments, whose complexity is correlated to its edge density. Hence, in the MBBR, total population, which in turn could be associated with an increase in the activities that generate forest fragmentation, may increase forest susceptibility to fire.
Total population per locality was an important variable in our model. This variable showed a positive influence on ignition occurrence, meaning a higher probability of ignition in the more populated areas, as was previously hypothesized. In regions where most fires are human-caused, this is a logical result, and in several other studies, total population, changes in population, and population density were found to positively relate to wildfire ignitions (Cardille, Ventura and Turner, 2001; Mercer and Prestemon, 2005). While these studies employed linear methods to investigate fire occurrence, studies using nonlinear methods often show that, after a threshold, fire occurrence tended to decrease with higher population densities. It appears that this is associated with declines in ignitions when human population density and development results in a reduction of forested area (Keeley, 2005; Syphard et al., 2007; Syphard, Radeloff, Hawbaker and Stewart, 2009). This pattern is not reflected by our results, probably because the sparse location of settlement within the MBBR, the persistence of forested areas, and the lack of localities with a high total population keep the relationship in the range of positive linear relation. Although the aforementioned studies have used population density, we used total population by localities and obtained a positive correlation. These results imply that fire managers must consider human influence together with biophysical characteristics such as those represented for the variable forest edge density when making decisions regarding the allocation of suppression and hazard mitigation resources. If human presence is not explicitly included in decision making, inefficient decisions may result, because fire occurrence is related to human presence on the landscape.
Our model explains roughly one-third of the total deviation in fire occurrence. Several reasons may be responsible for this limited explanatory power. First, it would be illusive to pretend that such a multifaceted phenomenon as the occurrence of forest fires in a highly complex topography can be explained by only a few variables. Other variables such as socio-economic aspects and forest management practices would be needed to develop a more complete understanding of the determinants of forest fire occurrence.
Second, because of the lack of direct data, proxies had to be used as predictors. We hypothesized that these proxies could influence fire drivers, which in turn might impact fire occurrence. From our analyses, we suggest that not all hypotheses regarding the proxies can be confirmed. For example, distance to pasture and distance to agriculture on fire occurrence.
Hosmer and Lemsow (2000) established that ROC is a threshold-independent measure of model discrimination, where a value of 0.5 indicates no discrimination. When 0.7 ≤ ROC < 0.8, discrimination is acceptable, and 0.8 ≤ ROC < 0.9 denotes an excellent discrimination. Alencar et al. (2004) developed empirical functions to relate occurrence of understory forest fire to landscape features, throughout regression analyses, separating for years of the El Niño Southern Oscillation (ENSO), and non-ENSO years. The ROC of the model that best predicted understory fires in both ENSO and non-ENSO years, was 0.778 and 0.751 respectively. Alonso-Betanzos and others (2003) developed a fire prediction model based on neural networks. They assessed AUC ROC curves for several years, ranging from 0.757 to 0.915. Lozano, Suárez-Seoane and de Luis (2007) assessed several spectral indices derived from satellite image data to modeling fire occurrence probability throughout logistic regression. The ROC values varied from 0.70 to 0.85 for all the assessed regression models. As mentioned earlier, our result of 0.71 indicates that the model has acceptable discrimination.
To discern the relative weight of each variable more accurately, we suggest integrating other analytical methods. Ouyang, Yung Han and Tong (2016) have found that machine learning algorithms, like random forest or conditional inference trees, showed greater predicting power of AUC than logistic regression. These also prevented bias detected in classification and regression tree algorithms. However, to achieve this, they used over 1,000 sample points. On the other hand, the three-dimensional kernel density estimator, a nonparametric descriptor tool, has been used to draw up smoothed maps which represent the continuous spatial density distribution of forest fires, including its time evolution (Tonini, Pereira and Parente, 2017).
Finally, we must underscore the lack data on fuel loads. Fuel loading is one of the three basic elements (oxygen, heat and fuel) that fires need to ignite. Fuel quantification can be done through direct measures in periodic forest inventories. The inclusion of this factor may be essential in modeling fire occurrence in the study area. It is clear that models will not improve considerably unless fuels loads are quantified, particularly in response to fire management practices.
Conclusions
This study highlights the importance of examining multiple causes of forest fire occurrence in the emblematic MBBR. Fires at the MBBR are not randomly spread throughout the landscape but are mediated by the condition of the remaining forest and present human activity. Specifically our study shows that human influences, such as land use change and the creation of forest edges, is strongly related with fire occurrence. The mapped probabilities of fire occurrence show the places where such condition is present. Forest fragmentation combined with slope and total population, correlates with a high likelihood of forest fires. These factors should be considered as fundamental elements in the design of forest fire prevention programs.
To increase the certainty of the model it would be necessary to incorporate other explanatory factors which were not analyzed here (e.g., fuel loadings, soil humidity, etc.) in combination with a longer period of fire occurrence records. This would produce a more comprehensive understanding of the factors underlying fire occurrence patterns and their spatial relationships. Nevertheless, the probabilities estimated here can be useful in the development of a fire management plan at MBBR.