Combining precipitation forecasts and vegetation health to predict fire risk at subseasonal timescale in the Amazon

Current forecast systems provide reliable deterministic forecasts at the scale of weather (1–7 days) and probabilistic outcomes at the scale of seasons (1–9 months). Only in recent years research has begun transitioning to operational settings to provide numerical predictions for a lead time of 2–4 weeks, a timescale known as subseasonal. The Subseasonal Experiment (SubX) multi-model ensemble mean precipitation forecast (2017–2021) for days 8–14 (week-2 forecast) is used as a covariate in logistic regression models to predict fire risk in the Amazon. In a complementary experiment, a vegetation health index (VHI) is added to SubX precipitation forecasts as a predictor of fires. We find that fire risk can be skillfully assessed in most of the Amazon where fires occur regularly. In some sectors, SubX week-2 precipitation alone is a reliable predictor of fire risk, but the addition of VHI as a predictor results both in (a) a larger portion of the Amazon domain with skillful forecasts and; (b) higher skill in some sectors. By comparing two sectors of the Amazon, we find that the added information provided by VHI is most relevant where the mosaic of land covers includes savannas and grassland, whereas SubX precipitation can be used as the sole predictor for week-2 fire risk forecast in areas where the mosaic of land cover is dominated by forests. Our results illustrate the potential for using numerical model forecasts, at the subseasonal timescale, in combination with satellite remote sensing of vegetation to obtain skillful fire risk forecasts in the Amazon. The operationalization of the methods presented in this study could allow for better preparedness and fire risk reduction in the Amazon with a lead time greater than a week.


Introduction
The use of statistical methods to transform numerical model outputs into calibrated seasonal predictions is a well-established field in climate research with applications in agriculture (Esquivel et al 2018, Shelia et al 2019, Fernandes et al 2020, Pons et al 2021, water resources (Muñoz et al 2010, Robertson et al 2014, Najafi et al 2021 and public health (Muñoz et al 2017, Landman et al 2020, Colón-González et al 2021. In ecosystems management, one type of application consists of combining statistical methods with numerical predictions to assess or forecast seasonal fires (Aragão et al 2007, Murdiyarso and Adiningsih 2007, van der Werf et al 2008, Field et al 2009, Fernandes et al 2017. Sea surface temperature (SST) in the Atlantic and the Pacific basins are good indicators of anomalous precipitation, at various timescales, in the Amazon (Fu et al 2001, Espinoza et al 2013, Fernandes et al 2015, Marengo and Espinoza 2016. Thus, SST forecast and observations can be used to predict fire occurrence in the Amazon (Chen et al 2011, Fernandes et al 2011 and other fire prone regions (Ceccato et al 2014, Chen et al 2016, Turco et al 2018 with a lead time of 1-4 months. At the other end of the prediction timescale lies fire weather forecasts (1-7 days), varying in spatial coverage from global (Field et al 2015, Di Giuseppe et al 2020 to regional, where models are calibrated to better represent the environmental characteristics of the region (IDEAM 2015, Aliaga Nestares et al 2018, Setzer et al 2019, Martins et al 2020. While both timescales provide decision-makers with much needed information, seasonal forecasts can be too long a horizon for the positioning of resources and the fire weather can be too short. In between lies the 2-to-4 week subseasonal timescale, for which numerical prediction has recently been established by research-to-operation efforts such as the S2S Prediction (Vitart et al 2017) and the Subseasonal Experiment (SubX) (Pegion et al 2019).
We use the SubX model outputs, given the statistically significant skill of precipitation forecast over all months in South America and the facilitated data retrieval and manipulation offered by the International Research Institute for Climate and Society (IRI) Data Library (Blumenthal et al 2014).
We develop a hybrid (dynamical-statistical) model using the NextGen prediction system methodological approach (Munoz et al 2019, WMO 2020 to predict fire risk. NextGen involves designing, calibrating, building ensembles, and verifying objective forecasts by first identifying decision-relevant variables followed by the analysis of the physical mechanisms, sources of predictability and suitable candidate predictors. Here, the decision-relevant variable is fire risk, defined as the likelihood of active fires occurrence as a function of precipitation and vegetation health conditions. This definition might be considered 'event risk' in disaster or hazards literature (Sarewitz et al 2003), as it does not explicitly consider the vulnerability of specific communities to fire. The Visible Infrared Imaging Radiometer Suite (VIIRS) (Hillger et al 2013, Schroeder andGiglio 2017) active fires are modeled as a function of the SubX multimodel ensemble mean (MME) week-2 precipitation forecast, for the 2017-2021 period. A complementary experiment includes the satellite-based NOAA Vegetation Health Index (VHI, Kogan 1997Kogan , 2001 as an additional predictor. VHI is updated in real-time and it reliably detects ecosystem sensitivity to meteorological conditions at various timescales (Andujar et al 2017. We then take a closer look at the skill of the statistical models in two sectors of the Amazon, each representative of a distinct fire season. Lastly, we explore the role of contrasting land covers on fire risk models' performance in each of these two sectors.

VIIRS active fires (VNP14IMG)
Fire occurrence is estimated using the active fire product from the VIIRS instrument on board of the Suomi National Polar-orbiting Partnership satellite. The daily VIIRS product at a resolution of 375 m (VNP14IMG) consists of a fire mask with ten classes. Three identify low, nominal and high confidence levels of active fire detection while the other seven are classes unrelated to fires, such as water and clouds (table 1 of Schroeder and Giglio 2017). We resampled all variables to an intermediary grid of 0.22 • latitude-longitude resolution within which, every active fire of nominal and high confidence is counted. The resampled daily active fires are then averaged to other timescales as needed (sections 2.5.1 and 2.5.2). Although the focus of our study is the Amazon biome, the full domain of the analyses comprises of tropical South America (81.25 • W-44.25 • W and 19.5 • S-11.5 • N). The VIIRS data is available from January 2012 to current.

Subseasonal precipitation forecast (SubX)
SubX (Kirtman et al 2017, Pegion et al 2019 consists of a multimodel publicly available database of historical re-forecasts and real-time forecasts. The historical re-forecasts start in 1999 and end between 2014 and 2016 depending on the model. Real-time forecasts are available from six groups and include forecasts for at least 32 days beyond the initialization date (table 1). The real-time forecasts start at various times in 2017. By the 5 August initialization date, five of the six models produce forecasts except for NCEP-CFSv2 which starts on 15 October. All models maintain up-to-date forecasts on a weekly (or more frequent) basis that can be accessed through IRI Data Library (Blumenthal et al 2014). Although combined historical re-forecast and real-time forecast provide data for over 20 years, we use only output from the real-time forecasts, beginning in 2017, as the number of ensemble members in each model differs from historical re-forecast to real-time periods and a multi-model mean calculated over 1999-2021 would be inconsistent. The SubX models' spatial resolution of 1 • latitude-longitude, was linearly interpolated to 0.22 • latitude-longitude grid to match the resolution of aggregated VIIRS data (section 2.1). Spatially downscaling numerical model outputs can introduce noise in interpolated precipitation (Wood et al 2004), especially during the wet season. This should be less of an issue during the dry season. To test this assumption, we compared the distribution of a sample consisted of SubXPr multimodel mean week-2 forecast at its original resolution (1 • latitudelongitude) with the linearly downscaled resolution (0.22 • latitude-longitude) sample during one dry season in southern Amazon. Both samples include every gridcell within the domain of coordinates 76 • W-57 • W and 13 • S-4 • S and the 13 weeks of the 2018 dry season (July-September). Using a two-sample Kolmogorov-Smirnov test, we tested the null hypothesis that both samples came from the same continuous distribution. The test result does not reject the null hypothesis at the 1% significant level.
The SubX real-time precipitation forecast (Sub-XPr) is arranged to reproduce 52 week-2 (days 8-14) forecast period per year by: (a) Calculating, for each SubX model, the ensemble members mean week-2 precipitation forecast corresponding to 52 annual non-overlapping weeks, and a final sample of 208 elements per model (52 week-2 forecasts per year times 4 years. Mid August 2017-mid August 2021). See table S1 for details. (b) The MME mean consist of week-2 forecasts averaged among the six SubX models.
Constructing a multi-model mean of week-2 realtime forecasts is a challenging task. The different models are initialized on different days of the week and at different intervals. To take advantage of Wednesday initializations in three of the six models, we fix non-overlapping week-2 real-time forecasts (average of forecast days 8-14) from Thursday (day 8) to Wednesday (day 14). For the other three models, week-2 real-time forecast varies from average forecast days 9-15 to average forecast days 12-18 depending on the SubX model and the forecast week (table 2).

Vegetation health index (VHI)
Fire ignition in the Amazon results from human activities, but the intensity and spread can be regulated by natural processes, such as vegetation water stress (Asner and Alencar 2010). One indicator of vegetation stress is the satellite-based NOAA-STAR VHI calculated from the vegetation condition index (VCI), a proxy for vegetation moisture, and temperature condition index (TCI), a proxy for thermal effects (Kogan 2001) VHI is available at 4 km spatial resolution, which was also re-gridded to the 0.22 • latitude-longitude resolution. Near real-time weekly-averaged VHI is released on the day of the week corresponding to 8 January of each year. In 2018, for example, the product is updated every Monday (table 2) for a total of 52 global maps per year. For VHI to be used as a predictor, the data needs to be available by the weekly initialization date of the SubX multi-model mean precipitation forecast, and this varies within the years of study. From 2012 to current, the indices are calculated using the VIIRS instrument (Kogan et al 2015. Table 2. Example of initialization dates used to produce a multi-model mean week-2 forecast for target date 18-24 January 2018, shaded blue in the table. The green cells mark initialization dates for each model. The X marks the Wednesday (10 January) initialization for a week-2 (average of forecast days 8-14) forecast and the closest prior dates for models with initializations on days other than Wednesday. The actual lead time for each forecast is shown below the target week-2 forecast period. The dates of observed VHI and VIIRS active fires corresponding to the example SubX precipitation target week are also marked in the The 2017-2021 VHI is obtained for dates corresponding to the available SubX week-2 precipitation forecasts.

Land cover-Mapbiomas
We use Mapbiomas Land Cover Collection-6 Brazil (Souza et al 2020) in two areas of the Amazon chosen for a closer evaluation of fire risk model skills. We chose the year 2019 land-cover as representative of the 2017-2021 period.

Active fire season
The analysis relating precipitation, fires and vegetation health is conducted for the peak fire season. The nearly 10 years of VIIRS data (January 2012mid-August 2021) is used to determine the climatological most active fire trimester at each gridcell by calculating monthly averages of daily VIIRS active fire data (section 2.1) and then calculating the climatology for the annual 12 overlapping trimesters (JFM, FMA, etc). The trimester with highest fire activity (FA) is then identified (figure 1).

Weekly data
Once the most active fire trimester is identified, the VIIRS, SubXPr and VHI data for the 13 weeks within that trimester are fitted in a logistic regression model. If a gridcell's most active fire trimester is September-November, data for the season's 13 weeks of years 2017, 2018, 2019 and 2020 (total of 52 weeks = 13 weeks times 4 years) is fitted in the statistical model. For March-May, for example, the data sample contains the trimester's 13 weeks of 2018, 2019, 2020 and 2021 given the SubX multi-model forecast is available from August 2017 to August 2021.

Statistical analysis of fire probability
Fires are used in the Amazon as a tool to clear debris from deforestation and to manage agropastoral land.
They can occur concentrated in a few weeks as burning is planned for times when vegetation is dry. Thus, the timeseries of weekly FA may show a mix of high activity interspaced by no fires. In the majority of our study domain, weeks with no fires occur at least during 20% of the time in the dry season (figure S1) and a two-mode behavior (fire/no-fire) is then fitted to predictors using logistic regression.
To assess the potential for predicting weekly FA, we choose logistic regression expressed as: where P(FA = 1) is the probability of fire occurrence (FA = 1) modeled as a linear combination of predictors X n with slopes β n and an intercept α. The probability of fire occurrence is modeled as a function of non-overlapping week-2 SubXPr and VHI. The FA binary response is coded 1 for values greater than 0 and 0 for no fire within a 0.22 • latitude-longitude gridcell.
The logistic regression models were fit using Matlab 2020a fitglm function (MathWorks 2020) for a binomial distribution. In one experiment, FA was fit to SubXPr and in the other FA was fit to Sub-XPr and VHI. The experiments are henceforth called FA_f (SubXPr) and FA_f (SubXPr&VHI), respectively. The logistic regression models were crossvalidated using the leave-one-out option, which corresponds to leaving out one week of paired datasets at each iteration. The assessment of gridcell's logistic regression model fit is evaluated by determining what portion of observed positives is classified as positive (sensitivity) and what portion of observed negatives is classified as negative (specificity). A plot of sensitivity against 1-specificity describes the receiver operating characteristics (ROC, Wilks 2011) and goodness of fit can be assessed as the area under the ROC curve (AUC). An AUC of 0.5 represents discrimination expected by chance to predict positive outcomes. AUC above 0.5 mean that the positive outcome has a higher predicted probability than the negative subject, and the logistic regression model is considered skillful. A value of 1 stands for a perfect model.

SubXPr and VHI as predictors of fire probability
We discuss our results for a subset of gridcells with intermediate and high values of AUC (>0.6). In some sectors, using SubXPr week-2 forecast as the sole predictor (FA_f(SubXPr) results in AUC > 0.7, notably in the southwestern Brazilian Amazon marked with a box in figure 2(a) and in eastern equatorial Amazon. AUC values show improvement for experiment FA_f (SubXPr&VHI), especially in the southern Amazon where an extensive area shows AUC > 0.6 ( figure 2(b)). An exception is the eastern equatorial Amazon, where precipitation alone results in better performance both in absolute values and in number of gridcells with high AUC.
The results in figure 2 illustrate the potential for relying on SubX multi-model mean precipitation forecasts to determine the probability of fire occurrence in the Amazon with a lead of 8-14 days (week-2). Moreover, by including the state of vegetation as an indicator of fire risk, the logistic regression models perform better over a larger sector of tropical South America.

Model performance in southern and northern Amazon
To further explore fire occurrence response to SubXPr and VHI, we selected two sectors of the Amazon with a cluster of high AUCs for each of the experiments. In the southern Amazon (73.6 • W-68.5 • W, 9.8 • S-6.5 • S), where the main climatological fire season is August-October, high AUC is found for experiment FA_f(SubXPr) (figure 2(a)). In the northern sector (63 • W-59.2 • W, 1.7 • S-4.8 • N), where the climatological fire season is January-April, a cluster of high AUC is found for experiment FA_f (SubXPr&VHI) ( figure 2(b)).

Forecast time series
We calculate, for each spatial domain, how many gridcells predicted the correct outcome among those with AUC > 0.6. A successful fire forecast means fire probability higher than 0.5 in a gridcell that detected fire or fire probability below or equal 0.5 in a gridcell that detected no-fire. The gridcells with a 'hit' are counted and a ratio between the number of correct hits and the total number of gridcells with AUC > 0.6 within the domain is calculated for every forecast week (figure 3). Given the first SubXPr forecast is available for the week of 17-23 August 2017, the weeks used in the cross-validated models within the southern domain (fire season August-October) start in  (table S1).
Within the southern domain the number of gridcells with AUC > 0.6 are 154 and 161 for FA_f (SubXPr) and FA_f (SubXPr&VHI), respectively. This similar number of skillful models can be observed in figures 2(a) and (b), where the gridcells with AUC > 0.6 are shaded. Collectively, both experiments result in similar performance with an average of 75.6% of gridcells correctly predicting fire risk in experiment FA_f (SubXPr) and 74.3% in the FA_f (SubXPr&VHI) experiment ( figure 3(a)).
In the northern domain the number of gridcells with AUC > 0.6 are 76 and 127 for experiments FA_f (SubXPr) and FA_f (SubXPr&VHI), respectively. On average, 74% of gridcells correctly forecast outcomes for the FA_f(SubXPr) experiment and 69% for FA_f(SubXPr&VHI) ( figure 3(b)). Relatively, the models in experiment FA_f(SubXPr) outperform those of experiment FA_f(SubXPr&VHI), but the latter represents a considerably larger overall area-more gridcells-with skillful forecasts. The average 74% correct hits from FA_f(SubXPr) corresponds to 56 gridcells, whereas the 69% from FA_f(SubXPr&VHI) corresponds to 87 gridcells, indicating that the addition of VHI as predictor provides a more useful set of information in the northern domain, especially for extreme VHI values (see figure 4 discussion). In grasslands, common in the northern domain (section 3.2.3), fires occur regularly, and fire risk will be mainly determined by whether it rains or not. However, if at the start of and during the dry season VHI is 'neutral to highwetness' , VHI becomes a valuable indicator of firerisk as absence of rain will not immediately translate in higher fire-risk given the high-wetness state of vegetation. The high VHI state in 2021 (figure S2) would explain the models' performance in that year ( figure 3(b)), when SubXPr alone did a poorer job in predicting fire risk than SubXPr and VHI together.

Fire risk sensitivity to SubXPr and VHI
To test fire risk sensitivity to SubXPr and VHI in two domains identified in section 3.2, we fitted a logistic regression model to one continuous sample comprised of data from gridcells within each domain and all 52 weeks. Our objective is to evaluate if one single logistic regression model fitted to the aggregated data from each domain exhibits characteristics similar to the gridcells' individual models, namely low sensitivity to VHI in the southern domain and high sensitivity in the northern domain. The southern domain has the spatial dimension of 25 (longitude) per 17 (latitude) cells for a total of 425 at each timestep. A sample derived from this domain contains a maximum of 22 100 elements (425 gridcells times 52 weeks). In the northern domain the spatial dimension is 18 (longitude) per 32 (latitude) resulting in 576 gridcells at each timestep for a maximum of 29 952 elements (576 gridcells times 52 weeks). The sample sizes for the southern and northern sectors contain in fact 13 185 and 14 201 elements, respectively, once all the gridcells with an undefined fire season (figure 1) and occasional missing data are removed. We then reproduce the experiment FA_f(SubXPr&VHI) using the aggregate covariates for each domain and present Figure 3. Ratio between the number of gridcells with correct week-2 fire risk forecast and the total number of gridcells for which a skilfull model exist (AUC > 0.6) within the southern (a) and northern (b) domains. A correct outcome means the logistical model forecasted a probability greater than 0.5 for fire occurrence in gridcells where fires were detected and lower or equal to 0.5 for gridcells with no-fire.  the logistic regression models' coefficients, standard error, statistical significance, and AUC in table 3. The SubXPr and fire relationship sensitivity to VHI is evaluated from setting two extreme VHI scenarios. Low VHI is the average of the 10% lowest values in the aggregated time series used to fit the statistical models and high VHI is the average of the 10% highest values. This is done for both domains using equations (5) and (6) Fire probability as a function of SubXPr for low VHI (vegetation stress) shows higher values in general than that describing fire probability for high VHI (no-stress/high-wetness) conditions in the southern domain ( figure 4(a)). For precipitation of 3 mm d −1 there is a 45% chance of observing fires under high VHI conditions, increasing to 56% for low VHI, indicating low fire sensitivity to vegetation conditions regardless of precipitation amounts.
In the northern domain, a more robust shift of probabilities is observed, especially for lower values of precipitation (<7 mm d −1 ) as seen figure 4(b), describing a prominent role of vegetation health in determining fire risk. At 3 mm d −1 , fire risk goes from 43% for high VHI to 69% for low VHI. Beyond 10 mm d −1 , the two VHI scenarios converge to describe a dominant role of precipitation on fire risk regardless of the state of vegetation.

Land cover
These two sectors of the Amazon are characterized by different land covers as estimated by the Mapbiomas Land Cover classification for year 2019 (figure S3). We resampled the 30 meter resolution land cover map (Souza et al 2020) to the VHI 4 km original grid, by retaining the class with the highest frequency to determine the land covers most influential for the VHI estimate. The southern domain contains 98% forest and nearly 2% pasture (figure 5(a)), whereas in the northern domain vegetation is more varied and comprised of 77% forest, 20% grassland and 1.6% pasture ( figure 5(b)).
The vegetation distribution in the two sectors, provides a plausible explanation for the distinct fire risk sensitivity to VHI. In VHI, greenness normalized vegetation index (NDVI) is as a proxy for vegetation moisture conditions (equation (1)). However, reduced cloud cover and more light toward the end of the dry season can result in a green-up (higher NDVI) of the Amazon forest (Huete et al 2006). At the spatial scale of VHI (4 km), this green up could mask NDVI of other land covers within a gridcell, reducing the power of VHI as an indicator of fire risk in land covers that are more susceptible to burning but border green forest patches. Even if VHI can describe forest water stress, that does not necessarily mean susceptibility to fires. Understory humidity and proximity to fire sources, also play an important role in determining whether fires occur in forest land covers.
In contrast, savannas and grasslands are more sensitive to short-term water stress and VHI responds to meteorological drought, or recovery from it, within a couple of weeks to a month (Sehgal et

Discussion
While atmospheric scientists focus on the drivers of fires, fire ecologists focus on the atmospheric demand for moisture (e.g. vapor pressure deficit) and the potential for vegetation to become fuel (Fu et al 2021).
Our study aims at reducing this gap by combining observed states of vegetation with precipitation forecasts. Using a SubX multimodel ensemble mean of week-2 precipitation forecasts, we find that fire probability can be skillfully predicted in large portions of the Amazon. By pairing week-2 SubXPr with VHI to predict fire probability, we find areas in which logistic models' skills improve upon those using only week-2 SubXPr forecast as a predictor. It is also evident that a greater number of gridcells in the domain of study present skillful statistical forecast for the experiment using both VHI and SubXPr as predictors of fire risk. This is especially true in the Brazilian arc of deforestation and in Bolivia that are characterized by land covers that are more susceptible to fires (grassland, fallow, and agropastoral land). In areas where land cover is dominated by forests, the role of vegetation health as a predictor of fire risk is modest. That is the case of the southern Brazilian Amazon, where fire risk is mainly driven by SubXPr and the probability of fire occurrence increases only about 10% for very low (water stressed) to very high VHI across different precipitation forecasts.
In contrast, in the northern Amazon where savannas and grasslands are more prevalent, VHI status helps determine fire risk in the weeks ahead. This is seen in the larger number of gridcell with skilfull fire risk forecast for experiment FA_f (SubXPr&VHI) compared to FA_f (SubXPr). The sensitivity test also indicates that fire risk increases by as much as 28% from high VHI (non-stressed) to low VHI for lower precipitation values indicating the more prominent role of vegetation status in determining fire risk.

Conclusions
Our study shows the potential for using precipitation subseasonal forecasts in combination with vegetation health in a hybrid dynamical-statistical model to forecast fire risk. The operationalization of the methods presented in this study have the potential to inform better practices for fire prevention and mitigation in the Amazon with lead time beyond one week, which is a crucial timescale for fire management resources allocation. Continuous improvement of statistical models' accuracy can potentially be achieved by using satellite vegetation conditions at higher spatial resolution, which would allow for better discrimination between forests and non-forest land covers, and by including more years of SubX precipitation forecast.

Acknowledgments
We acknowledge the agencies that support the SubX Project and the modeling groups for producing and making available their model output, including the National Centers for Environmental Prediction (NCEP), the NCEP Environmental Modeling Center, the National Aeronautics and Space Administration Global Modeling and Assimilation Office, the National Center for Atmospheric Research and University of Miami Rosenstiel School of Marine and Atmospheric Science, the Naval Research Laboratory, the National Oceanic and Atmospheric Administration Earth System Research Laboratory, and Environment and Climate Change Canada. We wish to thank Columbia University's International Research Institute for Climate and Society (IRI) Climate Data Library for disseminating the SubX data (https://iridl.ldeo.columbia.edu/SOURCES/.Models/. SubX/; DOI: 10.7916/D8PG249H).
K F and M B, wish to acknowledge the support of the SERVIR-Amazonia, a joint development initiative of the National Aeronautics and Space Administration (NASA) and the United States Agency for International Development (USAID). A G M was partially supported by NOAA Grant: NA18OAR4310275 and ACToday, the first Columbia World Project. We also thank the anonymous reviewers for the constructive comments and suggestions.