*Article* **Evaluating a New Relative Phenological Correction and the Effect of Sentinel-Based Earth Engine Compositing Approaches to Map Fire Severity and Burned Area**

**Adrián Israel Silva-Cardoza 1,†, Daniel José Vega-Nieva 2,†, Jaime Briseño-Reyes 2,\*, Carlos Ivan Briones-Herrera 2, Pablito Marcelo López-Serrano 2, José Javier Corral-Rivas 2, Sean A. Parks 2,3 and Lisa M. Holsinger 2,3**


**Abstract:** The remote sensing of fire severity and burned area is fundamental in the evaluation of fire impacts. The current study aimed to: (i) compare Sentinel-2 (*S2*) spectral indices to predict field-observed fire severity in Durango, Mexico; (ii) evaluate the effect of the compositing period (1 or 3 months), techniques (average or minimum), and phenological correction (constant offset, *c*, against a novel relative phenological correction, *rc*) on fire severity mapping, and (iii) determine fire perimeter accuracy. The Relative Burn Ratio (*RBR*), using *S2* bands 8a and 12, provided the best correspondence with field-based fire severity (*FBS*). One-month *rc* minimum composites showed the highest correspondence with *FBS* (R<sup>2</sup> = 0.83). The decrease in R<sup>2</sup> using 3 months rather than 1 month was ≥0.05 (0.05–0.15) for *c* composites and <0.05 (0.02–0.03) for *rc* composites. Furthermore, using *rc* increased the R<sup>2</sup> by 0.05–0.09 and 0.10–0.15 for the 3-month *RBR* and *dNBR* compared to the corresponding *c* composites. *Rc* composites also showed increases of up to 0.16–0.22 and 0.08–0.11 in kappa values and overall accuracy, respectively, in mapping fire perimeters against *c* composites. These results suggest a promising potential of the novel relative phenological correction to be systematically applied with automated algorithms to improve the accuracy and robustness of fire severity and perimeter evaluations.

**Keywords:** Sentinel-2; post-fire severity; initial fire assessment; soil burn severity; fire perimeter; image compositing; vegetation phenology

#### **1. Introduction**

The evaluation of burned area and fire severity is fundamental for understanding a variety of ecological processes (e.g., [1–4]), including fire emissions and carbon cycling (e.g., [5,6]), post-fire erosion (e.g., [7–9]), tree mortality (e.g., [10–12]), post-fire recovery trajectories (e.g., [13–16]), and ecosystem resilience (e.g., [17,18]). In spite of the great progress in fire severity and burned area evaluations using medium-resolution sensors, several challenges remain for the development of more standardized, automated methods of post-fire assessment. Among others, these challenges include: (i) the assessment of spectral indices, particularly for Sentinel-2, to map burned area perimeters and measure fire severity and fire effects by strata, namely vegetation and soil; (ii) the evaluation of the

**Citation:** Silva-Cardoza, A.I.; Vega-Nieva, D.J.; Briseño-Reyes, J.; Briones-Herrera, C.I.; López-Serrano, P.M.; Corral-Rivas, J.J.; Parks, S.A.; Holsinger, L.M. Evaluating a New Relative Phenological Correction and the Effect of Sentinel-Based Earth Engine Compositing Approaches to Map Fire Severity and Burned Area. *Remote Sens.* **2022**, *14*, 3122. https:// doi.org/10.3390/rs14133122

Academic Editors: Leonor Calvo, Elena Marcos, Susana Suarez-Seoane and Víctor Fernández-García

Received: 18 May 2022 Accepted: 21 June 2022 Published: 29 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

effect of compositing periods and techniques, particularly on automated platforms such as Google Earth Engine (e.g., [19]), on the accuracy for mapping initial assessments of fire severity and burned area perimeters; (iii) the quantification of the effect of non-fire-related temporal and spatial variations in spectral indices [20–22] and the development of simple, readily automated methods of phenological correction that account for the spatial and temporal heterogeneity of such effects [23].

#### *1.1. Evaluation of Spectral Indices to Map Fire Severity and Perimeter, Particularly for Sentinel-2*

Regarding this challenge, several spectral indices using medium-resolution sensors have been developed [3]. These include several formulations of *NBR* indices based on NIR/SWIR bands (e.g., [24–27]), and more recently, spectral indices including Sentinel-2 red-edge bands such as *CIre* [28,29] or *BAIS2* [30]. Nonetheless, the selection of the most appropriate indices from Sentinel-2 to map fire severity and perimeters is still an open research question, with a large variety of results depending on the ecosystems analyzed (e.g., [8,9,31–35]). Furthermore, while the majority of fire severity research has been developed in temperate to boreal ecosystems, there are comparatively fewer studies in areas with different climatic and vegetation characteristics. This is the case in Mexico, a very ecologically diverse country with high fire activity (e.g., [36–39]), where previous studies evaluating spectral indices, validated with field data, to map fire severity and burned area perimeters, are lacking.

#### *1.2. Evaluation of Compositing Techniques and Period for Fire Severity and Perimeter Mapping* 1.2.1. Comparison of the Effects of Compositing Techniques for Fire Severity and Burned area Mapping

Pixel-based compositing methods can provide a new analysis paradigm instead of depending on a single scene [40], resulting in an unprecedented opportunity for standardizing and automating fire severity and burned area evaluations [27,41–43]. Nevertheless, the image composite techniques analyzed in the literature have largely varied, depending on the study goal, location, and period analyzed. For example, mean composites have been used for extended assessments of fire severity [19,41,44,45] and mean, median, or medoid [46,47] composites have been also used to map burned area and recovery [47–50]. Likewise, minimum, or maximum composites [26,27,42,51–53] have been used for mapping fire perimeters and trends of burned area [54] or fire severity [55]. Because of this variety of approaches, there is a need for a formal evaluation of the effect of compositing techniques (i.e., average or minimum) and compositing periods on their accuracy in mapping fire perimeters and fire severity. Such an evaluation is particularly lacking for immediate postfire assessments of fire impacts (hereafter termed "fire severity" following [1], [2] or [54], which propose the use of this term for initial severity assessments), which are fundamental for post-fire erosion emergency planning (e.g., [8,9,56]) but have received less attention from compositing analysis studies.

1.2.2. Assessment of the Influence of Compositing Period for Mapping Initial Fire Severity and Burned Area Perimeter

The post-fire image assessment period can largely influence what is being observed (e.g., [1,57]). Non-fire variations in spectral indices can be caused by natural oscillations in vegetation and soil moisture (e.g., [58–60]) or vegetation processes such as senescence, regrowth, and phenological changes in grasses, shrubs, and trees (e.g., [61,62]), both for broadleaves (e.g., [63]) and conifer species (e.g., [64]). Therefore, the values of a spectral index can vary considerably from slight changes in image dates, resulting in substantially different burn severity maps [21,22,54]. This can affect image accuracy in mapping fire perimeters [62] and fire severity [19,58]. Moreover, it could also influence the widely reported lack of transferability of fire severity and burned area thresholds between seasons and ecosystems [20,65–67].

However, although evidence of variability in the assessment of burn severity driven by image selection has been shown in many studies, paradoxically, little attention has been paid to quantifying the temporal and spatial variations in spectral fire severity indices [20,21,54] and how this impacts their suitability for mapping fire perimeters [23] and severity [22]. In addition, in contrast to the more widely analyzed Landsat satellite, there are relatively few studies examining the temporal variations in Sentinel spectral indices' accuracy in mapping burned area [68] and fire severity [22]. Furthermore, most of those relatively scarce studies have used paired-image techniques, contrasting with a scarcity of studies addressing temporal and phenological variations in fire severity spectral indices' accuracy using composite images [19,22]).

#### *1.3. Development of a Spatially Variable Automated Phenological Correction*

In spite of the frequent documentation of the above-mentioned phenological biases, most of the literature has generally failed to include any measure to quantify or correct for such non-fire phenological changes to standardize fire severity assessments [20,69]. The most notable exception is the constant offset method [41,57,67,70,71]. This approach consists in subtracting from the uncalibrated *dNBR* image a phenological constant offset, generally calculated as the average *dNBR* from the unburnt area surrounding the fire perimeter [70,72,73]. Because of the relative simplicity of this method, some steps have been taken towards its automatization for Landsat [19,41,74] and Sentinel images [75]. Nevertheless, in order to be representative, the average offset should be calculated in similar vegetation types to those where the fire occurred [41,74,76]. This strong assumption presents a challenge in heterogeneous fuels, making the constant offset method not advisable in those cases, as it would be compromised by the spatial heterogeneity of phenological effects [41,45]. Still, fuel heterogeneity is very common (e.g., [23]), resulting in frequent spatial and temporal variations in the phenological effects (e.g., [54,77]).

Although some attempts to propose spatially variable techniques of phenological correction have been made (e.g., [22,54,77]), in their current states, the majority of those approaches require supervision to either manually select control pixels (e.g., [54,77]) or to manually create phenologically matched composites [22]. In this sense, to the best of our knowledge, the use of the pre-fire Normalized Burn Ratio [78] (*NBRpre*), a widely used proxy of biomass [70], has not been yet tested to automatically stratify and sample for spatial phenological variations in unburned vegetation. Such a systematic stratification of phenological effects could lead to a standardized, automated, novel relative phenological correction procedure for analyzing and minimizing the spatial heterogeneity of phenological variations. This approach could overcome some of the limitations of the offset method, while benefiting from its relative simplicity and ease of automation.

Consequently, the present work aimed to evaluate and map fire severity and burned area perimeter through paired and composite spectral indices of S2-MSI images, the latter calculated with different periods, techniques, and phenological corrections (constant offset, *c*, against a novel spatially variable relative phenological correction, *rc*). Specifically, the goals of the study were: (i) to evaluate the correspondence of paired Sentinel-2 (S2) images with field data from wildfires in the temperate to semi-arid forests of Durango, Mexico, to select the best S2 spectral index to map the initial assessment of fire severity in the analyzed ecosystems; (ii) to analyze the effect of Google Earth Engine S2 composite period length (1 or 3 months), technique (average or minimum), and phenological correction (*c* and *rc*) for estimating field-observed fire severity; and finally, (iii) to quantify the effect of composite period length, technique, and phenological correction on the accuracy in mapping the burned area perimeter for each wildfire.

#### **2. Materials and Methods**

#### *2.1. Study Area*

This study was conducted in the Mexican physiographic province of *Sierra Madre Occidental* in Durango State (SMO), between 26◦59 50N, 22◦20 26S latitude North and 103◦42 26E, 107◦14 12W longitude West. The study area is dominated by conifers of the genera *Pinus*, *Juniperus*, *and Cupressus*, as well as grasslands and shrublands dominated by *Acacia* and *Prosopis* species [79,80]. The climate ranges from semi-dry temperate to semicold subhumid and subhumid temperate according to the Köppen climate classification adapted for Mexico [81]. The average annual temperature ranges from 5 to 18 ◦C, with average annual precipitation ranging from 300 to 1500 mm [82]. The Durango SMO region is characterized by surfaces of large mesas that extend over the central part of the state, as well as mountain areas with slopes ranging from 20 to 80%. Altitude ranges from 800 to 3200 m [83]. We selected three fires that occurred between April and May 2019, covering a gradient in fuel types and climates ranging from semi-arid (Site 1) to humid (Site 2) and semi-humid (Site 3) (Table 1, Figure 1). Fire regimes in the study area include both highfrequency, low-severity fires in the denser, wetter pine and oak forests and lower-frequency, higher-severity fires in the more open, semi-arid forests [84]. Severe fires can affect species of timber value, such as pine stands, as well as tree and shrubby species used for charcoal production, such as *Quercus* or *Prosopis* species.

**Table 1.** Summary of wildfires selected in the Sierra Madre Occidental, Durango.


<sup>1</sup> [85], fuelbed types [84] (Figure 1): TU = grass or shrub understory mixed with timber litter, TL = timber litter, GS = grass and shrublands, and GR = grasslands.

**Figure 1.** Location of the three 2019 wildfires analyzed in Durango State, Mexico, where (**a**) = Site 1, (**b**) = Site 2, (**c**) = Site 3, AGAF = aggregated active fire perimeter [75], fuelbed types (scale: 1:250,000) [84]: 1.1.GR1 = warm–humid tall grass, 1.2.GR2 = temperate–humid short grass and shrublands, 1.3.GR3 = semi-arid short grass and shrublands, 2.1.GS1 = semi-arid grass and shrublands, 2.2.GS2 = arid grass and shrublands, 3.1.SH1 = arid tall shrublands, 3.2.SH2 = arid short shrublands, 4.1.TU1 = subhumid temperate timber understory with dense shrubs, 4.2.TU2 = humid temperate–cold timber understory with dense pasture, 4.3.TU3 = subhumid timber understory with sparse grasslands, 4.4.TU4 = subhumid warm to semi-arid timber understory with sparse grass and shrublands, 5.1.TL1 = humid temperate short-needle coniferous timber litter, 5.2.TL2 = humid long-needle coniferous timber litter, 5.3.TL3 = subhumid to humid sclerophyll broadleaf timber litter, 5.4.TL4 = seasonally dry humid broadleaf timber litter, 7.NB = not burnable (bare ground), 8.NB1 = not forestry areas (e.g., urban development or agricultural), and 8.1.NB2 = water bodies.

The selected wildfires included a heterogeneous-severity fire (low–medium severity to crown fire patches) (Site 1), with an estimated burned area of 23,809 ha (the largest wildfire in Durango state for 2019) [85]. Fuels in Site 1 ranged from semi-arid grass and shrublands in the eastern part of the fire to a more subhumid to temperate conifer and broadleaf understory and litter in the west (Figure 1a). Site 2 was dominated by pine and oak timber litter (Figure 1b), affected by a low-to-medium-severity fire in a humid climate in the west of the SMO. Finally, Site 3 was a medium-to-high-severity fire (with frequent crown fire patches) in a subhumid and temperate–cold pine timber understory mixture with shrublands in NW Durango (Figure 1c) [84].

#### *2.2. Field Sampling of Fire Severity*

Fire severity field sampling was performed during the months of November and December 2019 after fires at the three locations. Field and satellite imagery dates were selected to represent the early conditions of fire severity immediately after the fire (i.e., the initial assessment) in order to minimize the fire obscuration effects caused by rapid regrowth that have been documented for similar semi-arid ecosystems (e.g., [61,71]). We measured a total of 115 circular plots (30 m in diameter), distributed in a stratified sampling scheme, measuring homogeneously burned patches—stratified by *dNBR* [72]—within the fire perimeter (85 plots), as well as unburned patches (30 control plots). At each fire location, we also stratified sampled plots by vegetation types (Figure 1) and pre-fire *NBR* (e.g., [86,87]), both for burned and unburned plots. The center of each plot was georeferenced using a high-precision GPS (Spectra Precision *MobileMapper*© 60, with an average precision of <2 m) [88].

The field protocol for fire severity was based on the methodology for temperate forest described by Silva-Cardoza et al. [86,87]. This field protocol is an adaptation of different published guides for fire severity and mapping evaluation, i.e., Key and Benson [57], Parson et al. [89], Vega et al. [90,91], and Rodriguez et al. [92]. Both vegetation and soil strata were measured (Table 2). For every tree with a diameter at breast height (DBH) > 7.5 cm, we measured the mean scorch height, total height, crown base height, crown diameter, and DBH. The scorched crown volume (%), defined as the sum of the estimated percent assessed for the scorched and consumed crown of each tree, was visually estimated [93]. At the time of field data collection, the overstory stratum was defined as every individual with DBH > 7.5 cm [94]. Once we processed field data, the overstory was separated by 75th percentile heights values on each plot sampled into canopy (>75th percentile) and subcanopy (<75th percentile) [95]. Understory strata were defined as small trees and shrubs with DBH < 7.5 cm (Table 2), given the availability of allometric equations which can be used to estimate biomass for the main species in the study area [96]. For these strata, we deliberately did not include severity measurements in herbaceous plants, given their ephemeral nature and fast recovery in the analyzed ecosystems, among other challenges documented in the literature [97], and focused our evaluation on the effects on small trees and shrubs. For those, we quantified frequencies by basal diameter classes 0–2.5, 2.5–5, and 5–7.5 cm and measured the average total height, crown base height, crown diameter, and estimated scorched crown volume for a representative sample for each basal diameter class.

**Table 2.** Strata levels used in the field inventory performed in this study.


DBH = diameter at breast height, h = total tree height (m) and P75 = 75th percentile height per individual for each plot.

The assessment of fire impacts on soil—hereafter named soil burn severity (SBS) following Parson et al. [89], Keely [98], Vega et al. [91], and Sobrino et al. [8]—was performed on 9 square quadrats (30 × 30 cm) per plot, at every 5 m in 3 transects of 15 m (at 0, 120, and 240◦ azimuth) from the plot center. In each quadrat, the cover, depth, and intact/charred state of the remaining litter and duff were measured. The color and depth of ash layer, the degree of surface fine root consumption, and the soil structure alteration were also evaluated at each quadrat. To classify SBS, the following levels were considered, following Vega et al. [91] and Sobrino et al. [8]: 0: unburned; 1: burned litter but limited duff consumption; 2: litter and duff layers totally charred covering unaltered mineral soil; 3: litter and duff completely consumed (bare soil), but mineral soil organic matter not consumed and structure unaffected; 4: bare soil, surface mineral soil organic matter consumed and soil structure affected; 5: bare soil and mineral soil structure and color largely altered. For each quadrat, the percentage of cover by SBS level (0–5) was registered. The weighted SBS average by cover was calculated for every quadrat, and the plot-level SBS was finally obtained as the average of all the measured quadrats per plot. The SBS levels (1–5) were scaled from 0 to 100, named *soil burn severity index* (*SSI*), by multiplying by 20 to homogenize the levels with the vegetation strata severity scale.

Based on the analyzed strata, we calculated the following *field-based severity indices* (FBSs):

$$FSI = \frac{\sum(\text{CCS} + \text{SCS} + \text{UCS} + \text{SSI})}{4} \tag{1}$$

where: *FSI* = Fire Severity Index, *CCS* = canopy crown scorch volume (%), *SCS* = subcanopy crown scorch volume (%), *UCS* = understory crown scorch volume (%), and *SSI* = soil burn severity index.

The structure of the *FSI* index (0–100) is relatively similar to the first proposed Composite Burn Index (*CBI*) [57]), which averages the severity observed in all evaluated strata. In order to account for possible correlations in the crown damage between canopy and subcanopy strata, as well as potential correlations between understory severity and SBS (e.g., [90]), we tested a simplified *FSI2* index, using overstory strata (all trees with DBH > 7.5) and grouping the understory and soil fire severity as a surface stratum, following Equation (2):

$$FSI\_2 = \frac{OCS + \frac{UCS + SSI}{2}}{2} \tag{2}$$

where: *FSI2* = Fire Severity Index 2, *OSC* = overstory crown scorch volume (%), *UCS* = understory crown scorch volume (%), and *SSI* = soil burn severity index.

We additionally evaluated a field severity index weighted by the cover of each stratum, named *weighted FSI*, resembling the *Geometrically structured CBI* (*GeoCBI*) proposed by De Santis and Chuvieco [99], by multiplying each *FSI* stratum by its corresponding cover:

$$wFSI = \text{CCS} \ast \text{CC} + \text{SCS} \ast \text{SC} + \text{LCS} \ast \text{LIC} + \text{SSI} \ast \text{SOC} \tag{3}$$

where: *wFSI* = weighted Fire Severity Index, *CSC* = canopy crown scorch volume (%), *CC* = canopy cover (%), *SCS* = subcanopy crown scorch volume (%), *SC* = subcanopy cover (%), *UCS* = understory crown scorch volume (%), *UC* = understory cover (%), *SSI* = soil severity index, and *SOC* = soil cover (%), calculated as 100 − (*CCC* + *SCC* + *UCC*).

Finally, we also tested a simplified version of *wFSI*, named *wFSI2*, with a similar structure to *FSI2*, considering overstory and surface strata (soil and understory), weighted by their corresponding cover:

$$wFSI\_2 = OSC\*OC + \left(\frac{USC + SSI}{2}\right) \* \left(100 - OC\right) \tag{4}$$

where: *wFSI2*= weighted Fire Severity Index 2, *OSC* = overstory crown scorch volume (%), *OC* = overstory cover (%), *UCS* = understory crown scorch volume (%), and *SSI* = soil severity index.

#### *2.3. Remotely Sensed Data*

2.3.1. Aggregated Active Fire Perimeters

In order to locate the areas of fire occurrence to retrieve Sentinel images, we used monthly aggregated active fire (AGAF) perimeters [75], downloaded from the Fire Danger Forecast System of Mexico (Sistema de Predicción de Peligro de Incendios Forestales de Mexico, SPPIF) (http://forestales.ujed.mx/incendios2/, accessed on 16 May 2022), [100,101]. AGAF perimeters are calculated in the SPPIF by applying a convex hull algorithm to active fire data from the Moderate Resolution Imaging Radiometer Spectrum (MODIS, 1 km) [102] and the Visible Infrared Imaging Radiometer Suite (VIIRS, V1375 m) [103], as described in Briones-Herrera et al. [75,104]. Historical monthly and 10-day AGAF perimeters for Mexico are available for download from the SPPIF for the period 2012–current day [101]. They are updated in the SPPIF at every MODIS and VIIRS satellite pass [104], provided in near real time by the *Comisión Nacional para el Conocimiento y Uso de la Biodiversidad* (CONABIO) [105].

#### 2.3.2. Sentinel-2 Data

Data download and analysis for fire severity and burned area mapping involved three stages, as summarized in Figure 2: (i) *Section I* aimed to analyze the correspondence of the field-based severity indices (*FBSs*) and field severity data by stratum, with single-date paired images. The goal of this first analysis was to select the spectral index (*SI*) with the highest correspondence against *FBS* data for further evaluation in the subsequent stages. (ii) *Section II* aimed to evaluate the effect of the image compositing period (1 or 3 months), technique (average or minimum), and phenological correction (constant or relative, *c or rc*) for predicting field fire severity. The evaluated composites were downloaded from *Google Earth Engine* (*GEE*) using the selected *SI* from section I. (iii) *Section III* aimed to evaluate the *GEE* composites' accuracy in burned area perimeter mapping, by means of a kappa and overall accuracy analysis of burned/unburned points for each individual fire (Figure 2).

**Figure 2.** Summary of the study methodology for analyzing Sentinel-2 imagery performance for: (i) evaluating spectral indices from paired dates to predict field-observed fire severity (*Section I*); (ii) comparing image composites of varying length (1 or 3 months), composite technique (average or minimum), and phenological correction (*c* or *rc*) for predicting field fire severity (*Section II*) and for (iii) mapping burned area perimeter (*Section III*). Where: AGAF: aggregated active fire perimeters (2.3.1), *TOA*: top-of-atmosphere, *dSI*: differenced spectral index (Equation (29)), *RSI*: relative spectral index (Equation (31)), *c:* constant phenological correction (Equation (30)), *SI*: spectral index, *FBSs*: field-based fire severity indices, *FSI:* field severity index (Equation (1)), *FSI2*: field severity index 2 (Equation (2)), *wFSI*: weighted field severity index (Equation (3)), *wFSI*2: weighted field severity index 2 (Equation (4)), *AA*, *AM*, *MM*: compositing technique (Table 4), *rc:* relative phenological correction (Equation (32)), and *GCI:* Google Earth Engine composite index.

2.3.3. Section I: Comparison of Spectral Indices from Paired Images against Fire Severity

For the analysis of paired images, we downloaded single-date pre- and post-fire L1C images from the United States Geological Survey (USGS) Earth Explorer data portal (https://earthexplorer.usgs.gov/, accessed on 12 February 2022). We selected the images with the lowest cloud cover, closest to the dates before and after fire occurrence based on *MODIS* and *VIIRS* active fires dates for each fire (Table 3). The digital values of L1C-level images were transformed into apparent reflectance using top of atmosphere (TOA) correction, based on the Dark Object Subtraction (DOS) method, proposed by Chávez [106]. Earth Explorer images were processed in the Semi-Automatic classification plug-in proposed by Congedo [107] on QGIS software v3.16.16-Hannover [108].

Spectral Indices Analyzed from Paired Sentinel-2 Imagery

The Sentinel-2 spectral indices analyzed for fire severity discrimination are listed in Table 3. For all of the tested spectral indices *SI* (Table 3), we calculated the difference in spectral values between pre- and post-fire dates, *dSI*, following Equation (5):

$$dSI = \left(SI\_{pre} - SI\_{post}\right) \times 10^3 \tag{5}$$

where *dSI* = differenced spectral index [109] from paired pre- and post-*SI* images (*SIpre* and *SIpost*, respectively).

In order to account for phenological and weather-related influences between the preand post-fire images (e.g., [41,70,72,73]), we calculated the phenologically corrected spectral indices, *dSIc*, for all the considered indices, following Equation (6). For the constant "offset" *c* [73] calculation, we used a buffer from 3 to 5 km from the aggregated active fire perimeter, following Briones-Herrera et al. [75].

$$dSI\_{\mathcal{E}} = dSI - \mathcal{c} \tag{6}$$

where *dSIc* = constant phenology − corrected dSI and *c* = constant phenological correction, calculated as the average *dSI* value in the unburned vegetation.


**Table 3.** Spectral indices derived from Sentinel L1C and L2A imagery calculated in this study.

<sup>1</sup> Sentinel-2 MSI spectral bands: B1 = ~443 nm, B4= ~665 nm, *red-edge bands (\*re):* B5=~704 nm, B6 = ~740 nm, B7 = ~783 nm, *NIR:* B8 = ~833 nm, *NIR narrow (\*n):* B8A = ~865 nm, *SWIR:* B11= ~1613 nm, B12= ~2202 nm. Spatial resolution bands = 20 m (B1, B4, and B8 were resampled) [116]. <sup>2</sup> References with "-" are index modifications tested in this study.

In addition, for the *NBRi* indices considered (Equations (12)–(16), Table 3), we calculated the *Relativized Burn Ratio* (RBR) following Equation (31) [70]:

$$RSI\_{ic} = dNBR\_{ic} \;/\left(NBR\_{iprc} + 1.001\right) \tag{31}$$

where *RSIic* = *RBR* with constant phenological correction, *dNBRic*: differenced *NBRi* (Table 3) with constant phenological correction (Equation (6)), *NBRi pre*: pre-fire *NBRi* (Table 3) and *i* = 1, 2, 3, 4*n*, and 5*n* (i.e., *NBR1*, *NBR2*, *NBR3*, *NBR4n*, and *NBR5n*, Table 3).

#### 2.3.4. Section II: Comparison of GEE Composites for Mapping Fire Severity

Sentinel-2 composites were generated and downloaded from Google Earth Engine (GEE), a multi-petabyte catalog of satellite imagery and geospatial datasets which is a cloudbased platform [117]. GEE composite indices were derived from Sentinel-2 processing level 2A: bottom-of-atmosphere [GEE dataset ID: ee.ImageCollection ("COPERNICUS/S2\_SR")] imagery. SCL and Q60 bands were used to select the "optima" pixels (i.e., free of clouds, shadows, water, and/or snow) within a given period in each pre- and post-fire window.

#### GEE Composite Periods and Techniques Analyzed

Composite indices were generated with varying composite lengths of 1 or 3 months for the best *dSI* and *RSI* selected in *Section I*. Dates immediately before and after each fire, based on MODIS and VIIRS active fire data, were used to define each composite search period (Table S2). We tested average and minimum compositing since they are the most widely used approaches in the literature (e.g., [19,27,28,41,42,44,45,48]). For each 1- and 3- month composite, three combinations of *average* (*A*) and *minimum* (*M*) for the pre- and post- fire images (*AA*, *AM*, and *MM*) were tested (see Table 4). GEE code for calculating and downloading the composites was based on the modification of Parks et al. [41] for Sentinel-2 by Briones-Herrera et al. [75], initialized with active fire perimeters.

**Table 4.** Composite techniques analyzed for GEE Sentinel-L2A imagery.


Constant and Relative Phenological Correction

For each of the 6 combinations of composite length (1 or 3 months) and technique (*AA*, *AM*, or *MM*, Table 4), two phenological corrections were compared, resulting in a total of 12 composites calculated for each *dSI* and *RSI* analyzed. The first correction method tested was the constant phenological correction *c* [73]. The phenological corrections were calculated for each of the composites analyzed by subtracting the *c* offset, following Equations (30) and (31). The constant offset *c* was calculated in *GEE* following Parks et al. [41] and using the 3 to 5 km buffer of the aggregated active fire perimeter, as proposed by Briones Herrera et al. [75], to calculate the offset as the average *dSI* from the unburned buffer.

The constant phenological correction method applied the same correction *c* to every pixel, therefore assuming that phenological changes in *dSI* for the unburned vegetation were homogeneous. Additionally, in order to consider the potential spatial heterogeneity in unburned *dSI* changes, we calculated the average *dSI* (e.g., *dNBR*) for every pre-*SI* (e.g., *NBRpre*) value in the unburned 3–5 km buffer. This allowed us to assess variations in vegetation phenological changes (i.e., unburned *dNBR*) between vegetation with varying characteristics, since *NBRpre* values are frequently used as proxies of total biomass or tree cover [70]. To calculate a spatially variable phenological correction, we assigned the average unburned *dSI* (e.g., *dNBR*) to each pre-*SI* (e.g., *NBRpre*) pixel for the entire image extent, resulting in a map of relative phenological correction (*rc*). Such *rc* maps represented the expected phenological change (i.e., unburned *dNBR*) in the absence of fire for each type of vegetation, as stratified by *NBRpre* values. Finally, we calculated the relative corrected *rc* indices, *dSIrc*, using Equation (32):

$$dSI\_{rc} = dSI - rc \tag{32}$$

where *dSIrc* = relative corrected *dSI*, *dSI* = uncorrected differenced spectral index [109], and *rc* = relative phenological correction map.

Based on *dSIrc*, we also calculated its correspondent in the relative version *RSIrc* (Equation (33)) for all 6 of the composites analyzed.

$$RSI\_{\rm irre} = dNBR\_{\rm irre} / \left(NBR\_{\rm iprc} + 1.001\right) \tag{33}$$

where *RSIirc*= *RBR* with relative phenological correction, *i* = 1, 2, 3, 4*n*, and 5*n* (i.e., *RBR1rc*, *RBR2rc*, *RBR3rc*, *RBR4nrc*, and *RBR5nrc*, *dNBRiRc*: differenced *NBRi* (Table 3), with relative phenological correction (Equation (32)), *NBRi pre* : pre-fire *NBRi* (Table 3) and *i* = 1, 2, 3, 4*n*, and 5*n* (i.e., *NBR1*, *NBR2*, *NBR3*, *NBR4n*, and *NBR5n*, Table 3).

This resulted in a total of 24 composites analyzed for the two phenological corrections tested: *c* and *rc*.

#### *2.4. Predicting Fire Severity from Spectral Indices and Composites*

We evaluated linear and non-linear (exponential and power) regression models to predict paired-image spectral indices (*section I*) from field-measured severity variables by stratum (*OSC*, *CSC*, *SSC*, *USC*, and *SSI*) and from field severity indices, i.e., *FSI*, *FSI2*, *wFSI*, and *wFSI2* (see Equations (1)–(4)). For example, for a linear model, we tested the following regression (Equation (34)):

$$y = \beta\_0 + \beta\_1 \times \text{FBS} \tag{34}$$

where *y* = Sentinel-2 spectral index (*SI*), *FBS* = field-based severity per strata, i.e., *OSC*, *CSC*, *SSC*, *USC*, and *SSI*, and field severity indices, i.e., *FSI*, *FSI*2, *wFSI*, and *wFSI*<sup>2</sup> (Equations (1)–(4)).

Models were fitted using the "*lm*" and "*nls*" packages from the R software [118]. Model goodness of fit was evaluated based on the squared correlation coefficient between the predicted and observed values (R2) and their *Root Mean Square Error* (*θ*, RMSE) [119]. The best-performing spectral index was selected for further analysis in *Section II*, again using linear and non-linear (exponential and power) regression models, in this case to predict GEE composite indices of the selected absolute and relative spectral index (*dSI* and *RSI*) from the analyzed field severity indices. Composites' correspondence with field fire severity was also evaluated using model R2 and RMSE.

#### *2.5. Fire Severity Thresholds*

We defined the following fire severity categories based on field fire severity indices: Extreme (>90), High (60–90), Moderate (30–60), Low (30–10), and Very low/unburned (<10). The proposed fire severity categories considered both CBI protocol categories [72] and previous research in post-fire pine mortality in Mexico [92,120] and in temperate pine ecosystems elsewhere [10]. The threshold for High fire severity was based on such previous studies, which observed higher post-fire pine mortality when approximately >2/3 of the tree crown volume was scorched. Conversely, lower fire-caused tree mortality is expected at the Moderate fire severity level, which can be in turn more susceptible to pest attack, in contrast with low mortality either by fire or pests expected at the *Low* fire severity level [10]. Finally, to support post-fire emergency soil erosion hazard assessments, we differentiated an Extreme fire severity level, corresponding to areas of full tree canopy consumption, due to crown fire. In these areas, the absence of scorched litter covering the bare mineral soil can result in higher soil erosion risk [7,8,91]. This is particularly marked when crown fire is accompanied by high soil organic matter consumption, consistent with high soil burn

severity levels [9,90,121]. Based on these field fire severity categories, the best fit models of the best-performing spectral index and compositing technique were used to obtain the corresponding spectral index thresholds to map fire severity levels.

#### *2.6. Section III: Comparison of GEE Composites for Mapping Burned Area Perimeter*

In order to evaluate the burned/unburned (BUR/UNB) separability of the 24 analyzed composites for each wildfire (i.e., Sites 1–3), the *SI* values of the GEE composites were extracted to a sample of BUR/UNB points. We used filtered VIIRS active fires with *SI* > 0 inside the active fire perimeter of each wildfire as burned points (BUR). The same number of data sample points were randomly located inside the unburned buffer (UNB). Agricultural areas, agricultural–forest interface areas (downloaded from the SPPIF, accessed on 17 February 2021 http://forestales.ujed.mx/incendios2/#), bare soil, human settlements, water bodies, and ice/snow [122] were not considered for analysis. This resulted in a total of 5144, 634, and 1708 burned points for sites 1, 2, and 3, respectively, and an equal-sized sample of unburned points. For each wildfire, we tested the 99th, 95th, 90th, 85th, 80th, and 75th percentiles of *SI* values as candidates for the burned area threshold for the BUR sample, as well as the 1st, 5th, 10th, 15th, 20th, and 25th *SI* percentiles for the UNB sample. For every fire and for each candidate sampled composite *SI* threshold value, we calculated the overall accuracy, kappa coefficient, sensitivity (i.e., omission errors), and specificity (i.e., commission errors) [123] to classify burned/unburned data, using the '*caret*' package [124] (https://cran.r-project.org/web/packages/caret/, last accessed on 23 March 2022) from the *R* software [118]. The spectral index value with the highest kappa and overall accuracy was set as the optimal threshold to discriminate unburned from burned areas [23] for each composite and wildfire in question.

#### **3. Results**

#### *3.1. Section I: Correspondence between Field-Based Severity Variables and Spectral Indices from Earth Explorer (Paired Images)*

For the field severity variables and spectral indices analyzed, higher R<sup>2</sup> values were obtained with linear models, which are shown in Table S3. The highest correspondence with spectral indices was observed for *FSI* and *FSI2*, which performed very similarly, followed by *wFSI2*, generally outperforming *wFSI* (Table S3). *RBR5nc* and *RBR3c*, using bands 8a/12 and 8/12, respectively (Equations (14) and (16) [31]), had the first- and second-highest correspondence with the analyzed field measures of fire severity (an R2 of approximatedly 0.85 for both *FSI* and *FSI2*) and by individual stratum as well (Table S3). The third ranking spectral index was *dBAIS24nc* (Equation (10) [30]), with an R2~0.83 for both *FSI* and *FSI2* (Equations (1) and (2), respectively). *RBR5nc*, *RBR3c*, and *dBAIS24nc* also performed consistently, with the highest correspondence with overstory, canopy, and subcanopy severity and very similar performance between them for these strata (R2 = 0.80–0.78) (Table S3). On the other hand, lower but significant correlations were observed for the understory and soil strata, with an R<sup>2</sup> of 0.69–0.67 (Table S3). For the understory stratum, *RBR5nc* and *RBR3c* were also the best-performing *SI*, with an R2 of approximately 0.69, followed by *RBR1c* (using bands 11/12 [26], Equation (12), Table 4) and *dBAIS24nc*, with an R<sup>2</sup> of approximately 0.66 for both indices (Table S3). For soil burn severity, the highest correspondence was observed for *RBR5nc* (R<sup>2</sup> of 0.685), closely followed by *RBR3c* (R<sup>2</sup> = 0.678) (Table S3). The absolute fire severity spectral indices based on *NBR*, such as *dNBR5n* (8a/12 bands), consistently ranked higher than red and *red-edge* versions of *dNDVIc* for all strata but showed a lower performance compared to their relative versions [70] for all of the analyzed fire severity indices and strata.

#### *3.2. Section II: Comparison of Google Earth Engine Composites with Varying Period, Technique, and Phenological Correction for Predicting Field-Observed Fire Severity*

#### 3.2.1. Correspondence between GEE Composites and Field-Based Severity Indices

Based on *Section I*, we selected the composite with the highest correspondence index, *RBR5nc* (hereafter, *RBR*) for further testing in GEE. We also selected for further testing its corresponding absolute version, *dNBR5nc* (hereafter, *dNBR*), in order to quantify the differences between absolute and relative indices with the same spectral bands and because of the wide utilization of this latter index [67]. Similar to *Section I*, higher R<sup>2</sup> values were again observed with linear models. The correspondences between field fire severity indices and the composites for the selected indices *RBR* and *dNBR* with varying composite lengths (1 or 3 months), techniques (*AA*, *AM*, or *MM*), and phenological correction (*c* and *rc*) are shown in Table 5. In agreement with the results in *Section I*, *RBR* composites showed better performance than *dNBR.* The greatest association levels were obtained for the relative phenological correction indices (i.e., *RBRrc* and *dNBRrc*) compared with their analogous constant phenological correction (\**c*) (Table 5). These gains using *rc* composites were higher for the 3-month composites, where the R<sup>2</sup> values increased by 0.05–0.09 and 0.10–0.14 for *RBR* and *dNBR*, respectively (Table 5).

**Table 5.** Goodness of fit and coefficients of linear models between spectral indices from composite images and the field-based severity index *FSI2* for the three wildfires of Durango.


*β*0, *β*<sup>1</sup> = intercept and coefficient of linear regression model (Equation (34)), R2 = coefficient of determination, *θ* = root mean square error, *FSI2* = Fire Severity Index 2 (Equation (2)), Suffix: c = phenological correction [73], rc = relative phenological correction (suggested in this study), composite technique: *AA* = average in preand post-fire composites; *AM* = average (pre) and minimum (post); *MM* = minimum (pre and post), time window: 1 and 3 months, Diff. R2 (rc-c): observed difference in R2 between rc and c composites for corresponding technique and index, and Diff. R2 (1-3): observed difference in R2 between 1-month and 3-month composites for corresponding technique and index. Changes in R2 > 0.05 are marked in bold.

For all the analyzed compositing techniques and indices, the R<sup>2</sup> between field severity and composites was higher for the 1-month period than for the corresponding 3-month period (Table 5). The decrease in the R<sup>2</sup> with increasing composite length from 1 to 3 months was ≥0.05 for all of the *c* composites (0.05, 0.12, and 0.10 and 0.08, 0.15, and 0.14 decreases in the R2 for *AA*, *AM*, and *MM* for *RBRc* and *dNBRc*, respectively). Interestingly, *rc* composites showed, in contrast, a ≤ 0.05 decrease in the R<sup>2</sup> (0.03–0.04) for all *RBRrc* and *dNBRrc* composites with increasing composite length from 1 to 3 months (Table 5). Regarding the compositing technique, the highest correspondence was obtained using the minimum in pre- and post-fire (*RBRrc\_MM1* with R2= 0.832 and *RBRc\_MM1* with R2= 0.820), closely followed by *AM* (R<sup>2</sup> of 0.830 and 0.817 for *RBR5nrc* and *RBR5nc*, respectively). Conversely, composites using the average in both pre- and post-fire images had the lowest R<sup>2</sup> for the 1-month composites, with an R<sup>2</sup> = 0.815 (*RBRrc*) and 0.779 (*RBRc*). *AA* composites also showed lower R2 values than *AM* or *MM* composites for the 1-month *dNBR* composites (Table 5).

#### 3.2.2. Fire Severity Thresholds for Each GEE Composite Technique

Based on the best-performance regression models obtained in Table 5, we determined composite thresholds categorizing *FSI2* (%) as follows: Extreme (>90), High (60–90), Moderate (30–60), Low (30–10), and Very low (<10) (Table S4). As expected, there were some variations in the thresholds between composites of different lengths, compositing techniques, and phenological corrections applied. Interestingly, the thresholds were higher for the 1-month composites compared to lower thresholds for the 3-month composites. For example, the *Extreme* fire severity thresholds for the *RBRrc* composites, using *AA* and *AM*, would be 478 and 506, respectively, for the 1-month composites, compared to values of 392 and 405 (approximately 100 *RBR* lower), respectively, for the 3-month composites. Similar magnitude differences were observed for all fire severity categories, for both *RBR* and *dNBR*, for all the compositing techniques applied.

#### 3.2.3. Evaluating Phenological Variations for each GEE Composite Period and Technique

The differences between the average and minimum for 1- and 3-month *NBR* preand post-fire composites are illustrated in Figure 3 for two selected wildfires: Site 1 (Figure 3a–h) and Site 2 (Figure 3i–p). For the post-fire composites, lower *NBR* values, more clearly capturing the high-severity areas, can be observed for the minimum composites for both wildfires and compositing periods (Figure 3f,h,n,p). In contrast, the average post-composites showed a more diluted response, less clearly capturing high-severity areas in both wildfires (Figure 3b,d,j,l). This is particularly noteworthy for the post-composites in Site 2 (Figure 3j,l), where the average composites very weakly captured the burned surface, achieving positive (green and/or wet vegetation) post-*NBR* values, in contrast to a more distinct observation of burned areas in the corresponding minimum composite (Figure 3n,p).

The process of phenological correction, by sampling the unburned *dNBR*, is illustrated in Figure 4 for the 3-month *MM* composites of the three analyzed wildfires, Site 1 (Figure 4a–d), Site 2 (Figure 4e–h), and Site 3 (Figure 4i–l). The left column (Figure 4a,e,i) shows the uncorrected *dNBR*, together with the active fire perimeter buffer, used to sample the phenological changes in the *NBR* in the unburned area adjacent to the fire. An inspection of the *dNBR* within the unburned buffer suggests spatial variation in the phenological *NBR* changes, particularly for Site 1 and 2 (Figure 4a,e, respectively). For example, higher *dNBR* values can be observed in the southwest of Site 1 fire's unburned buffer (Figure 4a). Similarly, higher *dNBR* values can also be seen in the south of Site 2 s unburned buffer (Figure 4e), in contrast with lower values in the northern part of the unburned area. A comparison against the corresponding *NBRpre* maps for these sites (Figure 3g,o, respectively) reveals that such areas with higher phenological change corresponded to the drying or senescence processes of areas with high *NBRpre* (i.e., dense, humid forests, shown in blue in Figure 3g,o), in contrast with lower phenological change (*dNBR*, Figure 4a,e)

observed in the low-biomass (low-*NBR*) areas which were already dry in the pre-fire image (Figure 3g,o).

**Figure 3.** Pre- and post-fire *NBR* GEE composites with varying techniques average (**a**–**d**) and (**i**–**l**) and minimum (**e**–**h**) and (**m**–**p**) and different time windows: 1 month (left columns) and 3 months (right columns)) for Sites 1 (**a**–**h**) and 2 (**i**–**p**).

**Figure 4.** Phenological corrections on *dNBR* images, where *dNBR\** = without any correction (**a**,**e**,**i**), corr\_rc = relative correction map (**b**,**f**,**j**), *dNBRc* = constant-offset-corrected *dNBR* (**c**,**g**,**k**), and *dNBRrc* = relative-phenology-corrected *dNBR* (**d**,**h**,**l**) for sites 1 (**a**–**d**), 2 (**e**–**h**) and 3 (**i**–**l**). Very low/unburned thresholds (shown in green) were obtained for each fire in Section III (Table 6).

**Table 6.** Overall accuracy, kappa index, sensitivity, and specificity of the best-performing *RBR* burned area thresholds, for each analyzed composite technique in 1–3 months and different phenological corrections (*c* and *rc*).


Composite technique: *AA* = average in pre- and post-fire composites, *AM* = average (pre) and minimum (post) and *MM* = minimum (pre and post), Perc. = best-performing *RBR* percentile, \* percentile of the unburned sample, OA = overall accuracy, Sens. = sensitivity, Spec. = specificity, time window: 1 and 3 months, phenological correction: *RBRc* = *RBR* with constant phenological correction [73], *RBRrc* = *RBR* with relative phenological correction (proposed in this study), Diff OA (rc-c) = difference in overall accuracy between *rc* and *c* composites, Diff kappa (rc-c) = difference in kappa coefficient between *rc* and *c* composites, Diff. Sens. (rc-c): difference in sensitivity between *rc* and *c* composites, and Diff. Spec. (rc-c): difference in specificity between *rc* and *c* composites. Differences in overall accuracy, kappa coefficient, sensitivity, or specificity between *rc* and *c* composites higher than 0.05 are shown in bold.

The average unburned *dNBR* values for every *NBRpre* value, used to calculate the relative correction, are shown in Figure 5 for the three analyzed wildfires, for every composite technique and period tested. Changes in the *NBR* in the unburned areas were substantially larger in the 3-month composites (Figure 5b,d,f), reaching values of up to 500–600 average *dNBR* in the unburned areas. In contrast, significantly lower unburned *dNBR*, of up to 100–200 values, were observed in the 1-month composites (Figure 5a,c,e), for all the wildfires and composite techniques analyzed. Both composites using the minimum on the post-*NBR* (*AM* and *MM*, showed in blue and green in Figure 5, respectively) behaved similarly, suggesting a higher influence of the post-fire image than the pre-fire image on the observed phenological change for the studied wildfires.

**Figure 5.** Average *dNBR* by *NBRpre* in the unburned buffer for sites 1 (**a**,**b**), 2 (**c**,**d**), and 3 (**e**,**f**) for 1-month and 3-month composites (**a**,**c**,**e** and **b**,**d**,**f**, respectively), where composite technique: *AA* = average in pre- and post-fire composites, *AM* = average (pre) and minimum (post), and *MM* = minimum (pre and post).

Interestingly, for both 1 and 3 months, unburned *dNBR* was mainly positive (i.e., drying) for *AM* and *MM* composites in the majority of fires and periods analyzed (shown in blue and green in Figure 5). In contrast, for composites using the average in both preand post-fire images (*AA*, shown in red in Figure 5), the unburned *dNBR* varied from negative (i.e., greening or increased moisture) in low-biomass (low-*NBR*) fuels to positive (i.e., drying) in higher-biomass (high-*NBR*) fuels. For those *AA* composites, averaging

positive and negative *dNBR* values in a constant offset calculation could result in the underestimation of larger positive or negative phenological changes that occurred for each fuel type, not reflecting the spatial heterogeneity of phenological changes (Figure 5). In the case of composites with the minimum in the post-image *AM* and *AM*, a constant average offset could largely underestimate the higher drying or senescence present in unburned high-biomass fuels, particularly for the longer periods (Figure 5b,d,f). In addition, the smaller phenological changes that occurred in the low-biomass fuels (Figure 5) could also be overestimated by an average constant offset in this case.

By assigning these observed unburned average *dNBR* values (Figure 5) to each *NBRpre* value, we created a spatially variable relative correction (*rc*), shown in Figure 4b,f,j. The relative correction *rc* maps for the 3-month *MM* composites shown in Figure 4 use the patterns analyzed in the unburned buffer (Figure 5) to spatially represent the expected phenological change in the absence of fire for the entire area of study. Consequently, for Site 1 and Site 2 (Figure 4b,f), relative correction *rc* maps reflect that, on average, for those sites, pixels with a higher *NBRpre* experienced a higher phenological increase in the *dNBR* in the period analyzed. Conversely, observed lower phenological changes are mapped for the low-biomass (low-*NBR*) areas in those sites and periods of study (Figure 4b,f). Relative-phenology-corrected *dNBRrc* maps, calculated by subtracting *rc* rasters (second column, Figure 4b,f,j) from uncorrected *dNBR*s (first column, Figure 4a,e,i) are shown in the last column of Figure 4d,h,l. In contrast to the constant phenological *dNBRc* maps (Figure 4c,g,k), *dNBRrc* composites seem to minimize confusion (i.e., commission errors) in the unburned buffer (as evaluated in *Section III*), while also maintaining a clear distinction of the patches of higher severity (Figure 4d,h,l).

#### *3.3. Section III: Evaluating the Effect of Composite Period Length, Technique, and Phenological Correction on the Accuraccy in Mapping the Burned Area Perimeter for each Wildfire* 3.3.1. Burned Area Thresholds for Each Composite and Wildfire

We selected the best-performing *RBR* value (highest kappa and overall accuracy) from the candidate burned and unburned percentiles tested as the best burned area perimeter threshold for each wildfire and composite technique, period, and phenological correction (i.e., *c* and *rc*) (Table 6, Figure S1). Both the overall accuracy and kappa peaked at the same *RBR* values, which were in turn selected as burned area thresholds (Table 6). Burned area perimeters estimated with 1-month composites consistently showed higher kappa and overall accuracy values than their corresponding 3-month composites; this decrease in kappa with increasing composite period was higher for *c* than for *rc* composites (Table 6, Figure S1). Regarding the compositing technique, composites using the minimum (*AM* and *MM*) largely outperformed those using the average (*AA*) in Site 2 (Table 6 and Figure S1). For both Site 3 and 1-month composites in Site 1, relatively similar results were obtained for both the average and minimum techniques (Table 6 and Figure S1).

Remarkably, *rc* composites had generally higher kappa and overall accuracy values than their corresponding *c* composites (Table 6 and Figure S1). The highest increases in kappa were observed for Site 2 (where kappa values raised by 0.16–0.22 and 0.06–0.13 for 3 months and 1 month, respectively) and Site 1 (with kappa increases of up to 0.04–0.07 for the 3-month composites). For Site 2, those improvements in accuracy were generally a result of both gains in sensitivity of 0.09–0.16 (i.e., lower omission errors) and increases of up to 0.11 in specificity (lower commission errors) for the *rc* composites (Table 6). For the 3-month composites in Site 1, accuracy increases were mostly caused by improved specificity values of up to 0.10 (i.e., decreased commission errors). In contrast, for Site 3, both phenological corrections and all the composite techniques performed relatively similarly, together with a lower observed reduction in kappa values from 1 to 3 months.

The performance of *rc* and *c* composites for burned area mapping for the 1- and 3 month *MM* composites is shown in Figure 6, using the burned area thresholds obtained for each composite and wildfire analyzed (Table 6). Higher commission errors in the sampled unburned points (shown in yellow in Figure 6) can be observed for the *c* composites, particularly for the 3-month composites in Site 1 and 2 (Figure 6c,g), compared to a lower commission observed for their corresponding relative correction *rc* composites (Figure 6d,h). In addition, higher omission errors (shown in blue in Figure 6) for *c* composites in Site 2 can also be observed at both 1 (Figure 6e) and 3 months (Figure 6g), in contrast with lower omission errors for the *rc* composites during both periods (Figure 6f,h, respectively). Contrarily, slightly higher omission errors, mainly in low-severity areas, were observed with *rc* 3-month composites in Site 1 (Figure 6d), which nevertheless had a higher overall classification accuracy and kappa because of the large improvement in commission errors (Table 6).

**Figure 6.** Burned area (shown in red) estimated for site 1 (**a**–**d**), 2 (**e**–**h**) and 3 (**i**–**l**), using the corresponding burned area threshold from Table 6, for 1-month (left columns) and 3-month (right columns) *MM* composites using constant phenology correction *c* [73] and relative phenological correction *rc* (proposed in this study). Sampled points for burned areas (active fires) and for unburned areas are shown in gray and green within the active fire perimeter and the unburned buffer, respectively. Omission and commission errors in the sampled burned and unburned point sample are shown in blue and yellow, respectively.

#### 3.3.2. Fire Severity Mapping

We finally mapped fire severity (using *FSI2* thresholds from Table S4), also considering the wildfire-specific burned area thresholds for the Very low/unburned category (see Table 6). Fire severity maps for the three analyzed wildfires are shown in Figure 7 for the 1-month and 3-month *MM* composites, using constant and relative phenological corrections (*c* and *rc*). Notably, relative phenological corrections, in addition to allowing for a clearer burned/unburned delineation as previously shown in Figure 6 for these same composites, also maintained a clear distinction of high- and extreme-severity patches, even slightly improving the correspondence with field measured plots of fire severity (Table 5), which are shown in circles in Figure 7.

**Figure 7.** Fire severity mapping with *MM* technique and different time windows (1 and 3 months) and constant (*c*) or relative phenological correction (*rc*) for wildfires in Site 1 (**a**–**d**), Site 2 (**e**–**h**), and Site 3 (**i**–**l**), using *FSI2* thresholds (see Table S4). The corresponding burned area threshold (Very low/unburned, shown in green) was assigned for each wildfire (see Table 6) according to the time window and composite technique.

#### **4. Discussion**

#### *4.1. Field-Based Fire Severity Metrics and Paired Images S2 Spectral Indices Correspondence*

The observed higher performance of *RBR5nc* and *RBR3c* (using bands 8a/12 and 8/12, respectively) as the best predictors of field-based fire severity (Table S3), agree with previous Sentinel-based studies of fire severity (e.g., [9,34,125]) and burned area (e.g., [35,126]). The superior performance of the relative index *RBR* compared to its absolute version *dNBR*, both for paired images and composites evaluated (Table 5 and Table S3), also supports the findings of previous studies such as Fernández-Manso [29], Botella-Martínez [127], Arellano et al. [128], Quintano et al. [129], Parks et al. [41,70], Holsinger et al. [19], and Fernández et al. [9], where relative indices showed a higher correlation against field-based severity data. In our study, the strong heterogeneity in fuel types (Figure 1) and prefire biomass (Figure 3) was probably related to this better performance of the relative spectral indexes.

The higher performance of fire severity indices (i.e., *FSI*, *FSI2*) than individual strata (i.e., *OCS*, *UCS*, and *SSI*) supports previous observations where composite indices of fire/burn severity showed a higher R<sup>2</sup> with spectral indices than the individual strata evaluated [61,130]. For the temperate to semi-arid pine/oak study sites analyzed, our hypothesis that there could be correlations between canopy and subcanopy crown scorch volume and between understory and soil burn severity was met; in fact, we observed a significant correlation in both cases (R Pearsons = 0.96 and 0.77, respectively). This could explain our observed similar, even slightly higher performance of *FSI2* against *FSI* and, more interestingly, our markedly higher R2 for *wFSI2* against *wFSI* (Table S3). This suggests some potential of *wFSI2* to compensate for the overemphasis of soil strata in sparsely vegetated plots, observed using *wFSI*, with the latter being in line with previous observations [131]. Further research in a variety of forest ecosystems, analyzing the presence or absence of correlations between strata and also against satellite indices (e.g., [132]), might be useful in better understanding the limitations of both satellite and current integrated field indices of fire severity, and it could provide suggestions for potential improvements for simplified or more representative field indices.

Furthermore, while it is clear that integrated field fire severity indices can improve correlations with satellite data, the question remains as to whether such integrated indices are representative of ecosystem effects (e.g., [1,2,98,133]). In this sense, although they were lower than field fire severity indices, significant correlations were observed for all strata (Table S1). Such observed high correlations with overstory strata has been extensively documented elsewhere (e.g., [61,66,130]). This relationship can be potentially useful when coupled with mortality prediction models (e.g., [11,92,120]), while also accounting for forecasted post-fire pest risk (e.g., [10]), among other factors. Nevertheless, because of the complex interactions behind tree mortality, the performance of remotely sensed crown scorch volume as an input for predicting tree mortality is a subject of ongoing research (e.g., [12]).

Our observed relatively lower correlation for the understory and soil strata (R<sup>2</sup> of 0.691 and 0.685, respectively, Table S3) supports previous observations of lower accuracies for such strata against *SI* (e.g., [34,61,66]), with this limitation being potentially more marked in ecosystems with relatively high tree density. In this sense, for our study, the relatively low tree density in the majority of the pine/oak forests measured might have allowed for the signal of burned soil to be partially captured by the Sentinel spectral indices. This would agree with the observations of Allen and Sorbel [130], who also found significant correlations of substrate severity scores (R = 0.74), although they were also lower than for canopy (R = 0.83) in open-woodland canopies and tundra forest ecosystems. In addition, our observed outperformance of the *NBR* against *NDVI*-based indices for SBS prediction agrees with the results of Sobrino et al. [8] in temperate forests in NW Spain. Additionally, the best performance of *RBR* in the current study supports the best performance observed by Fernández et al. [9] for this index to predict post-fire physical soil properties.

While the relatively open canopies of our studied ecosystems seem to reinforce the concept that the ash and burned soil spectral response might be partially captured in relatively low-tree-density forests, as previously suggested by both some empirical [8,9,130] and radiative transfer modeling studies (e.g., [134]), this might not hold true under more dense canopy conditions (e.g., [66]). In this sense, several authors have stated that passive remote sensing imagery is more likely to detect fire effects on tree crowns than on the ground (e.g., [34,61,135]), together with its limitations to separate strata effects in ecosystems with complex structures [136]. This could be related to a relative inability of passive remote sensing to "pass through" over outer vegetation strata [34], this latter limitation being potentially more marked in ecosystems of relatively high tree density or cover. Accordingly, beyond using satellite reflectance alone, future studies could further explore the role of other ancillary variables such as land surface temperature from thermal bands, which has shown potential for predicting both soil burn severity [137,138] and field indices of fire severity such as CBI [139,140]. In addition, the role of topographic variables [141] or fire intensity and duration [142], with the latter possibly being inferred from active fire detections [143], could be also explored to improve our ability to predict and map fire effects on soil.

#### *4.2. Evaluation of GEE Composites for Mapping Fire Severity and Perimeter*

#### 4.2.1. Effect of Compositing Technique on Fire Severity and Burned Area Mapping

Higher correspondence with field fire severity (Table 5) and increased burned area classification accuracy (Table 6) were observed for the composites using a minimum in the post-fire image, particularly for the *MM.* In addition, for the equivalent evaluation period and phenological correction, the fire severity thresholds for *MM* were higher than for their equivalent *AA* composites (Table S4). Both results are possibly related to their better potential to detect initial post-fire effects, by selecting the driest lower-biomass images both in pre- and post-fire periods (e.g., Figure 3h,p), and therefore being potentially less influenced by weather oscillations or regrowth effects in the composite period. In contrast, for average composites, the initial post-fire severity was more diluted, likely by forest recovery responses (e.g., Figure 3b,d) and possibly also by increased fuel and soil moisture with post-fire precipitation (Figure 3j,l), which was undesirable for burned area perimeter and initial fire severity detection. This supports studies that have proposed the use of the minimum or maximum values (depending on the spectral index) for burned area perimeter mapping (e.g., [27,51,144]). For example, Roteta et al. [26] reported that minimum *NBR* composites retained the pixel values when the burned signal was considered to be strongest, corroborated by our comparisons of post-fire minimum versus average *NBR* composites (Figure 3). Additionally, it confirms suggestions from previous studies that warned about the limitations of average composites for initial assessments in forests that revegetate rapidly after fire (e.g., [41,45]). In contrast, for a different goal, such as mapping extended assessments and long-term ecosystem responses, average composites of longer compositing periods might be valuable, particularly in forest ecosystems of slower postfire recovery (e.g., [19,41,45]). These potential differences support the observations of Holsinger et al. [19], who suggested that the best approach for measuring fire severity may depend on the context and conditions of each fire, as well as the measurement objectives. Future studies in Mexico should focus on performing extended assessments of burn severity, potentially capturing delayed mortality, which might be particularly relevant for firesensitive species such as wet tropical forests [84].

#### 4.2.2. Effect of Compositing Period on Fire Severity and Perimeter Mapping

Shorter-length (1-month) composites showed higher fire severity thresholds (Table S4) and higher correspondence with field fire severity (Table 5) than 3-month composites for the fires analyzed. Our observed decrease in the *dNBR* and *RBR* fire severity thresholds with increasing time confirms, for Sentinel-2 composites, previous observations from Landsat compositing studies (e.g., [19,41,45]) or from paired Landsat imagery (e.g., [54,58,61,130]), where longer evaluation periods also generally resulted in lower fire severity thresholds and lower correspondence with fire severity [19]. Furthermore, this better agreement of fire severity with imagery immediately after a fire seems to support the previous suggestions to, ideally, use an image captured within a few weeks of field plot measurements, so that the image reflects the conditions at the time of field measurements (e.g., [19,58,72]). In addition, the use of MODIS and VIIRS active fires to define the composite search window starting immediately post-fire (e.g., [19,47]) was confirmed to be useful, in the light of the fast decrease in the *dNBR* with increasing time since fire. In our study, an initial obscuration of fire effects, likely related to forest recovery, was already captured in the third month after fire, being particularly more marked when using average composites and constant offset corrections (Table S4, Figures 6 and 7). This likely agrees with similar rapid obscuration effects by regrowth in the first few months after a fire in fast-recovery ecosystems, such as semi-arid and tropical areas in the southwestern and southeastern United States and grasslands areas in the CONUS [61,71], or Mediterranean climate locations in southern Europe (e.g., [54]). We consequently recommend that, for the evaluation of initial fire severity in the semi-arid to temperate pine/oak ecosystems in Northern Mexico and similar fast-growing ecosystems, images taken immediately after fire and short compositing periods (<3 months) should be favored.

For our study area, images taken directly after a fire (1 month) were also particularly advantageous for fire perimeter mapping (Table 6, Figure 6). This supports various studies that have preferentially used early post-fire images for perimeter delineation, given the relatively clearer separation of burned/unburned areas in the fresh burn scar before ecosystem recovery (e.g., [54,62,68,145,146]). This need for short-term imagery has to be balanced with image availability, for which multi-sensor approaches [49,129,147] should be further explored [22]. Nevertheless, cloud-free, post-fire, phenologically matching images might not be readily available in the first month (e.g., [21,62]), particularly in in areas with persistent cloud cover. Consequently, in addition to defining the best season and time period for compositing, continued research in improved phenological correction techniques will likely continue to be necessary for sound burned area and fire severity evaluations.

Interestingly, the decrease in the R<sup>2</sup> for predicting field fire severity with increasing composite period was only higher than 0.05 for all the constant offset correction *c* composites. In contrast, *rc* composites using the proposed relative phenological correction showed two advantages for fire severity mapping: first, they provided a higher correspondence with field fire severity, particularly for the 3-month composites. Secondly, unlike *c* composites, *rc* composites experienced a much smaller decrease in the R<sup>2</sup> with increasing compositing periods to 3 months against corresponding 1-month *rc* composites. This suggests the better performance and potentially greater stability of the proposed relative phenological correction to minimize the observed decrease in correspondence against field fire severity with increasing observation periods. This should be corroborated in future studies testing this novel relative phenological correction technique in other locations.

#### 4.2.3. Effect of Phenological Corrections on Fire Severity and Burned Area Mapping

The observed values in the unburned *dNBR* were much larger than the range of −100 and 100 *dNBR* proposed by Key and Benson [57] or Picotte [74] and Picotte et al. [67] to be sampled in the unburned buffer for a constant offset calculation, particularly for the 3-month period (Figure 5). Such observed strong spatial variability in the unburned *dNBR* would violate the assumptions for a reliable application of the constant offset method [57,74]. This limitation is not uncommon: in their evaluation of the 18,497 images manually processed by MTBS analysts in the 1984–2014 period, Picotte [74] found that close to one third of the calculated offset values did not meet the recommended range for a reliable constant offset application (offset average and standard deviation exceeding +/− 50). This common unsuitability of available imagery, limiting reliable constant phenology correction because of frequent strong, spatially variable phenological effects, reinforces the need to account for the temporal and spatial variability in phenology changes. In contrast

to a constant offset calculation approach, where such strong variability in the range of the unburned *dNBR* can limit the average offset reliability, the relative correction method allows for the inclusion of the spatial variability in phenological changes in unburned areas, and therefore should not limit sampling values exceeding 100 or −100 *dNBR.* In fact, such variability can and should be incorporated into a spatially variable phenological correction to take into account the actual changes observed in different fuel types.

The relative phenological correction improved both the correspondence with fieldobserved fire severity (Figure 7) and burned area perimeter accuracy (Table 6 and Figure 6), mainly because of reduced commission, but also because of reduced omission errors in some sites. Improvements in specificity (reduced commission) in sites 1 and 2, which had large phenological effects in the unburned fuels, particularly higher during the longer compositing periods (Figure 5), were likely a consequence of the ability of *rc* to minimize the confusion of unburned stressed or senescent vegetation with burned areas by considering a fuel-specific (*NBRpre*) phenological change that better reproduces the spatial variability in those changes than a constant phenological correction.

Regarding omission errors, for Site 2, which had the strongest phenological effects, sensitivity was also largely improved using *rc*, possibly as a consequence of this clearer separability between burned and stressed vegetation for this site compared to the constant offset method. In contrast, the observed slight increase in omission errors for *rc* against *c* for the 3-month composites in Site 1 seems to suggest that, in some locations, in order to minimize commission errors, the *rc* perimeter thresholds might omit some of the lowestseverity areas. This is a common limitation using remote sensing imagery (e.g., [97]) and is relatively more acceptable than commission errors, since those low-severity areas are expected to be mildly impacted by fires and require fewer or no post-fire rehabilitation interventions (e.g., [10,90]).

In any case, omission errors should be treated with caution, given the spatial uncertainties inherent to the use of filtered active fire observations as proxies for burned area locations. This approach, already proposed by Roteta et al. [27], among others, for calibrating burned area thresholds for Sentinel and Landsat, has the advantage of reducing the need for extensive ground datasets of burned areas and of being susceptible to being automated in near real time using active fire observations to obtain fire-specific burned area thresholds [148]. Nonetheless, the relatively coarse resolution of active fires requires caution, particularly in smaller fires (with a lower number of observations than the large fires analyzed here), where the technique might have limited applicability because of this spatial uncertainty. Future studies could alternatively test the use of Landsat active fire detections [149] or aerial infrared imagery [150], where available, to calibrate locally adaptive burned area thresholds [148] from the percentile distribution of the spectral values of S2 indices, extracted from those detections at much finer spatial resolutions.

We observed a high degree of spatial and temporal variation in phenology, as measured by *rc*. In this sense, the most widely used method of phenological correction, the constant offset, had previously been documented to slightly improve correlations with CBI when analyzing several fires (e.g., [41,76]). Conversely, some studies reported decreased correlations with field CBI using constant offset correction methods compared to uncorrected images (e.g., [58]) and attributed this to the dynamic nature of the ecosystems studied, such that a single average *dNBR* value in unburned areas could not provide a reference point between pre- and post-fire images. This was also acknowledged by Parks et al. [41], who warned that, in their current formulations, constant offset methods would not be appropriate under heterogeneous fuel conditions, which can unfortunately be very common (e.g., [23,74]). Moreover, on an individual fire level, a constant offset correction is not expected to improve neither correlations with field fire severity (e.g., [130]) nor to change the individual fire perimeter delineation accuracy. It is obvious that simply subtracting a constant value to a raster image would not alter the relative differences between pixel values (e.g., burned/unburned) at an individual fire level. In contrast, spatially heterogeneous methods, such as the relative correction method, apply a different phenological correction

to every fuel type, as defined by the *NBRpre*, by including the observed heterogeneity in phenological changes (Figures 4 and 5), resulting in the potential for improving accuracies for fire severity and burned area mapping (Table 6, Figure S1).

Interestingly, the higher increases in accuracy using *rc* composites for fire severity and fire perimeter mapping, observed for the 3-month composites, corresponded with the higher variation in the *dNBR* values observed in the unburned area for those longer periods (Figure 5). In addition, larger improvements using *rc* against c composites for burned area perimeter mapping also corresponded to sites where observed phenological changes were also larger, such as Site 1 and 2, while accuracies and fire severity and perimeter maps (Figures 6 and 7) were similar for sites with smaller observed phenological changes, such as Site 3. This seems to suggest that relative phenological correction methods could be routinely applied both for fire severity and fire perimeter delineation, resulting in a similar performance to the current constant offset method when phenological changes are small. More importantly, when large and spatially heterogeneous phenological changes are present, such spatially variable phenological correction could improve both fire severity and perimeter delineation accuracy. Because of the high number of images with large and/or spatially variable phenological effects (31% of all MTBS images analyzed by [74]) and the frequent lack of immediate post-fire imagery [62], we believe that the potential of the proposed *rc* method to improve both fire severity and perimeter delineation can be significant.

Previous approaches have used supervised methods that require a manual selection of control plots by subjectively interpreted fuel type based on expert knowledge (e.g., [77]) or the manual creation of phenology-matched composites [22], which currently require the large involvement of analyst manual supervision [23]. Unlike those approaches demanding manual supervision, the novel relative correction method proposed in this study uses a systematic approach by automatically stratifying by *NBRpre* as a proxy of fuel type and biomass and quantifying and correcting for the phenological changes for each of those spatially variable fuel categories. All the processes in the *rc* calculation can be easily automated and readily incorporated into existing operational GEE codes, very similar to the current automatization of constant offset [41,75], but only by automatically stratifying its calculation by the *NBRpre*, with this layer also already being calculated automatically in Parks et al. [41] GEE code. Including an automated *rc* calculation into existing GEE burned area tools could improve the accuracy of operational large-scale automated assessments of fire severity [19,41,44,45] and burned area [26,27,42]).

Additionally, another advantage of the proposed *rc* method is that, unlike methods that assume a baseline from previous years (e.g., internannual differencing [23]), this approach neither requires a long time series as a baseline, nor does it make any assumptions about the magnitude of the phenological changes from those previous periods. Instead, it quantifies the observed temporal and spatial variations in the unburned fuels and uses this observed change to normalize the image with respect to observed unburned phenological changes. In this sense, this method can capture phenological changes of varying magnitudes (Figure 5), potentially reflecting a variety of greening/drying conditions between different fuels, specific to each location and assessment time. The portability of the *rc* method suggests the potential for improved analysis in a variety of locations and assessment periods, as it minimizes non-fire-induced changes in vegetation dynamics. This could improve the transferability of fire severity thresholds, again, because it minimizes fuel-specific variations in phenological and weather influences. Further research in spatial and temporal variations of phenological corrections across a large variety of ecosystems should be encouraged and ultimately incorporated into operational burned area and fire severity mapping efforts, potentially improving their accuracy and robustness.

#### **5. Conclusions**

For the first time, the current study analyzed the accuracy of Sentinel-2 spectral indices to map fire severity and burned area perimeter in Mexico. It included an evaluation of Google Earth Engine (GEE), composite length, and technique and tested a novel relative phenological correction (*rc*) against the constant offset (*c*) method. The *RBR* index, using bands 8a and 12, showed the highest correspondence with both field fire severity indexes and evaluated individual strata (i.e., canopy, subcanopy, understory, and soil), for a variety of sites with marked vegetation heterogeneity under a climatic gradient from semi-arid to temperate.

One-month GEE composites, using the minimum metric for post-fire images, showed a higher correspondence with field severity and higher fire perimeter delineation accuracies than 3-month composites or using the average. Interestingly, the decreases in fire severity and perimeter accuracy with increasing compositing period from 1 to 3 months were lower for the composites using the *rc* phenological correction, which showed greater stability compared to a constant offset correction. Therefore, the *rc* method could be potentially useful to minimize such effects when immediate post-fire, cloud-free imagery is not available.

Relative phenologically corrected composites showed both improved correspondence with field fire severity and higher accuracy for fire perimeter delineation for most sites and periods analyzed. This latter improvement was more significant in sites and periods where phenological effects (as quantified by changes in the *dNBR* in the unburned area) were stronger. Because strong and heterogeneous fuel phenology effects have been frequently observed, the *rc* method could have the strong potential to significantly improve the accuracy of fire severity and burned areas if routinely applied in post-fire evaluations.

The proposed *rc* approach, using the *NBRpre* to systematically stratify fuels and sample the spatially variable temporal variations in the spectral indices of the unburned fuels, goes beyond the constant offset method in incorporating fuel heterogeneity but shares its relative simplicity of being fully automated. Including such a spatially variable phenological correction into existing automated tools, such as *GEE*, could contribute to deriving improved, systematic large-scale evaluations of fire severity and perimeters in Mexico and elsewhere.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs14133122/s1, Figure S1: Kappa value across RBR values tested as candidate burned area thresholds for each wildfire (sites 1,2,3). Where: Phenological correction: c = constant offset phenological correction [73], rc = relative phenological correction (proposed in this study), Composite technique: AA = Average in pre and post fire composites (Figure 7a,b), AM = Average (pre) and Minimum (post) (Figure 7c,d), MM = Minimum (pre and post) (Figure 7e,f) and Time-window: 1 month (left column) and 3 months (right column).; Table S1: Acquired Sentinel 2 L1C level images from Earth Explorer summary; Table S2: Sentinel-L2A imagery period window for GEE composites indices calculation, Table S3: Coefficient of Determination (R2) and RMSE (θ) of univariate linear models between spectral indices (Table 3) (paired single day images) and field-based severity variables in the three wildfires of Durango, Mexico, Table S4: GEE Composites Index fire severity thresholds obtained using equations from table for FSI2 (%) levels where: Very Low/unburned (<10), Low (10–30), Moderate (30–60), High (60–90) and Extreme (>90).

**Author Contributions:** A.I.S.-C.: formal analysis, methodology; D.J.V.-N.: conceptualization, methodology, formal analysis, writing—original draft; J.B.-R.: programming; C.I.B.-H.: programming; P.M.L.-S.: methodology; J.J.C.-R.: methodology; S.A.P.: programming, writing—review and editing; L.M.H.: programming, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** Funding for this study was provided by CONAFOR/CONACYT Project "CO-2018-2-A3- S-131553, Reforzamiento al Sistema Nacional de Predicción de Peligro de Incendios Forestales de México para el pronóstico de conglomerados y área quemada (2019-2022)", for the enhancement of the Forest Fire Danger Prediction System of Mexico to map and forecast active fire perimeters and burned area, funded by the Sectorial Fund for forest research, development and technological innovation "Fondo Sectorial para la investigación, el desarrollo y la innovación tecnológica forestal". Additional support was provided by the USDA Forest Service, Rocky Mountain Research Station, and Aldo Leopold Wilderness Research Institute. The findings and conclusions in this publication are

those of the authors and should not be construed to represent any official USDA or U.S. Government determination or policy.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** MODIS and VIIRS active fire data used in the study can be accessed from FIRMS: https://firms.modaps.eosdis.nasa.gov/active\_fire/ (accessed on 26 March 2022). Aggregated active fire perimeters for the study period analyzed are publicly available through the Forest Fire Danger Prediction System of Mexico: http://forestales.ujed.mx/incendios2/ (accessed on 26 March 2022). Sentinel-2 single-date L1C images can be freely downloaded from the United States Geological Survey (USGS) Earth Explorer data portal (EE, https://earthexplorer.usgs.gov/, accessed on 12 February 2022). Monthly composite indices of Sentinel-2 processing level 2A: bottomof-atmosphere (GEE dataset ID: ee.ImageCollection ("COPERNICUS/S2\_SR")) imagery can be freely accessed from *Google Earth Engine*.

**Acknowledgments:** We would like to thank CONAFOR's personnel for their support during the current study and their continuous feedback regarding the development of the *Sistema Nacional de Predicción de Peligro contra Incendios Forestales* de México.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

