Introducing Spatially Distributed Fire Danger from Earth Observations (FDEO) Using Satellite-Based Data in the Contiguous United States

: Wildﬁre danger assessment is essential for operational allocation of ﬁre management resources; with longer lead prediction, the more e ﬃ ciently can resources be allocated regionally. Traditional studies focus on meteorological forecasts and ﬁre danger index models (e.g., National Fire Danger Rating System—NFDRS) for predicting ﬁre danger. Meteorological forecasts, however, lose accuracy beyond ~10 days; as such, there is no quantiﬁable method for predicting ﬁre danger beyond 10 days. While some recent studies have statistically related hydrologic parameters and past wildﬁre area burned or occurrence to ﬁre, no study has used these parameters to develop a monthly spatially distributed predictive model in the contiguous United States. Thus, the objective of this study is to introduce Fire Danger from Earth Observations (FDEO), which uses satellite data over the contiguous United States (CONUS) to enable two-month lead time prediction of wildﬁre danger, a su ﬃ cient lead time for planning purposes and relocating resources. In this study, we use satellite observations of land cover type, vapor pressure deﬁcit, surface soil moisture, and the enhanced vegetation index, together with the United States Forest Service (USFS) veriﬁed and validated ﬁre database (FPA) to develop spatially gridded probabilistic predictions of ﬁre danger, deﬁned as expected area burned as a deviation from “normal”. The results show that the model predicts spatial patterns of ﬁre danger with 52% overall accuracy over the 2004–2013 record, and up to 75% overall accuracy during the ﬁre season. Overall accuracy is deﬁned as number of pixels with correctly predicted ﬁre probability classes divided by the total number of the studied pixels. This overall accuracy is the ﬁrst quantiﬁed result of two-month lead prediction of ﬁre danger and demonstrates the potential utility of using diverse observational data sets for use in operational ﬁre management resource allocation in the CONUS.


Introduction
Wildfires can result in tremendous economic loss across the United States. For example, the year 2018 was the largest fire year on record, resulting in approximately $3 billion in suppression costs (National Interagency Fire Center (NIFC) [1]). Fire suppression cost depends on multiple factors, including vegetation type and weather conditions and a successful fire suppression will reduce the area burned [2]. Contributing to these high suppression costs are logistical costs of relocating resources (personal comm. with Ed Delgado of NIFC). As such, there is a need to predict where fires will occur in order to optimize the timing and distance of relocation. The current fire danger predictions that inform resource allocation are primarily based on meteorological forecasts and expert judgement ( [3]). There are daily and weekly predictions produced using the National Fire Danger Rating System (NFDRS; [4], as well as 1-, 2-, 3-and 4-month Wildland Fire Potential Outlooks (FPOs). The US NFDRS [4] generates 1-day ahead fire danger predictions. The predictions use regional reanalysis meteorological input [5,6]. FPOs report fire potential in qualitative categories, such as below-normal, normal and above-normal. These qualitative metrics are derived using expert knowledge informed by meteorological forecasts and fuel maps. In both predictions, there is a reliance on meteorological forecasts, which have several limitations. For example, meteorological inputs rely largely on ground-based observations (e.g., weather stations), which have uneven spatial distribution and sparse sampling in remote regions [7]. In addition, the meteorological projections are often downscaled in wildfire analysis, especially in diverse regions with complex fire behavior [8]. Furthermore, the meteorological forecasts may not provide reliable estimates beyond 10 days [9]. Real-time remotely-sensed data, on the other hand, do not have these limitations and can potentially be used to improve the accuracy of wildfire danger predications and help fire managers make effective resource allocation decisions prior to wildfire occurrence.
A number of past studies have developed models to predict wildfire danger using hydrologic and ecologic variables. For example, one study has looked at seasonal predictions of wildfire burned area as below-normal, normal, and above-normal in the western United States using the Palmer Drought Severity Index (PDSI; [10]. In another study, the Southern Oscillation Index (SOI) of El Niño-Southern Oscillation (ENSO) has been used to predict the historical wildfire season severity in Everglades National Park in southern Florida with three-month and one-year lead [11]. NOAA AVHRR-based Normalized Difference Vegetation Index (NDVI) data were used to develop a model of dynamic fire danger in the Mediterranean regions [12]. NDVI has also been used to predict fire danger on Tenerife Island during the period 1992-1995 [13]. PDSI and logistic regression technique were used to predict the number and location of large fire events (> 400 ha) one month ahead of time [14]. Building on these studies, PDSI has also been used to forecast fire activity in Canadian and Alaskan forests during the fire season. In 2011, PDSI and Sea Surface Temperature (SST) were used to predict seasonal forest fire activity in Canada [15]. In 2014, an algorithm for predicting above-normal summer wildfire activity in southern Europe was developed based on meteorological drought index Standardized Precipitation Index (SPI) [16]. Furthermore, Actual Evapotranspiration (AET) as a proxy for fuel amount and water deficit (WD) as a proxy for fuel moisture have been used for predicting fire-burned area [17]. Vapor Pressure Deficit (VPD) has shown robust relationships with the variability of burned area in western United States [18] and has been found useful for early detection of drought onset [19]. In 2017, model-assimilated pre-season (January-April) soil moisture from NASA Gravity Recovery and Climate Experiment (GRACE) observations of Terrestrial Water Storage showed that soil moisture in months preceding the fire season is strongly related to both the number of wildfires and wildfire burned area [20]. Building on studies that linked VPD and soil moisture to wildfire burned area, a recent study used satellite-derived VPD and soil moisture to predict burned area in the nine US Geographic Area Coordination Centers (GACCs) [21].
Despite previous studies, gaps have remained in producing algorithms for operational use for allocating fire management resources. The previously developed models generally do not cover the entire CONUS, are not delivered in temporal scales needed for national resource allocation, and often rely on indirect hydrological observations. Unlike previous studies that do not deliver wildfire prediction models for entire CONUS [11,18,22], this study introduces a spatially distributed wildfire prediction model for CONUS. This is especially important as CONUS is the primary scale of national resource allocation. Furthermore, most previously developed models are not at temporal scales of operational resource allocation. The previous works focus on seasonal [10] or pre-season wildfire danger predicting [20], which are not temporal scales of operational resource allocation.
In this study, we developed a two-month lead wildfire danger model. Although our model exceeds the preferred one-month target lead time needed for operational resource allocation, it is still Remote Sens. 2020, 12, 1252 3 of 19 at a shorter lead time, of compared to the previous models. The two-month lead time was selected to account for potential satellite data latency. Finally, most of the previous wildfire prediction models rely on indirect observations of hydrology (e.g., PDSI drought indicator uses precipitation and temperature from reanalysis). This study, however, uses direct satellite hydrological observations for predicting fire danger. Remote sensing observations have generally less data latency compared to reanalysis and are important for operational deployments of resources. In addition, satellites produce crucial observations where access to ground-based hydrologic observations is scarce. In this study, we used remotely-sensed VPD, surface soil moisture (SSM) and enhanced vegetation index (EVI) for predicting wildfire danger in the CONUS using a spatially distributed framework and two-months' prior observations from satellites. The objective is to fill the existing gaps by developing Fire Danger from Earth Observations (FDEO), an algorithm for two-month lead operational prediction of fire danger defined as expected area burned as a departure from "normal" (defined as the mean over historical record). This definition aligns with how the decision support system defined fire danger (i.e., preparedness levels).

Datasets
FDEO uses three types of remote sensing data from Moderate Resolution Imaging Spectroradiometer (MODIS), Atmospheric Infrared Sounder (AIRS), GRACE, and two other datasets of historical fire data and land cover (Figure 1) to construct a statistical model. The first dataset is EVI from MODIS satellite (MOD13C2 [23]). The MODIS EVI product has been reported over land globally since 2000 in 1-km and 500-m pixel resolutions every 16 days [24]. EVI is a modified version of NDVI. NDVI was initially introduced by [25] and is based on the differences in pigment absorption features in the red (~0.660 µm) and near-infrared (NIR~0.860 µm) regions of the electromagnetic spectrum. EVI is more responsive than NDVI to canopy structural variations, including leaf area index, canopy type, plant physiognomy, and canopy architecture [24,26]. Studies show that EVI improves the sensitivity in high biomass regions, while NDVI tends to saturate in those regions [24]. We mapped EVI data using linear interpolation onto 0.25-degree spatial resolution. We also generated an accumulated EVI metric since the start of the record. For this purpose, we calculated the sum of EVI since the beginning of the record.
The second remote sensing dataset used in this study is VPD, derived from the AIRS satellite instrument. In this study, we use monthly Level 3 Version 6 VPD data at a spatial resolution of 1 degree (2002-present) [27,28]. We resampled this data using linear interpolation to 0.25-degree spatial resolution.
The third remote sensing dataset is SSM, derived from the physically-based Catchment Land Surface Model (CLSM) that assimilates ground and space-based meteorological observations [29][30][31]. One of the primary datasets used for generating SSM is terrestrial water storage anomalies (TWSA) from the Gravity Recovery and Climate Experiment (GRACE; [32]. Although GRACE data have relatively low spatial resolution (>150,000 km 2 ), its assimilation into CLSM produces SSM estimates at 0.25-degree spatial resolution and in monthly temporal resolution. The data record has been available since April 2004.
For model training and validation, we used the National Fire Program Analysis Fire-Occurrence Database (FPA FOD), which provides verified historical fire burned area information in the United States [33]. This database excludes prescribed and agriculture fires, reporting only wildfire information from federal, state, and local agencies. There are several attributes associated with this dataset, including location, discovery date, final fire size, and fire cause [33,34]. The data available are from 1992 to 2013. These data were aggregated into 0.25 × 0.25 degree cell spatial resolution at monthly time intervals. The operational needs for a fire danger metric require predicted burned area as a departure from normal. As such, we developed a climatology of burned area to provide a benchmark for the  The second remote sensing dataset used in this study is VPD, derived from the AIRS satellite instrument. In this study, we use monthly Level 3 Version 6 VPD data at a spatial resolution of 1 degree (2002-present) [27,28]. We resampled this data using linear interpolation to 0.25-degree spatial resolution.
The third remote sensing dataset is SSM, derived from the physically-based Catchment Land Surface Model (CLSM) that assimilates ground and space-based meteorological observations [29][30][31]. One of the primary datasets used for generating SSM is terrestrial water storage anomalies (TWSA) from the Gravity Recovery and Climate Experiment (GRACE; [32]. Although GRACE data have relatively low spatial resolution (> 150000 km 2 ), its assimilation into CLSM produces SSM estimates at 0.25-degree spatial resolution and in monthly temporal resolution. The data record has been available since April 2004.
For model training and validation, we used the National Fire Program Analysis Fire-Occurrence Database (FPA FOD), which provides verified historical fire burned area information in the United States [33]. This database excludes prescribed and agriculture fires, reporting only wildfire information from federal, state, and local agencies. There are several attributes associated with this dataset, including location, discovery date, final fire size, and fire cause [33,34]. The data available are from 1992 to 2013. These data were aggregated into 0.25 × 0.25 degree cell spatial resolution at monthly time intervals. The operational needs for a fire danger metric require predicted burned area as a departure from normal. As such, we developed a climatology of burned area to provide a benchmark for the spatial distribution of monthly burned area. We developed the climatology from the FPA FOD for the period of 2004-2013 ( Figure 2). Predictive models were developed for each land cover type. Land cover type data were obtained from the 2011 USGS National Land Cover Database (NLCD 2011) [35]. The original database has a spatial resolution of 30 meters. A reclassification was made to nine categories: Evergreen, Shrubland, Deciduous, Herbaceous, Wetland, Developed/Urban, Barren Land, Planted/Cultivated, and Mixed Forest [20]. We resampled the land cover map to 0.25 × 0.25 degrees using the dominant land cover type in each grid and used only vegetated land cover types-Evergreen, Shrubland, Deciduous, Herbaceous, and Wetland-which are critical to wildfire occurrence. Finally, we used the time period of 2004-2013 to develop the predictive models.

Converting the Hydrologic and Vegetation Variables to Standardized Indices (SI)
Depending on the ecosystem, climate-fire relationships are dependent on either drought or anomalously high moisture conditions in the months prior to the fire occurrence [36]. Standardization provides a baseline for monitoring both stressed and non-stressed hydrologic conditions. In this study, we therefore created standardized indices (SI) of SSM, VPD, and EVI. The standardized indices are derived based on the non-parametric approach introduced by [37]. The methodology uses the Gringorten plotting position as: where Var is the remote sensing variable (VPD, SSM, or EVI), i is the rank of data sorted from either the smallest (SSM and EVI) or the largest (VPD), and n is the sample size. Afterwards, the standardized indicators are derived as: where Φ −1 is the inverse normal distribution function with the mean of zero and standard deviation of one and SI is the standardized drought index. A negative SI is an indication of below average SSM, EVI, or above average VPD and thus stressed hydrologic (ecologic) conditions, while a positive SI, is an indication of above average SSM, EVI, or below-average VPD and thus non-stressed hydrologic (ecologic) conditions. One advantage of such standardized indicators is that they can be derived in different time scales e.g. one month, three months, or six months, and compared spatially.

Model Development
FDEO generates two-month lead expected area burned as a departure from the "normal". While the output of the model is similar to that developed in previous literature [10,36,38], we developed a wildfire prediction model for each land cover type using regression between the normalized burned area (y) and the best-fit (highest coefficient of determination) single explanatory variable (x). This method demonstrates the use of previously established land cover-dependent relationships between hydrological variables and area burned based on retrospective analyses into an algorithm for operational prediction. In our algorithm, we built models based on operation prediction requirements, such as: (1) predictor variables aligned with decision-making processes, (2) explanatory variables based on remote sensing hydrological variables that are openly and consistently observed globally, and (3) the latency between time of observation and availability for use in a decision-support tool. There were 16 explanatory variables tested; four for each of the remote sensing variables: EVI, accumulated EVI, VPD, and SSM from two months (to account for data latency between observation and availability) prior to burned area: (1) the standardized value by removing the mean and the standard deviation, (2) 1-month SI, (3) 3-month SI, and (4) 6-month SI. As shown in Figure 3, the histogram of cumulative burned area per best SI is plotted in the left column. The horizonal axis is the standardized index (SI) and each row represents one land cover type. The normalized burned area (right column) is the binned cumulated burned area (left column) divided by the number of samples (middle column) in each SI bin. We then derived the best fitted regression lines for each normalized burned area histograms (red lines). The methodology used in this study relies on the regressed relationship between normalized burned area and the best single explanatory variable out of the 16 variables tested. Rows represent different land cover types. The left column is the histogram of cumulative burned area, the middle column is the number of samples in each bin. The right column is the normalized burned area histogram. The horizonal axis is the standardized index (SI). We show deciduous, shrubland, herbaceous, evergreen, and wetland land cover, as each relied on a different explanatory variable. The red lines are the best fitted lines for the normalized burned area histograms.
To categorize observations and predictions as below normal, normal, and above normal, we first subtracted the average historical burned area ( Figure 2) from the burned area observations and predictions data in each grid cell, respectively, to derive burned area anomalies ( Figure 4). Next, we used the empirical cumulative distributions, separately for both the observed and predicted burned area anomalies ( Figure 5), to bias adjust the prediction distributions and determine the respective thresholds for categorizing each class. The Gringorten plotting position function (Equation 2) was used to calculate the empirical cumulative distributions. The 33 rd and 66 th percentiles were used as thresholds classifying below normal, normal, and above normal, three statistically-similar ranges. The methodology used in this study relies on the regressed relationship between normalized burned area and the best single explanatory variable out of the 16 variables tested. Rows represent different land cover types. The left column is the histogram of cumulative burned area, the middle column is the number of samples in each bin. The right column is the normalized burned area histogram. The horizonal axis is the standardized index (SI). We show deciduous, shrubland, herbaceous, evergreen, and wetland land cover, as each relied on a different explanatory variable. The red lines are the best fitted lines for the normalized burned area histograms.
To categorize observations and predictions as below normal, normal, and above normal, we first subtracted the average historical burned area ( Figure 2) from the burned area observations and predictions data in each grid cell, respectively, to derive burned area anomalies ( Figure 4). Next, we used the empirical cumulative distributions, separately for both the observed and predicted burned area anomalies ( Figure 5), to bias adjust the prediction distributions and determine the respective thresholds for categorizing each class. The Gringorten plotting position function (Equation (2)) was used to calculate the empirical cumulative distributions. The 33rd and 66th percentiles were used as thresholds classifying below normal, normal, and above normal, three statistically-similar ranges.  To evaluate the performance of the model against the observation, we used two statistical metrics: Root Mean Squared Error (RMSE) and overall accuracy. RMSE was used to evaluate the probabilities of prediction versus the observation. RMSE is defined as:   To evaluate the performance of the model against the observation, we used two statistical metrics: Root Mean Squared Error (RMSE) and overall accuracy. RMSE was used to evaluate the probabilities of prediction versus the observation. RMSE is defined as: To evaluate the performance of the model against the observation, we used two statistical metrics: Root Mean Squared Error (RMSE) and overall accuracy. RMSE was used to evaluate the probabilities of prediction versus the observation. RMSE is defined as: Remote Sens. 2020, 12, 1252 where LC is the land cover type, n is the sample size, p i is the prediction probabilities, and o i is the observation probabilities. p i and o i are empirical cumulative probabilities of prediction and observation derived from Figure 5. We calculate RMSE at each time step and for each land cover type.
To evaluate the categorical wildfire prediction models against the observation, we apply the overall accuracy metric [39]. This metric is used to assess the performance of discrete models. For the model associated with each land cover type, we calculate overall accuracy defined as: where LC is the land cover type. The higher the overall accuracy, the better the model performance, compared to the observation. Similar to RMSE, overall accuracy is calculated at each time step and for each land cover type. We will also show the performance of the categorical models in each of the below normal, normal and above normal categories using the confusion matrices.

Results
In Figure 2, the burned area climatology for the United States shows the spatial progression of the burned area from the East to the West through the entire fire season (January through December). Although more burned area occurs in the West later in the fire season, there is consistent wildfire in the eastern and southern US year-round. In the southwest, fire activity starts around March and peaks in June. In the Northwest, the majority of fires happen in July and August. In the Southern and Northern California, there is wildfire activity from May through October, peaking in October and July or August, respectively. In the Northern Rockies, fire activity peaks in the summer months: July and August. Not only does the peak month for burned area vary regionally, but also peak burned area. For example, larger burned area occurs in the West, while smaller burned area occurs in the East.
The best model for each land cover type was a quadratic regression and the explanatory variable differed by land cover (Figure 3). The best SI for Deciduous Forest, Shrubland, and Herbaceous land cover types is 1-month SI SSM (R 2 = 0.72), 1-month SI EVI (R 2 = 0.79), and 3-month SI VPD (R 2 = 0.78), respectively. The best SI for Evergreen and Wetland land cover types are 3-month SI VPD (R 2 = 0.68) and 3-month SI SSM (R 2 = 0.66), respectively. As observed in Figure 3, in Deciduous Forest, Herbaceous, Evergreen Forest, and Wetland land cover types, the relationship between the Standardized Drought indicators and wildfire burned area is negative. In other words, the likelihood of wildfire burned area increases with low SI values (i.e., hence stressed vegetation conditions). In Shrubland, however, the relationship between SI EVI and burned area is positive. High SI EVI values (i.e., more green in the months prior) are associated with higher wildfire danger. Figure 4 shows the observation and prediction anomalies of August and January 2013. When predictions were converted to burned area anomalies (Figure 4), there is spatial alignment between the observed and predicted below-normal burned area anomalies; however, there exist inconsistencies between observed and predicted anomalies representing above-normal and normal. Explicitly, several regions show around normal fire observation activity, while the prediction map shows above normal burned area. This is specifically the case for January 2013. As seen, there are vast areas in the western parts of the US where the observed anomalies are around normal, while the predicted anomalies are above normal. In order to adjust for the bias in over-predicting burned area as above normal, we used the observed and predicted empirical cumulative distributions ( Figure 5) to find the relative thresholds of the lower 33% (below-normal), 33-67% (normal), and greater than 67% (above normal). The bias-adjusted thresholds were necessary to align the similar regions of the distribution. Figure 5 only shows the bias-adjusted thresholds for August 2013. For example, in the wetlands, there is a very steep slope in the normal region of the distribution; the bias adjusted thresholds enable classification of the same region of the distribution, despite a slightly lower slope in the predicted burned area anomalies.
For the final results, we produce maps of the best ( Figure 6) and worst (Figure 7) burned area predictions, both as a probability and as categorical burned area anomaly. August 2013 was one of the best predicted months with an overall accuracy of 75%. Figure 6 shows consistent patterns between the observed and predicted burned area anomaly probability and categorical maps. The middle regions of the US show high probabilities of burned area for both the observed and predicted, while the southeast and northwest show consistent patterns of low probability of burned area anomaly. In the northeast, prediction and observation probability as well as categorical maps show above-normal wildfire activity. As shown, the categorical maps match the spatial patterns of the probability maps. Although overall accuracy (52%) was lower in January 2013, Figure 7 shows that the probabilities and categorical maps of January 2013 show more consistent patterns after the bias-adjustment. While the bias-adjusted observation and prediction maps (Figure 7) show significant improvement compared to the unbiased maps (Figure 4), there are still minor inconsistencies in the wildfire pattern as well as wildfire danger probability values between the observed and predicted maps. For example, the observed probability map shows very high fire danger in most parts of shrubland regions of Nevada, Idaho, Washington, Colorado, and Wyoming, where the prediction map indicates less fire danger.
For the final results, we produce maps of the best ( Figure 6) and worst (Figure 7) burned area predictions, both as a probability and as categorical burned area anomaly. August 2013 was one of the best predicted months with an overall accuracy of 75%. Figure 6 shows consistent patterns between the observed and predicted burned area anomaly probability and categorical maps. The middle regions of the US show high probabilities of burned area for both the observed and predicted, while the southeast and northwest show consistent patterns of low probability of burned area anomaly. In the northeast, prediction and observation probability as well as categorical maps show above-normal wildfire activity. As shown, the categorical maps match the spatial patterns of the probability maps. Although overall accuracy (52%) was lower in January 2013, Figure 7 shows that the probabilities and categorical maps of January 2013 show more consistent patterns after the biasadjustment. While the bias-adjusted observation and prediction maps (Figure 7) show significant improvement compared to the unbiased maps (Figure 4), there are still minor inconsistencies in the wildfire pattern as well as wildfire danger probability values between the observed and predicted maps. For example, the observed probability map shows very high fire danger in most parts of shrubland regions of Nevada, Idaho, Washington, Colorado, and Wyoming, where the prediction map indicates less fire danger. The spatial patterns of RMSE and Overall Accuracy (Figure 8) are consistent with those in the probabilities and categorized maps. Generally, the RMSE was lower and overall accuracy was higher in August 2013, compared to January 2013. The map shows that the RMSE for August 2013 has been relatively low for the majority of the regions in the US. The average RMSE for the entire CONUS is 0.2, with the lowest of 0.15 (for Herbaceous), and a highest RMSE of 0.3 (for Wetland). Once converted and categorized as below normal, normal, and above normal, the prediction has relatively high overall accuracy, showing a 75% average for the CONUS, 49% minimum (for Wetland), and 82% maximum (for Herbaceous) overall accuracy in August 2013. There is some spatial variability to the overall accuracy with slightly higher values in the Western half of the US, as compared to the Eastern half. For January 2013, the average RMSE for the entire US is 0.35, with the minimum and maximum of 0.30 (for Wetland) and 0.42 (for Herbaceous), respectively. The overall accuracy of January 2013 is 45% for the entire US, with 37% minimum (for Herbaceous) and 55% maximum (for Wetland) overall accuracy. The spatial patterns of RMSE and Overall Accuracy (Figure 8) are consistent with those in the probabilities and categorized maps. Generally, the RMSE was lower and overall accuracy was higher in August 2013, compared to January 2013. The map shows that the RMSE for August 2013 has been relatively low for the majority of the regions in the US. The average RMSE for the entire CONUS is 0.2, with the lowest of 0.15 (for Herbaceous), and a highest RMSE of 0.3 (for Wetland). Once converted and categorized as below normal, normal, and above normal, the prediction has relatively high overall accuracy, showing a 75% average for the CONUS, 49% minimum (for Wetland), and 82% maximum (for Herbaceous) overall accuracy in August 2013. There is some spatial variability to the overall accuracy with slightly higher values in the Western half of the US, as compared to the Eastern half. For January 2013, the average RMSE for the entire US is 0.35, with the minimum and maximum of 0.30 (for Wetland) and 0.42 (for Herbaceous), respectively. The overall accuracy of January 2013 is 45% for the entire US, with 37% minimum (for Herbaceous) and 55% maximum (for Wetland) overall accuracy.   To investigate the model performance over the entire time-series, we calculated the RMSE and overall accuracy for the 2004-2013 period. The time-series of average RMSE and average overall accuracy (Figure 9) shows that when RMSE values are low, overall accuracy is high and when RMSE values are high, the overall accuracy is low. In particular, during the winter months of each year, there are less fires by which to develop the models and thus higher uncertainty and lower accuracy. RMSE values range between 0.19 and 0.42, while overall accuracy ranges between 34 and 75 percent.

Burned Area
On average across the entire time-series, we have an average overall accuracy of 53 percent and RMSE of 0.31. This average denotes the ability of the models to capture both the spatial and temporal patterns of burned area anomaly in the United States over the record 2004-2013. To investigate the model performance over the entire time-series, we calculated the RMSE and overall accuracy for the 2004-2013 period. The time-series of average RMSE and average overall accuracy (Figure 9) shows that when RMSE values are low, overall accuracy is high and when RMSE values are high, the overall accuracy is low. In particular, during the winter months of each year, there are less fires by which to develop the models and thus higher uncertainty and lower accuracy. RMSE values range between 0.19 and 0.42, while overall accuracy ranges between 34 and 75 percent. On average across the entire time-series, we have an average overall accuracy of 53 percent and RMSE of 0.31. This average denotes the ability of the models to capture both the spatial and temporal patterns of burned area anomaly in the United States over the record 2004-2013. To better investigate the performance of the model in each land cover type, we have further shown the evaluation statistics in Table 1 and Supplementary Tables S1-S6. We calculated mean and standard deviation (SD) of RMSE for probabilistic wildfire danger model and mean and standard deviation of overall accuracy for the categorical wildfire danger model. It can be revealed that RMSE for the entire CONUS is 0.31 ± 0.05 and overall accuracy for the entire CONUS is 52.8 ± 8%. The underestimation and overestimation errors are 23.6%. The results indicate that the model performance does not vary much between land cover types. Wetland shows slightly higher RMSE and lower overall accuracy compared to other land cover types (0.35% and 45.5%, respectively). This is likely due to the considerably smaller sample size of historical Wetland wildfires compared to other land cover types (See Table 1). We have also shown the overall accuracy of the model for each land cover type (Supplementary Tables S1-S5) and Supplementary Table S6 Tables  S1-S6 are shown in Table 1. As indicated, the model shows a generally similar underestimation and overestimation error. The highest difference between the underestimation and overestimation errors is for the Deciduous land cover type (5.5%). This indicates that the model does not show significant systematic bias toward underestimation or overestimation. The correlation between 2004-2013 CONUS RMSE and overall accuracy is 0.95. This relationship between overall accuracy and RMSE can be quantified as: where RMSE is the RMSE at month t, OA t is overall accuracy at month t. a and b are constants. The goodness of fit (R 2 ) is 0.89. The strong relationship between RMSE and overall accuracy may be used to estimate the model performance using only one metric. We emphasize, however, that each of these metrics present a unique measure of model evaluation. RMSE is a measure of accuracy for the probabilistic wildfire danger model, while overall accuracy is a measure of accuracy for our categorical wildfire danger model.

Discussion
Our model (FDEO) relies on a single model for each land cover type, despite the variant behavior of historical fire occurrence throughout the year (Figure 2), which likely contributes to the uncertainty in predictive skill (Figure 9). This is due to the fact that we use satellite observations from two months prior to fire occurrence to build our fire danger models by land cover class (Figure 3). For example, for all fires in the deciduous land cover class, we used VPD, SSM, and EVI from two months prior in the location of each burn-similar in methodology to multiple studies ( [21,38,40]). By doing this, the model is temporally consistent, i.e., the same land cover model works for predicting a fire burning in July from one that burns in October. This could explain some of the uncertainty observed nationally (Figure 9), as there may be contrasting controls on fire within the same land cover class depending on the time of year of the burn [41,42]. Further contributing to the uncertainty in predictive skill is the assumption that fire regimes are similar in different vegetation classes, which over-simplifies the complexity of the fire regime concept [43,44] and neglects the variability of fire regime within a given land cover class [45,46].
Despite the crude characterization of fire regime by land cover class, the ability of explanatory variables to predict by land cover was relatively robust (R 2 ≥ 0.66). Previous studies show that fire occurrence is the eastern deciduous forests is due to a combination of climate (stressed conditions) and non-climate (elevated non-native insects and pathogens) factors [47]. Looking at Figure 3, it is revealed that fire occurrence in deciduous forests is related to prior months' stressed hydrologic conditions in the sense that stressed soil moisture conditions are linked to higher wildfire danger. Previous research show that herbaceous, evergreen forests, and wetlands experience larger wildfire danger the more stressed the vegetation is. In Western United States, the wildfire danger in the mountainous higher-elevation forested regions (Rockies, Cacades, and Sierra) and the herbaceous Great Plains is primarily related to the prior-month's stressed hydrologic conditions [36]. In the southeastern evergreen pine forests, wildfire risk is associated with droughts and the flammability of dead and live fuels are expected to increase with increasing magnitude and frequency of droughts [48,49]. Our results also show that the fire danger in evergreen forests, herbaceous lands, and wetlands is directly related to prior-month drought conditions. The fire danger in evergreen forests and herbaceous lands is related to VPD stressed drought conditions, while wetlands' fire danger is linked to SSM drought conditions. On the other hand, previous research indicates that the fire danger in the shrubland-dominated southwestern arid and low-elevations regions of the Western United States (Arizona-New Mexico Mountains Semidesert, American Semidesert and Desert, Chihuahuan Semidesert, Colorado Plateau Semidesert, and Intermountain Semidesert) is more associated with the vegetation growth and greenness in the prior months [36]. Our results also support prior literature findings. According to Figure 3, the shrubland wildfire danger is dominantly associated with prior months' EVI and hence the vegetation greenness.
The relatively high predictive ability of the models by land cover class results in robust spatial characterization of burned area anomalies (Figure 4), which improves (Figures 6 and 7) after bias adjusting the classification of below normal, normal, and above normal ( Figure 5). The bias-adjustment ( Figure 5) is based on the assumption of similar fire regimes across each land cover type and is achieved by setting thresholds of below and above normal based on the cumulative distributions of observation and prediction anomalies in each of the land cover types. This method was selected primarily due to the limitations of year-long wildfire activity in each of the 0.25 degree pixels. For example, historical fire occurrence (Figure 2) shows that there is little/no fire activity in the winter months of December, January, and February in western CONUS. Thus, aggregating the probability of area burned by land cover enabled classification of below normal, normal, and above normal. It is possible that the thresholds in probability vary regionally ( [21,38]), in which case a bias-adjustment based on land cover by geographic region might provide alternate results. Further research needs to explore generating probability and categorical fire prediction maps from a regional point of view and to compare the fire prediction with our land cover based model.
Looking at the results and the evaluation metrics (Figure 8), the developed model shows high accuracy in predicting the fire danger across the CONUS (average overall accuracy of 52 percent and average spring and summer overall accuracy of 57 percent). Model performance results in lower RMSE and higher overall accuracy in spring and summer months (Figures 8 and 9). The relative better performance of the model in the spring and summer months is likely due to the higher number of actual historical fire occurrences in those months, which consequently bias the model calibration to higher accuracy in those months. Spring and summer months show low uncertainty (RMSE) and high overall accuracy, while fall and winter months show high uncertainty (RMSE) and low overall accuracy. We anticipate that the model performance will improve as we obtain more data for training and the fire season lengthens [50]. For example, the RMSE and overall accuracy maps from August and January of 2013 ( Figure 8) reveal that the model performs with more accuracy in the regions with higher historical fire activity (Figure 2). For August 2013, lower RMSE and higher overall accuracy is seen in the western parts of CONUS, while for January 2013, we see more accurate model performance in southeast parts of the CONUS. Another observation from Figure 8 is the relatively higher model performance (lower RMSE and higher overall accuracy) in the western parts of the CONUS during August 2013. This is likely also due to higher historical fire activity in the western parts of the CONUS in the summer (Figure 2). The higher model performance in the western parts of the CONUS is critical since major large fires occur in the area and such fires result in billions of dollars in suppression and major long-lasting losses, including property damage, degradation of habitat, and air quality reductions [38,49,51]. Further contributing to the model error, other explanatory variables such as soil temperature conditions may lead to higher accuracy of the models during the winter period [52]. Frozen surface conditions may lead to less fires during winter and in freeze/thaw transition periods. Future research needs to examine the potential of soil moisture conditions for predicting wildfire danger during the winter months.
Another limitation of our model is that it is a relatively coarse resolution of 0.25 degree. Since fire danger models are dependent on land cover types, the coarse resolution of the land cover map could have impacted the results. The 2011 NLCD land cover map [35] is originally at 30-m resolution.
However, we re-gridded the map into 0.25 degree using the most popular land cover in each 0.25-degree pixel to match the spatial resolution of the other satellite data inputs. The re-gridded land cover map may have impacted the results as well as RMSE and overall accuracy, especially in the southeastern of the US, where there is great land cover heterogeneity within a 0.25 degree [35]. Future work is needed to determine the sensitivity of the models to the spatial heterogeneity of land cover.
In addition to earth observations, humans play a major role in wildfire occurrence. Human-started fires accounted for 84% of all fires and 44% of total area burned in the United States from 1992 to 2012. Majority of wildfires in the western and southeastern parts of the United States are human-caused, whereas lightning-started fires dominated the mountainous regions of the western United States [34,53]. Previous studies show that anthropogenic-caused wildfires have been linked to a combination of human activity and suitable climatic conditions [54,55]. While climatic conditions play a major role in detecting human-caused fires, they may not correctly predict all human-caused fires. This could also limit FDEO model capability to accurately forecast fire danger.
In meeting our objective to develop an algorithm for two-month lead prediction of fire danger for use in operational fire management resource allocation, we verified our spatially distributed, quantified prediction of fire danger against the expert-derived NIFC Wildfire Potential Outlook ( Figure 10). We investigated the maps of observation, prediction, and NIFC categorical prediction maps in August 2013. NIFC generates one to four months' significant wildland fire potential prediction maps for the US every month. To generate these maps, NIFC convenes regional experts who use various meteorological inputs to inform how to draw perimeters of expected fire danger classified as below normal, normal, and above normal. Because expert knowledge is inherently subjective and varies by region, it is difficult to compare "above normal" fire danger in one region to another. While expert knowledge have shown to help wildfire prediction, it is associated with a wide range of limitations, including careful matching of experts to application, lack of transparency, and repeatability, determining how to best aggregate multiple expert responses and linguistic uncertainty [56,57].The spatially-distributed deterministic approach presented in this paper offers a means for objectively classifying below normal, normal, and above normal such that regional comparisons and demand for fire management resources can be compared. As seen in the figure, observation and prediction maps (left and middle plots) show consistent patterns across the US, while the NIFC prediction map does not match the observation well. In the south-eastern parts of the US, all three maps show below-normal fire activity. While NIFC Wildfire Potential Outlook maps are derived from meteorological forecasts, FDEO is dependent on direct satellite observations. With overall high performance of our framework in both probabilistic and categorical fire danger forecast, FDEO is a promising framework for assessing wildfire danger in the entire CONUS, which relies on near real-time remote sensing observations and leverages lead time desired for resource allocation management. Further research needs to assess the performance of FDEO in real-time wildfire danger prediction and in decision-making processes with collaboration with NIFC and USFS. management. Further research needs to assess the performance of FDEO in real-time wildfire danger prediction and in decision-making processes with collaboration with NIFC and USFS.

Conclusion
In this study, we introduced FDEO, which is a two-month lead wildfire danger model for potential operational fire management practice. FDEO uses remotely sensed VPD, SSM, and EVI to predict spatial patterns of wildfire potential in the contiguous United States. While previous studies have either used meteorological input for wildfire predictions or looked at the relationships between satellite data with wildfire activity in different regions of the United States separately, this study introduced a spatially-distributed framework for predicting wildfire danger in the entire CONUS using remote sensing data. In addition, this study focuses on two-month lead predictions of wildfire danger, which is a suitable lead time needed for operational resource allocation, as opposed to the seasonal or daily lead predictions of most studies. The two-month lead time was selected to account for data latency of satellite observations. The results show that the model has high accuracy in predicting fire danger, especially in the spring and summer months, and offers a promising method for producing two-month forecasted fire danger maps. Had we had more wildfire data record, the model could potentially be improved in the non-fire season (fall and winter months). Future work should expand the analysis to a global scale and test near real-time deployment of these algorithms. As our model relies only on remotely sensed inputs, the global prediction of wildfire potential could be easily achieved. Finally, this can potentially be used for operational fire management in remote regions, where access to ground-based observations is limited. Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1,

Conclusions
In this study, we introduced FDEO, which is a two-month lead wildfire danger model for potential operational fire management practice. FDEO uses remotely sensed VPD, SSM, and EVI to predict spatial patterns of wildfire potential in the contiguous United States. While previous studies have either used meteorological input for wildfire predictions or looked at the relationships between satellite data with wildfire activity in different regions of the United States separately, this study introduced a spatially-distributed framework for predicting wildfire danger in the entire CONUS using remote sensing data. In addition, this study focuses on two-month lead predictions of wildfire danger, which is a suitable lead time needed for operational resource allocation, as opposed to the seasonal or daily lead predictions of most studies. The two-month lead time was selected to account for data latency of satellite observations. The results show that the model has high accuracy in predicting fire danger, especially in the spring and summer months, and offers a promising method for producing two-month forecasted fire danger maps. Had we had more wildfire data record, the model could potentially be improved in the non-fire season (fall and winter months). Future work should expand the analysis to a global scale and test near real-time deployment of these algorithms. As our model relies only on remotely sensed inputs, the global prediction of wildfire potential could be easily achieved. Finally, this can potentially be used for operational fire management in remote regions, where access to ground-based observations is limited.