Developing Models to Predict the Number of Fire Hotspots from an Accumulated Fuel Dryness Index by Vegetation Type and Region in Mexico

: Understanding the linkage between accumulated fuel dryness and temporal ﬁre occurrence risk is key for improving decision-making in forest ﬁre management, especially under growing conditions of vegetation stress associated with climate change. This study addresses the development of models to predict the number of 10-day observed Moderate-Resolution Imaging Spectroradiometer (MODIS) active ﬁre hotspots—expressed as a Fire Hotspot Density index (FHD)—from an Accumulated Fuel Dryness Index (AcFDI), for 17 main vegetation types and regions in Mexico, for the period 2011–2015. The AcFDI was calculated by applying vegetation-speciﬁc thresholds for ﬁre occurrence to a satellite-based fuel dryness index (FDI), which was developed after the structure of the Fire Potential Index (FPI). Linear and non-linear models were tested for the prediction of FHD from FDI and AcFDI. Non-linear quantile regression models gave the best results for predicting FHD using AcFDI, together with auto-regression from previously observed hotspot density values. The predictions of 10-day observed FHD values were reasonably good with R 2 values of 0.5 to 0.7 suggesting the potential to be used as an operational tool for predicting the expected number of ﬁre hotspots by vegetation type and region in Mexico. The presented modeling strategy could be replicated for any ﬁre danger index in any region, based on information from MODIS or other remote sensors.


Introduction
The quantification of the influence of varying fuel dryness on temporal fire occurrence risk is critical for improving decision-making in strategic fire management planning (e.g., [1][2][3]), especially under growing conditions of vegetation stress associated with climate change, which can alter the length and severity of the fire seasons (e.g., [4][5][6][7]).
One of the most used operational indices that integrates meteorological and satellite information is the Fire Potential Index (FPI) [28]. The FPI combines remotely sensed estimates of vegetation greenness-as measured by 10-day Normalized Difference Vegetation Index (NDVI) composites-with daily estimates of dead fuel moisture content [42,43] for mapping fuel dryness conditions and associated fire risk and danger. The FPI fuel dryness index has been operationally used for fire danger monitoring and occurrence risk prediction in the United States of America (USA) [28,[33][34][35], Indonesia [44], and on the European continent (e.g., [29][30][31]45]), including studies of regional application in northern Spain [46][47][48]. Several of these works have highlighted the need to understand how the same values of a fuel dryness index such as FPI result in different patterns of fire occurrence under different bioclimatic regions and vegetation types.
In Mexico, Manzo-Delgado et al. [49,50] demonstrated the potential of the temporal evolution of NDVI-based indices as indicators of fuel drought and associated fire risk for the central region of the country. Some studies have evaluated the role of weather variables such as precipitation or temperature on fire occurrence risk (e.g., [51][52][53]), mainly at regional or local scales. The work of Pompa-Garcia et al. [54] explored the association between fire records and the Standardized Precipitation-Evapotranspiration Index (SPEI) drought index, finding that the relationships of drought and fire vary by regions in the country. The FPI system [28] was tested by Sepulveda et al. [55,56] in the region of Baja California, but its performance has not been tested nationally for fire occurrence prediction.
Fire monitoring efforts in Mexico country include the development of an interface for active fire hotspots [57] near-real time monitoring [58,59]. The monitoring of fire hotspots is used for strategic decision making in fire suppression by considering the number of active fire hotspots as indicators of the level of fire ignitions and potential for forest fire spread in the main regions of the country.
Previous analysis of monthly fire MODIS hotspot and NDVI-based fuel greenness dynamics for the main vegetation types in Mexico [60,61] suggested that the relationship between fire hotspots occurrence and fuel dryness varied by month, with similar values of fuel dryness resulting in a higher number of fire hotspots as the fire season advanced. These results suggested that the accumulated dryness from the previous months might have potential for explaining the temporal evolution in fire hotspot occurrence.
Based on the above-mentioned results, this study provides further insight into the potential benefit of using an accumulated fuel dryness to explain forest fire occurrence variability. The current study tested an Accumulated Fuel Dryness Index (AcFDI), to predict 10-day (dekad) number of MODIS active fire hotspots, expressed as a Fire Hotspot Density index (FHD). The AcFDI was calculated from the temporal accumulation of a remote-sensing-assessed fuel dryness index (FDI)-based on the structure of the FPI index [28]-for the main ecoregions and vegetation types of Mexico.
Overall, the considerably wide variety of vegetation domains, meteorological conditions, and fire seasonality patterns existing in Mexico offered a valuable opportunity to explore the relationships of fire occurrence and fuel dryness under a variety of vegetation types, ecoregions, and weather conditions.
The specific goals of the current work are (1) to estimate thresholds for fire occurrence for the fuel dryness index (FDI) for each vegetation type and region in Mexico; and (2) to develop models to predict 10-day MODIS Fire Hotspot Density (FHD) for the main vegetation types and regions in Mexico.

Study Area
The area of study was the whole Mexican territory. Figure 1 shows the main vegetation types in the country according to the most recent land use map (INEGI Land Use Map Series V, year 2011, 1:25,000 http://www.inegi.org.mx/geo/contenidos/recnat/usosuelo/) of the National Institute of Geography and Statistics (INEGI in Spanish). Given the well-documented variations in fire regimes seasonality in the country (e.g., [62][63][64][65]), five geographical regions were established: Northwest (NW), Northeast (NE), North-Centre (NC), Centre (C), and South (S) (Figure 1). The definition of regions was based on the North American Level 3 Ecoregions Map (EPA, https://www.epa.gov/ecoresearch/ecoregions-north-america), together with previous analysis of the temporal trends and spatial patterns of clustering in fire hotspots in the country [54,60,61]. Vegetation types were reclassified into the following categories: temperate forest (FOR), dry tropical forest (DTROPF), seasonally dry tropical forest (SDTROPF), seasonally wet tropical forest (SWTROPF), wet tropical forest (WTROPF), wetlands (WET), desert shrubby vegetation (DSHR), natural pastureland (NPAS), and agricultural croplands (AGR) (Figure 1). Agricultural lands were included in our analysis because of the relevance of agricultural burns, including escaped burns from agriculture and slash-and burn agriculture in the occurrence of forest fires in the country-particularly in the Centre and South regions (e.g., [54,61,66]). The desert shrubby vegetation of all regions and natural pastureland of the North-Centre arid region were excluded from our analysis because of the relatively low registered number of hotspots and fire suppression registers in these vegetation types and regions. This resulted in a total of 17 vegetation and region units ( Figure 1).

MODIS Active Fire Hotspots
MODIS active fire hotspots for the study period were compiled from the Forest Fire Early Warning System operated by the National Commission for Knowledge and Use of Biodiversity (CONABIO) (http://incendios1.conabio.gob.mx/). Active fires were obtained based on the Contextual Fire Detection Algorithm for MODIS [57]. The algorithm classifies MODIS pixels as possible active fires by applying minimum thresholds for brightness temperature and considering the context of the surrounding pixels [57,67].
To minimize false detections, active fire hotspot data were filtered by CONABIO following the protocol described by Cruz-Lopez [67]. This process includes the consideration of specific thresholds for brightness temperature [67], in order to avoid false detection by sun-light and very high albedo (e.g., [68]). Several additional masks were also considered in the filtering, such as a NDVI and vegetation mask to remove desert areas without vegetation, or a stable light mask, generated from DMSP satellite images (http://julius.ngdc.noa.gov:8080/production/BIOMASS/night.html), as described in [67].

Fire Hotspot Density Index (FHD) Calculation
For each of the 17 vegetation types and regions considered, Fire Hotspot Density was calculated for each 10-day (dekad) period by dividing the number of MODIS hotspots observed by dekad in each vegetation/region by the surface (km 2 ) of the corresponding vegetation/region considered, scaled to a Fire Hotspot Density index (FHD) as follows (Equation (1)): The scaling factor was chosen as 1/10 of a reference surface of 200.000 km 2 , which is approximately equivalent to the surface of the temperate forests of the Northwest or Centre regions that are the main regions of forest fire occurrence in the country. Therefore, the FHD index directly represents 1/10 of the observed number of hotspots by dekad on a reference surface of 200.000 km 2 .

Fuel Dryness Index (FDI) Inputs
The inputs for calculating FDI were daily images of the estimated moisture content of 100-h time-lag dead fuels (H 100 ) [42] and 10-day NDVI composites. The daily H 100 images were calculated in CONABIO following the Cervera-Taboada [69] calibration of the NFDRS algorithms of dead fuel moisture content [42,43] for Mexico. This is based on temperature and relative humidity from MODIS atmospheric profiles combined with precipitation data from the TRMM satellite [58,59,69]. CONABIO calculated the 10-day MODIS NDVI composite images with a spatial resolution of 1 × 1 km, using the harmonic time series methodology (HANTS) proposed by Roerink et al. (2000), cited by De Badts et al. [70], as described by Cruz-Lopez [67]. The study timeframe was from 2011 to 2015, based on the availability of historical H 100 images calculated daily by CONABIO. We computed the average FDI value for 10-day periods for each vegetation type and region in the period of study.

Fuel Dryness Index (FDI) Calculation
The FDI was calculated based on the Fire Potential Index (FPI) developed by Burgan et al. [28], which integrates NDVI information with daily maps of dead fuels moisture content. FDI followed the general index structure proposed by Burgan et al. [28] (Equation (2)): where LR is the live ratio obtained from 10-day NDVI images following Equations (3)-(5), and MR is dead fuels moisture ratio obtained following Equation (6). Following Burgan et al. [28], the index represents the fuel dryness, with values close to 100 when both live and dead fuel moisture values are close to the maximum fuel dryness values, and lowest values when the fuels have higher fuel greenness and higher dead moisture conditions. Whereas some works have suggested a potential improvement by using the Normalized Difference Water Index (NDWI) instead of NDVI for the FPI calculation (e.g., [33]), we chose to utilize NDVI as in the original model formulation and several related works (e.g., [30,31,34,35]). The reason for this was the operational availability of 10-day MODIS NDVI compounds from the near-real-time fire monitoring systems of CONABIO [58]. This availability is advantageous both for the historical period and also for future real-time implementation of the fuel dryness index calculation.
The first component of the index-live ratio-was calculated based on relative greenness values following Burgan et al. [28] (Equations (3) and (4)): where RG is relative greenness, calculated as: where NDVI 0 is the observed NDVI for each pixel every 10-days, NDVI min and NDVI max are the absolute minimum and maximum NDVI values for each pixel from the period of study, and LR max is the maximum live ratio value for each pixel.
In the first FPI formulation [37], LR max was determined as a function of the live and dead loads, but this original proposal resulted in overestimated FPI values for some areas in the U.S. A second empirical formulation was proposed by Burgan et al. [28] for the calculation of LR max in the U.S., scaling this index from 35 to 75% based on the NDVI max map. For the FDI index calculation in this work, based on the observed number of hotspots for every LR value in Mexican vegetation types, we scaled LR max following Equation (5): where LR max is the live ratio for a given pixel when the vegetation is at maximum greenness, and NDVI max is the historical maximum NDVI observed for a given pixel. The values 125 and 255 are the absolute minimum and maximum NDVI values observed for Mexico, following Burgan et al. [28]. The second component of the FDI-MR-represents the ratio of dead fuel moisture. In the formulation of Burgan et al. [28] for the FPI, MR is calculated as a ratio of dead fuels moisture to fuels moisture of extinction based on the NFDRS fuels map for the U.S. Because no fuel extinction moisture map is available for Mexico, we tested the following expression for scaling MR for the FDI (Equation (6)): where H 100 is the observed value of 100 h dead fuel moisture, and H 100max , H 100min are maximum and minimum observed H 100 values for each pixel.

Threshold FDI p Values by Vegetation Type and Region
We calculated the accumulated percentage of fire hotspots by FDI for each vegetation type and region for the whole study period. Based on the accumulated percentage of fire hotspots, we calculated the minimum values of FDI above which 85%, 90%, and 95% of the accumulated FHD values occurred (FDI p , with p = 85%, 90%, 95%), similar to the work of Sebastián-López et al. [31].

Accumulated Fuel Dryness Index (AcFDI)
We tested a simple accumulated fuel dryness index AcFDI, based on the sum of the FDI values above the threshold for fire occurrence (FDI p ) in the previous three months. The index was assumed to be zero at values below the FDI p threshold (Equation (7)).
The index was calculated for each vegetation type and region as follows: where AcFDI i is the Accumulated Fuel Dryness Index for dekad i, FDI i is the observed Fuel Dryness Index for each dekad i, n is the current dekad, and FDI p is the Fuel Dryness Index threshold for fire occurrence, calculated for each vegetation type and region as the FDI value above which the p percentage of hotpots occurs, with p = 85%, 90% and 95% (Section 2.6).
The nine-dekad period for accumulating fuel dryness above the threshold was selected based on previous indices of accumulated fuel dryness derived for tropical countries (e.g., [71]), together with the observation of the FDI and FHD temporal evolution plots, by observing the average time difference between the start of fuel dryness accumulation above the FDI threshold and the observed starts and peaks of the fire season from FHD observations.

Models for the Prediction of Fire Hotspot Density (FHD)
We tested linear and non-linear models for the prediction of FHD i from FDI i or AcFDI i . Because the most robust results were obtained from non-linear models, here we report only results for these non-linear models. However, results for all linear models tested are available upon request from the corresponding author. The non-linear models tested were (Equations (8) and (9)): where FHD i is Fire Hotspot Density index for dekad i (Equation (1)), FDI i is Fuel Dryness Index for dekad i (Equation (2)), AcFDI i is Accumulated Fuel Dryness Index for dekad i (Equation (7)), and a and b are model coefficients to be estimated.
In addition, we also tested the following non-linear model, including an autoregressive term to account for the effect of the fires that occurred in the previous dekad FHD i−1 on the predicted FHD i values at dekad i, as follows (Equation (10)): The c coefficient for the autocorrelation of FHD was estimated as the calculated correlation coefficient between FHD i and FHD i−1 values for each vegetation type and region.
The models were fitted using non-linear least-squares regression (NLS) and non-linear quantile regression (NLQ) [72], using the packages nls and nlrq, respectively, of the software R [73]. For the quantile regression, we tested as candidate percentile values ranging 5 percentile units from 55 to 95%.
These candidate models were evaluated by using the coefficient of determination for nonlinear regression (R 2 ) as a statistical criterion (see [74], pp. 419 and 424), together with root mean square error (RMSE) and model Bias. R 2 is defined as the square correlation coefficient between the measured and estimated values (ry iŷi ) and the root mean square error (RMSE), which can be defined as follows (Equations (11)-(13)): where y i andŷ i are the observed and estimated values of the dependent variable, respectively, n is the total number of observations used to fit the model, and p is the number of model parameters.
Because the models will be used for fire risk strategic decision making, conservative predictions would be preferred. Therefore, in addition to highest R 2 and lowest RMSE values, we established as a criteria for selecting as best models those where the bias ranged from −5 to −15, therefore providing conservative overestimations for most of the observed values, as preferred for cautious risk decision making.

Observed Fuel Dryness Index for Vegetation Types and Regions
As expected, the observed temporal evolution of FDI varied by vegetation type and region. The observed FDI temporal trends for temperate forests of four regions (Figure 2 top plot) and for four tropical forests in the South region (Figure 2 bottom plot) for the five years of study are shown in Figure 2.  (Table 1).

Threshold FDI p Values by Vegetation Type and Region
Based on the accumulated percentage of hotspot distributions, we calculated the threshold FDI p values above which accumulated percentages values p of 85%, 90%, and 95% of the fire hotspots occurred for each vegetation type. The FDI p values for each vegetation type are shown in Table 1. Thresholds for FDI 95 are shown together with temporal evolution plots of FDI i in Figure 2 for some vegetation types. These plots allow visualization of the trends of fuel dryness accumulation above the threshold for the AcFDI index (Equation (7)). Once the FDI goes down to values below the threshold, fires are predicted to cease, according to the first condition of Equation (7). Therefore, these plots also allow visual interpretation of the predicted ends of the fire season each year based on FDI. FDI values in forests of the Northwest and Northeast generally started to decrease and reached values below the threshold one month later (June-July) than forests of Centre and South regions (May-June). The date of the FDI decrease also varied between years, possibly reflecting the differences in the occurrence of rain events between years and regions.

Predicting Fire Hotspot Density (FHD) from AcFDI by Vegetation Type and Region
Non-linear models based on fuel dryness showed higher goodness of fit statistics compared to linear models. The R 2 values for non-linear models using AcFDI and FDI (Equations (8) and (9), respectively) are shown in Table S1. For all vegetation types and regions, higher or similar predictive capacity for estimation of FHD was obtained using AcFDI compared to FDI (Table S1). The highest improvements for FHD using AcFDI instead of FDI were observed for temperate and tropical forests of Centre and temperate forests of the Northwest (NW), where R 2 values increased from 0.19, 0.21, and 0.35 to 0.51, 0.50, and 0.74, respectively (Table S1). Smaller gains or marginal increases were obtained for tropical forests of the South (S) and Northeast (NE) regions.
The best models for predicting fire FHD by vegetation type and region and the corresponding goodness of fit coefficients are summarized in Table 2. The best fitted models were obtained with AcFDI, calculated with FDI p threshold values of p = 95% (FDI 95 ) from Table 1, using a non-linear expression, including autoregressive terms (Equation (10)), for all vegetation types.
The chosen percentile for the best fit percentile regression models ranged from 65 to 90, with most vegetation types having a best fit for percentile values of 70 to 80. As defined in the criteria for model selection, the bias of the selected models in Table 1 ranged from −5 to −15, which was chosen to represent a conservative overestimation for cautions risk decision making. The R 2 values of the best fit models ranged from approximately 0.5 to 0.7 (Table 2), with the highest R 2 values for temperate forests of the Northwest and Centre regions and for dry tropical forest of the Centre region. The lowest values for several tropical forests of the South and Northeast regions, and for the temperate forest of the arid North-Centre region. RMSE values ranged from approximately 20 to 80, with most vegetation types having RMSE values in the range 30-60, as shown in Table 2. The highest values of the b coefficient were obtained for temperate forests of the South, North-Centre, Northwest, Centre, and Northeast regions; with b values ranging from approximately 2 to 4.7. This suggests that for these vegetation types, increases in accumulated FDI index resulted in strong increases in fire occurrence. In contrast, most of the tropical forests had lower b coefficients, ranging between 1 and 2, suggesting a less pronounced effect of AcFDI on FHD for these more humid regions.
The value of the c coefficient for autoregression of FHD was highest in the temperate and tropical forests of the Centre region, suggesting a larger effect of previously existing fires on fire occurrence in these regions.
Examples of the predicted FHD values against observed MODIS FHD with the selected models are shown in Figure 3 for several vegetation types.
The observed temporal evolution of FHD shown in Figure 3 varied by region. The Centre and South regions had earlier starts of the fire season, in the months of January and February, compared to a later start of the fire season in the Northwest region, where fire activity generally started in the months of March and April. (Figure 3).  Table 2 for predicted values).
Interestingly, for most vegetation types and years, predicted FHD showed a gradual increase from 0 (no risk) to levels of 50 (medium fire risk) in the first months of the year, prior to the start of the fire season (Figure 3). This might be useful as an early warning of the levels of fuel dryness accumulation before the start of the fire season, as will be briefly discussed below in Section 4.3.
For the temperate forests of the Northwest, years 2011, 2012, and 2013 had high-to-medium fire occurrence as shown in the observed FHD plots (Figure 3a), whereas years 2014 and 2015 had low and very low fire occurrence, respectively. This matches with the observed FDI values (Figure 2), where lower FDI values were reached in these latter years, and FDI remained above the FDI threshold for a shorter period compared to the first above-mentioned three years. In this region, the AcFDI index explained a substantial part of the variability in FHD, with R 2 values of more than 0.7 for the models with only FDI (Table S1).
For the temperate forests of the Northeast, the highest FDI values (up to 80) were achieved in the years 2011 and 2013, together with large periods of FDI well above the threshold (Figure 2), corresponding to years of high and medium FHD, respectively (Figure 3b). The fires occurring in this region in the very dry year 2011 were the largest fires recorded in the country, according to CONAFOR records, cited in [54]. Years 2012 and 2014-when lower FDI values of up to 70 were observed ( Figure 2)-corresponded to low fire activity (Figure 3b), and year 2015-where the period above the threshold was extremely reduced-had very low fire activity.
In the case of temperate forests from the Centre, both high FDI values and high fire occurrence were reached in the very dry year 2011. In contrast, 2014 was both a year with wetter fuel conditions as noted by lower FDI (Figure 2) and low fire activity as measured by FHD (Figure 3c). For the year 2015, in spite of wet fuel conditions as noted in the FDI plot, the observed fire FHD was higher than expected based on only fuel dryness. This seems to suggest a strong anthropogenic influence resulting in a large number of ignitions in spite of wet fuel conditions. In this region, the c coefficient for autoregression of FHD was higher than in the Northwest, and the increase in R 2 was higher with the addition of the autoregressive term compared to modeling with only AcFDI (R 2 of 0.7 for the model including autoregression of FHD (Equation (10), Table 2) compared to a value of 0.5 for the model with only AcFDI (Equation (9), Table S1).
For the wetlands of the South (Figure 3d), very high fire occurrence was recorded for the very dry year 2011 and high FHD was observed for the dry year 2013. Although the maximum FDI values reached ( Figure 2) did not vary strongly between years, the period where the FDI remained above the threshold for fire occurrence was longer in those dry years with higher fire occurrence, resulting in higher AcFDI values. Those values were larger than the AcFDI values in 2014, when fuel dryness accumulated above the threshold for a short period, resulting in a year of low fire occurrence. For the year 2015, the FDI index and the associated FHD predicted from the AcFDI seems to capture both the decrease in month 06 (June) and the rapid increase in the following weeks.

Observed Fuel Dryness Index (FDI) for Vegetation Type and Regions
The marked variation in observed FDI by vegetation type and region agrees with that reported from the works with the FPI index, which has a very similar structure to the FDI used here (e.g., [31]), and to previous works in Mexico that have reported distinct NDVI temporal cycles for different vegetation types (e.g., [50]).
The higher observed ranges of the FDI for the temperate forests of the Northwest, Northeast, and Centre regions (Figure 2), compared to the temperate and tropical forests of the South region ( Figure 2), suggest that less-marked drying and wetting processes occurred in the wetter regions and vegetation types, such as the temperate forests of the south or in the evergreen wet tropical forests. This agrees with the observations of Huesca et al. [46] in northern Spain, who found more marked seasonal fluctuations in the FPI temporal trends in the drier ecoregions, compared to smaller seasonal fluctuations in wetter ecoregions.
The temporal evolution of FDI from this work matches with previous observation of the monthly MODIS NDVI fuel dryness dynamics in the country [60,61], where a similar gradient in the dates of fuel dryness decrease from Centre and South regions to the Northwest region was observed. This is possibly related to the different timing of precipitation occurrence between these regions (generally one month later in the Northwest) (Figure 2), which corresponds well with the timing of the observed FHD decreases (Figure 3).
Compared to the temporal monthly NDVI based fuel dryness temporal evolution data shown in previous works (Vega-Nieva et al. [60,61]), the observed evolution of the weekly FDI in the current work showed more sharp increases and decreases ( Figure 2). This is not surprising: the NDVI temporal evolution are rather gradual, possibly reflecting accumulated dryness effects on the vegetation reflectance, which has been frequently associated with live vegetation dynamics in the literature (e.g., [33,41]). In contrast, the integration of the daily dead fuel moisture content ratio in the FDI index as proposed by Burgan et al. [28] implies a higher sensitivity to daily events of precipitation or raising temperature on the most rapidly changing dead vegetation moisture component of the index [42]. This sensitivity of the FDI to precipitation events could be advantageous compared to the monitoring of NDVI alone for establishing the predicted dates of end of the fire season based on the FDI thresholds, as illustrated in Figure 2 and Section 3.2.

Threshold FDI Values by Vegetation Type and Region
The thresholds obtained from the FDI distributions are similar to those observed in previous studies on the similar FPI index (e.g., [29,31,33]). The variation of FDI thresholds by region, with higher threshold values observed for drier regions and lower values for the wetter regions ( Table 1), supports that reported by Sebastián-López et al. [31] in a study for a number of bioclimatic regions in Southern Europe, where the highest threshold values for fire occurrence were also observed for the drier regions.

Predicting Fire Hotspot Density (FHD) by Vegetation Type and Region
The observed temporal evolution of FHD showed differences by vegetation type and region, coinciding with previous monthly FHD observations [61]. The observation of a more anthropogenic influence of fire occurrence-possibly from agricultural burns-in the Centre and South regions supports observations from previous works (Pompa-Garcia et al. [54] and Vega-Nieva et al. [61]). In these two regions, agricultural activities account for 40-50% of fires suppressed, as opposed to less than 20% attributed to this cause in the Northwest and Northeast regions, according to CONAFOR's fire suppression records for the study period.
The utilization of an accumulated fuel dryness index combined with percentile regression resulted in conservative estimates of the start of the fire season (Figures 2 and 3). Whereas all the NLS models underestimated a large number of the FHD observations, the use of percentile regression was more conservative. It has been acknowledged that using regression analysis can be problematic for these types of phenomena, because models equally penalize over-and under-estimations, which might have different implications for decision making in fire suppression [33]. The use of percentile regression [72] seems to offer a potentially useful alternative for more conservative predictions.
As noted in Section 3.3., the observed gradual increases in accumulated fuel dryness and associated FHD predictions (Figures 2 and 3) in the first months of the year-raising gradually from levels of 0 (no significant number of fires expected) to levels of 50 (medium fire risk) before the fire season starts-could be useful as indicators of the evolution of fuel dryness and potential fire risk early in the year. Because the effect of the autoregression on FHD before the start of the fire season is negligible, this gradual increase is entirely explained by the accumulated fuel dryness since the last rain event at the first months of the year. The monitoring of this early gradual increase in fuel dryness and associated fire occurrence risk might be useful for supporting operational decisions on fire suppression resources allocation with some advance prior to the start of the fire season. This predictive capability might be further enhanced in future works with the use of weather forecasts, similar to the works of Roads et al. [75,76], Chen et al. [77], or Preisler et al. [1,78], who developed monthly to seasonal estimates of fire danger from gridded weather forecasts.
Our study corroborated the nonlinearity of the relationship between the similar index FPI and number of fires observed by Sebastián-López et al. [30,31], and the coincidence of months and years with highest FHD values with FPI peaks, including bimodal fire occurrence, observed by Huesca et al. [46].
The strategy of developing specific models for different vegetation types and regions addressed in this study is consistent with the approach used in other studies with the FPI index, where different models were constructed for dissimilar ecoregions [31], sub-geographic areas of similar climate, fuels and topography [34,35], or combinations of vegetation types and climatic ecoregions [47,48], because of the largely documented influence of vegetation types on fire occurrence risk (e.g., [79][80][81][82]).
Our results seem to confirm the hypothesis of FDI saturation in the warmer regions and seasons suggested by Sebastián-Lopez et al. [31] for the similar FPI index. This also agrees with the observations from Vega-Nieva et al. [60,61], who found that similar monthly NDVI-based fuel greenness values resulted in a higher number of fire hotpots as fire season advanced in Mexico, suggesting an accumulated dryness effect on fire occurrence. A possible effect of the propagation of existing fires (modeled in this study by autoregressive coefficients) or a "contagion effect" in anthropogenic fires might also contribute to that result.
All of the best fit models in our study included autoregressive coefficients for FHD. Similarly, Huesca et al. [47], comparing models for predicting number of fires with only autoregressive terms, with other models combining both FPI and autocorrelation for number of fires, concluded that both autoregressive terms and FPI significantly explained the observed number of fires.
Although compared to other studies, our approach considered an ample ensemble of data, referred to a varied spectrum of ecological and weather conditions, in future works it will be necessary to validate this initial modeling effort with additional information from other years.
On the other hand, new sensors for active fire hotspot detection at higher spatial resolutions such as the 375 m VIIRS hotspots from NPP satellite [83], which have been used operationally in the country by CONABIO since 2015, should be tested in future works, once the period is long enough and covers representative variations in fuel dryness and fire occurrence.
Another limitation of the current work is inherent to the use of number of fire hotspots as indicators of fire risk. Whereas active fire hotspot products are used on a daily and weekly basis for decision support of fire suppression planning by the National Forestry Commission CONAFOR in the country, the total number of active fire hotspots is not the only variable considered for fire suppression decision making and resource prioritization. Similar to some works that have explored both active fire hotspots and fire suppression records from forest agencies to compensate limitations inherent to both types of datasets (e.g., [27,84]), future works might explore relationships between fuel dryness, fire hotspots, and fire suppression records from CONAFOR to develop operational recommendations for fire suppression resource allocations. Furthermore, the consideration of spatial variables such as road or population distance or density, which have proven of great value for developing spatial [85][86][87][88] and spatio-temporal fire occurrence risk maps (e.g., [82]), should be explored in future works.

Conclusions
The work presents the first effort towards developing an operational model for the prediction of the expected number of fire hotspots expressed as a fire hotspot density index, using remotely sensed weather and NDVI information, at a national level in Mexico. The models used percentile regression to predict 10-day fire hotspot density from accumulated fuel dryness and autoregression with previously observed fire hotspots. The reasonably good accuracies for most of the vegetation types and regions obtained suggest the potential to be included in operational GIS tools to assist in fire management decision making in the country.
It would be desirable that future works address the modeling of spatial patterns of fire occurrence from fuel dryness together with the consideration of human factors. The models could be tested in combination with forecasted fuel dryness-based on weekly or seasonal weather forecasts-to orient fire suppression strategic planning.
The modeling strategy presented here for predicting the expected number of fire hotspots by vegetation types and ecoregions from accumulated fuel dryness and satellite observed fire hotspots of the previous days could be replicated for any fire danger index in any region, based on records of daily active fire hotspots and available fuel dryness indices.