Global Analysis of Burned Area Persistence Time with MODIS Data

: Biomass burning causes a non-permanent land cover change (burned area) through the removal of vegetation, the deposition of charcoal and ashes, and the exposure of soil; the temporal persistence of these changes is highly variable, ranging from a few weeks in savannas to years in forests. Algorithms for the generation of moderate-resolution (10–30 m) continental and global burned area maps have been prototyped in an effort to meet the needs of diverse users of ﬁre information. Nevertheless, moderate-resolution sensors have reduced the temporal resolution (e.g., to 16 days for Landsat), which could potentially lead to omission errors, especially in ecosystems where the spectral signal associated with burning disappears quickly and cloud cover limits the number of valid observations. This study presents a global analysis of the burned area persistence time, deﬁned as the duration of the spectral separability of the burned/unburned areas mapped by the MODIS MCD64 Global Burned Area Product. The separability was computed by analyzing time series of normalized burn ratio (NBR) from nadir BRDF-adjusted MODIS reﬂectances (MCD43 product). Results showed that, globally, the median burned area persistence time was estimated at 29 days, and 86.6% of the global area, as detected by MODIS, can only be detected accurately for up to 48 days. Thus, results indicate that burned area persistence time can be a limiting factor for global burned area mapping from moderate-resolution satellite sensors, which have a low temporal resolution (e.g., Landsat 16 days, Sentinel-2A/B 5 days). ), interquantile range, minimum and maximum of the distribution.


Introduction
Fire is a natural component of any ecosystem, and it has effects on vegetation, soil, water and atmospheric composition [1][2][3]. Fire contributes to the global carbon cycle through the emission of greenhouse gases and aerosols from biomass burning [4,5], and as an agent of ecological change influencing the global vegetation dynamics [6][7][8] and the surface energy budget [9,10]. Earth Observation (EO) data allow the analysis of biodiversity, dynamics of biomass, productivity, and disturbances spatially from regional to global scales [11]. The application of satellite data to studying vegetation fires has greatly enhanced the possibility of introducing the effect of fire disturbances in global models of climate and atmospheric composition and dynamics estimation [12]. Several studies have used post-fire albedo observations to investigate fire-induced vegetation change impacts on radiative forcing and local effects on climate [9,13,14]. Due to their significant effects on terrestrial ecosystems and atmospheric processes, fire disturbances are included in the list of the required Essential Climate Variables (ECV) in support of the work of the United Nation Framework Convention on Climate Change (UNFCCC) and the Intergovernmental Panel on Climate Change (IPCC) [15].
Burned area maps are defined as the primary variable of the Fire ECV [15], and are currently one of the main inputs for the estimation of atmospheric emissions due to biomass burning [4,15,16].
Distribution Function (BRDF)-adjusted reflectance product (MCD43A4) time series of burned pixels, as detected by the MODIS global burned area product (MCD64A1), within the same ecoregion and land cover were used to characterize the pre-and post-fire temporal variations with the passing of time due to dissipation charcoal and ashes, vegetation regrowth, and snow cover. Reflectance time series of unburned pixels within the same ecoregion and land cover were used to characterize the variations in reflectance due to vegetation phenology and other disturbances. Summary metrics were defined for reporting the results at the global scale, aggregating by global biome and realm.
The paper is organized as follows. Section 2 describes the satellite datasets used for the analysis, and the ancillary datasets used as stratification variables. Section 3 describes the methods, establishing a rigorous probabilistic framework for the definition of the persistence time of burned areas at different scales, and presenting the formulae needed for a global, multiyear analysis. Section 4 presents the results and Section 5 discusses the differences observed across land cover types, biomes and realms. The paper concludes with recommendations for future research and application.

MODIS Global Burned Area Product
The most recent Collection 6 MODIS Global Burned Area Product (MCD64A1) provides the estimated date of burn for the 500 m MODIS pixels that are classified as burned to within a calendar month [18,19]. The global monthly MCD64A1 data record from 2003 to 2016 was used, in order to include every full year in which both MODIS instruments-on the Terra and Aqua platforms-were operating. The algorithm [18] was designed to be extremely tolerant of cloud and aerosol contamination, which affected the Collection 5 MODIS 500 m burned area product [49]. The algorithm applies dynamic thresholds to composite MODIS Terra and Aqua imagery generated from a burn-sensitive spectral band index derived from MODIS 1240 nm and 2130 nm Terra and Aqua bands, and a measure of temporal variability. Cumulative MODIS 1 km active fire detections are used to guide the selection of burned and unburned training samples and to guide the specification of prior burned and unburned probabilities. The MCD64A1 product is distributed in the standard MODIS Level 3 10 • × 10 • Land tile format in the sinusoidal projection [50].

MODIS Global Nadir BRDF-Adjusted Reflectance Product
The MODIS Nadir BRDF-Adjusted Reflectance (NBAR) product provides estimates of the nadiral surface reflectance at local solar noon performing a Bidirectional Reflectance Distribution Function (BRDF) inversion of MODIS surface reflectance for each day using a 16-day moving window centered on the nominal product date [51,52]. The Collection 6, Level 3 daily MODIS Terra and Aqua combined product (MCD43A4) was used in this work. The MCD43A4 product is defined in the MODIS Level 3 Land tile format in sinusoidal projection at 500 m resolution. Per-pixel quality assessment information is provided by the MCD43A2 product, also defined in the MODIS Level 3 Land tile format in sinusoidal projection at 500 m resolution and, for each reflective band of the MCD43A4 product, reports the quality of the BRDF inversion for each pixel of the relative MCD43A4 product, and flags the presence of water and snow pixels.

MODIS Land Cover Product
The MODIS Land Cover Type product (MCD12Q1) provides five land classification schemes, which describe land cover properties derived from one year of observations from Terra-and Aqua-MODIS [53]. The Collection 5.1, Level 3 yearly land cover product is defined in the MODIS Level 3 Land tile format in sinusoidal projection at 500 m resolution. The International Geosphere and Biosphere Programme (IGBP) scheme, which identifies 16 land cover classes, including 11 natural vegetation classes, 3 developed and 2 non-vegetated land classes, and has a reported 75% overall land cover classification accuracy [53], was used in this work.

Terrestrial Ecoregions of the World
The Terrestrial Ecoregions of the World (TEOW) map is a biogeographic division of the Earth's terrestrial biodiversity into 867 ecoregions, which belong to 14 biomes and 8 realms [54]. Ecoregions are defined as biogeographic units containing a homogeneous population of natural communities (flora and fauna) sharing a large majority of species, dynamics, and environmental conditions.

Theoretical Basis: Definition of the Burned Area Persistence Time
The burned area persistence time was defined as the maximum time after a fire during which burned pixels can be unambiguously separated from unburned ones (i.e., the length of the period in which burned areas can be reliably mapped).
The spectral signature of vegetation is altered by burning events, and the magnitude of these changes is greater in some portions of the electromagnetic spectrum than in others. It has been demonstrated that the spectral bands more suitable for burned area mapping are the near infrared (NIR) and shortwave infrared (SWIR) for both coarse- [55][56][57] and moderate-resolution satellites [58][59][60][61]. The spectral changes induced by fire are non-permanent; charcoal and ashes are removed by atmospheric agents exposing the bare soil, and vegetation regrows over time [62]. Additionally, the spectral signature of unburned pixels also changes over time: vegetation phenology, senescence and other disturbances, such as forest management thinning and clear cuts, mortality due to insect outbreaks and land use conversion, cause spectral changes that can be confused with those due to fire [63,64]. For these reasons, burned area spectral separability is sensitive to the time elapsed since the burning event [19,44,65]; the detectability of burned areas is a function-among other factors-of the burned area persistence time.
A wide variety of burned area mapping algorithms have been tested and developed for moderate resolution sensors including techniques that exploit single image analysis [66][67][68][69][70] and multi-temporal analysis [33,59,66,[71][72][73][74][75]. Although there is no agreement on the optimal burned area mapping algorithm, spectral indices are widely used in many studies. Among them, the Normalized Burn Ratio (NBR [76]) is widely used for burned area detection with moderate resolution satellite data [67,69,71,[73][74][75]77,78]. The index range is between −1 and 1, and it is defined as: where ρ N IR is the reflectance at 0.8 µm and ρ SW IR is the reflectance at 2.1 µm. The difference between pre-fire and post-fire NBR (dNBR) is an indication of the spectral changes due to fire. dNBR can assume values between −2 and 2. Because NBR drops after a fire due to the removal of vegetation, soil exposure and charcoal and ash deposition [58], positive dNBR values are associated with burning, thus making dNBR thresholding a simple burned area detection strategy [75,76,79]. In this study, daily time series of dNBR were generated for burned and unburned pixels from the Nadir BRDF-Adjusted Reflectance Albedo product (MCD43A4) to minimize the variability in reflectance due to different acquisition view geometry. Missing data, and water-and snow-covered pixels were discarded from further analysis.
For burned pixels at a generic location i that were detected as burned on Julian day of the year jdoy in the MODIS global burned area product (MCD64A1), the post-fire dNBR yearly time series was defined as: where NBR i,jdoy−8 , is the pre-fire NBR value, obtained from the MCD43A4 product with the nominal date 8 days before the day jdoy of burning, and NBR i,jdoy+∆t is the value observed ∆t days after the fire, with ∆t ≥ 8. The 16-day minimum time difference between pre-fire and post-fire observation was required to account for the 16-day inversion moving window used in the MCD43A4 product and to ensure that the NBR pre-fire and post-fire values were respectively calculated only from the inversion of pre-fire and post-fire reflectances. Generally, the spectral signature of vegetation varies during the year depending on phenology (e.g., senescence, leaf-off) and stress conditions; as a result, the NBR of unburned vegetation is variable during the year. Equation (2) generates yearly time series of dNBR measuring only the changes from the pre-fire conditions, regardless of the initial NBR value.
Similarly, a generic unburned pixel j, observed synchronously with the burned pixel i of Equation (2), will describe spectral temporal variations determined by vegetation phenology and other non-fire-related factors: where NBR j,jdoy−8 is the NBR value of the unburned pixel j on the same pre-fire date of Equation (2), and NBR j,jdoy+∆t is the value observed ∆t days after the starting day jdoy. Considering a set of burned pixels ba, all burning on the same day jdoy, the individual trends of Equation (2) were aggregated in a probability distribution function: This is illustrated in Figure 1. Considering a set of burned pixels detected on the same day jdoy in one spatial stratum of the analysis ('Victoria Plains tropical savanna' ecoregion in Northern Australia, savanna land cover), the histogram of the distribution P B (dNBR) ∆t,jdoy is represented by the dashed red lines. Notably, the range of dNBR values assumed by the burned area increases as time from the fire (i.e., ∆t), indicating heterogeneous post-fire trajectories, and the median value decreases, indicating a tendency to return to pre-fire conditions. The post-fire spectral signature varies depending on atmospheric events and vegetation regrowth following the fire; for example, wind and rain can scatter the charcoal and ashes, or a new layer of grasses can grow rapidly after a fire; these effects can have great variability also within the same burned area, resulting in a wider range of post-fire dNBR values.
Similarly, considering a set of unburned pixels ua, the individual trajectories of Equation (3) were also aggregated into a probability distribution: In Figure 1, the histogram of the distribution of P UB (dNBR) ∆t,jdoy for a set of unburned pixels in the same analysis stratum, and observed from the same day jdoy as the burned pixels, is represented by the dashed blue line. The range of dNBR values assumed by the unburned pixels also increases with ∆t, while the median remains stable around 0, indicating that the set of the unburned pixels shows increasingly heterogeneous trends as time passes; while some pixels show trends that cannot be confused with burning (i.e., negative dNBR changes), others do (i.e., positive dNBR changes), likely due to senescence and phenology.
The separability of burned/unburned pixels was measured in terms of the overlapping portion of the P B (dNBR) ∆t,jdoy and P UB (dNBR) ∆t,jdoy distributions. As a result of the temporal variation, the overlap of the two distribution is minimal immediately after the fire, and increases as time passes, reflecting the fact that burned areas can be mapped reliably only for a limited time after the fire. The persistence time of the burned area was, therefore, estimated as the maximum number of days in which the two distributions have minimal or no overlap, and more specifically the number of days in which at least 90% of the burned area pixels have dNBR higher than 90% of the unburned pixels. Using this empirical threshold, the burned area persistence time was estimated as the number of days on which the two distributions overlapped by 20%. the dashed red lines. Notably, the range of dNBR values assumed by the burned area increases as time from the fire (i.e., Δ ), indicating heterogeneous post-fire trajectories, and the median value decreases, indicating a tendency to return to pre-fire conditions. The post-fire spectral signature varies depending on atmospheric events and vegetation regrowth following the fire; for example, wind and rain can scatter the charcoal and ashes, or a new layer of grasses can grow rapidly after a fire; these effects can have great variability also within the same burned area, resulting in a wider range of post-fire dNBR values. Figure 1. Illustration of the method for the estimation of the burned area persistence time. The plot is generated using observations of a set of burned pixels (ba) and unburned pixels (ub) in the Savanna land cover of the "Victoria Plains tropical savanna" ecoregion in northern Australia, starting on the day of detection of the burned pixels (jdoy = 329, November 25th). The red dotted lines show the distribution of dNBR of the burned pixels ( (dNBR) , ), observed at Δ = 10, 20, 30 and 40 days after burning; the blue dotted lines show the distribution of dNBR of the unburned pixels ( (dNBR) , ) observed at the same time. The red and blue solid lines represent respectively the 10th percentile of the dNBR distribution of burned pixels (i.e., (Δ ) ) and the 90th percentile of the dNBR distribution of the unburned pixels (i.e., (Δ ) ). The persistence time is calculated as the time * when the two distributions overlap by more than 20%, i.e., when (Δ ) is equal to (Δ ) ( * = 31 in this example).

Figure 1.
Illustration of the method for the estimation of the burned area persistence time. The plot is generated using observations of a set of burned pixels (ba) and unburned pixels (ub) in the Savanna land cover of the "Victoria Plains tropical savanna" ecoregion in northern Australia, starting on the day of detection of the burned pixels (jdoy = 329, November 25th). The red dotted lines show the distribution of dNBR of the burned pixels (P B (dNBR) ∆t,jdoy ), observed at ∆t = 10, 20, 30 and 40 days after burning; the blue dotted lines show the distribution of dNBR of the unburned pixels (P UB (dNBR) ∆t,jdoy ) observed at the same time. The red and blue solid lines represent respectively the 10th percentile of the dNBR distribution of burned pixels (i.e., F B (∆t) jdoy ) and the 90th percentile of the dNBR distribution of the unburned pixels (i.e., F UB (∆t) jdoy ). The persistence time is calculated as the time ∆t * jdoy when the two distributions overlap by more than 20%, i.e., when F B (∆t) jdoy is equal to F UB (∆t) jdoy (∆t * jdoy = 31 days in this example).
By definition, 90% of the burned area pixels have a dNBR value higher than the 10th percentile of the P B (dNBR) ∆t,jdoy distribution, which, expressed as a function of ∆t, is defined as: Similarly, 90% of the unburned pixels have a dNBR value lower than the 90th percentile of the P B (dNBR) ∆t,jdoy : Finally, the persistence time of the burned areas was estimated as the number of days ∆t * jdoy in which F B (∆t) jdoy is greater or equal to F UB (∆t) jdoy : In Figure 1, F B (∆t) jdoy and F UB (∆t) jdoy are represented by the solid red and blue lines, respectively, and the persistence time ∆t * jdoy is defined by the intersection of the two lines.

Global Implementation
The probabilistic method described in Section 3.1 was implemented globally using the 14 years of MODIS datasets described in Section 2 as shown in Figure 2. Because of the analysis scale, it was Figure 2. Illustration of the procedure for the estimation of the burned area persistence time. The input data for the analysis are the MODIS burned area and NBAR product. The MODIS land cover product and the TEOW ecoregions are used to define each spatial stratum defined by ecoregion and land cover. The estimated burned date extracted from the MODIS burned area product is used to allocate the pixels in the different 16-day periods (Section 3.2.1). The median persistence time for each ecoregion and land cover is computed as the median value of the average annual burned area as a function of the persistence time (Section 3.2.3).

Spatial Stratification
Burned areas have different spectral responses depending on the type and condition of fuels, and more broadly depending on the ecosystem and time of burning [9,45,80], the analysis was performed using spatial and temporal units sufficiently fine to capture such variability. The analysis was stratified spatially adopting a two-level stratification. The Terrestrial Ecoregions of the World (TEOW) [54] were used to define the first stratification level. Ecoregions are delineated according to endemic genera and families (higher taxa), distinct assemblages of species, and the influence of geological history, such as past glaciations, on the distribution of plants and animals [54]. Consequently, they are more likely to accurately reflect the distribution of species and communities than alternative stratifications derived solely from biophysical variables, such as rainfall and temperature, vegetation structure, or spectral signatures from remote sensing data [54].
The TEOW map is a broad simplification of the variability of habitats, and all contain a variety of vegetation types [54]; for example, boreal ecoregions include forests, shrublands, and grasslands. Consequently, the MODIS MCD12Q1 land cover product was used to define a second level of Figure 2. Illustration of the procedure for the estimation of the burned area persistence time. The input data for the analysis are the MODIS burned area and NBAR product. The MODIS land cover product and the TEOW ecoregions are used to define each spatial stratum defined by ecoregion and land cover. The estimated burned date extracted from the MODIS burned area product is used to allocate the pixels in the different 16-day periods (Section 3.2.1). The median persistence time for each ecoregion and land cover is computed as the median value of the average annual burned area as a function of the persistence time (Section 3.2.3).

Spatial Stratification
Burned areas have different spectral responses depending on the type and condition of fuels, and more broadly depending on the ecosystem and time of burning [9,45,80], the analysis was performed using spatial and temporal units sufficiently fine to capture such variability. The analysis was stratified spatially adopting a two-level stratification. The Terrestrial Ecoregions of the World (TEOW) [54] were used to define the first stratification level. Ecoregions are delineated according to endemic genera and families (higher taxa), distinct assemblages of species, and the influence of geological history, such as past glaciations, on the distribution of plants and animals [54]. Consequently, they are more likely to accurately reflect the distribution of species and communities than alternative stratifications derived Remote Sens. 2018, 10, 750 8 of 32 solely from biophysical variables, such as rainfall and temperature, vegetation structure, or spectral signatures from remote sensing data [54].
The TEOW map is a broad simplification of the variability of habitats, and all contain a variety of vegetation types [54]; for example, boreal ecoregions include forests, shrublands, and grasslands. Consequently, the MODIS MCD12Q1 land cover product was used to define a second level of stratification within each TEOW ecoregion. The land cover classes of the IGBP classification scheme were aggregated into three major classes of interest (Forest, Shrubland, Grassland & Savanna), or masked out and removed from the subsequent analysis (urban areas, croplands, and miscellaneous non-burnable surfaces) ( Table 1). The spatial strata {Eco, LC}, defined by the generic land cover class LC within the generic ecoregion Eco, in which at least 250 km 2 of burned areas were cumulatively detected by the MCD64A1 product over the entire 14 year study period (2003-2016) were considered, thus excluding land cover class/ecoregion combination with negligible fire activity.

Temporal Stratification and Temporal Extent of the Analysis
Fire activity is mainly driven by fuel abundance and structure and by climatic conditions that result in changes in fuel moisture and lightning activity [81]. Consequently, different regions of the globe have distinct fire seasons. At the global scale, the fire season peak months in the Southern hemisphere are August-September, primarily driven by the extensive burning in Southern Africa and Australia, while in the Northern hemisphere the peak months are December-January, due primarily to burning across Northern and Central Africa [82]. At the regional scale, fire seasons vary in terms of length, peak activity and numbers of peaks [83]. These differences are linked mainly to the different climatic and the pre-burning fuel conditions, which vary during the fire season [84,85]. Generally, early season fires burn living, photosynthetically active vegetation, and dead vegetation that has a higher moisture content. During the dry season, vegetation senescence and climatic conditions lead to dryer living and dead vegetation [86,87]. As a result, fire-induced spectral changes for a given ecosystem and land cover vary within the fire season [85,88,89].
In order to account for these variations, a temporal stratification grid was created by dividing the January-December calendar year into 23 periods of 16 days each (t) (the last period being 13 days). This 16-day temporal grid-which coincides with the return interval of the Landsat satellite-was Remote Sens. 2018, 10, 750 9 of 32 sufficiently fine to characterize the spectral temporal variations of burned areas at different times of the year, while significantly reducing the computational load of a daily temporal grid.
To further reduce the computational load, only burned areas occurring during the fire season of each ecoregion were considered. The fire season of each ecoregion was defined by adapting the method proposed by Archibald, Lehmann, Gómez-Dans and Bradstock [84]. Each of the 16-day periods of the year was ranked based on the average annual burned area detected by the MCD64A1 product, and the fire season was defined as the union of the periods in which 90% of the annual burned area is detected. This definition requires no assumption on the continuity of the fire season and is applicable to ecoregions having a variety of temporal patterns of fire activity (Figure 3).
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 31 burned area is detected. This definition requires no assumption on the continuity of the fire season and is applicable to ecoregions having a variety of temporal patterns of fire activity (Figure 3).

Sampling Design
The MODIS MCD64A1 global burned area product detects, on average, ~4.0 × 10 6 km 2 of burned areas per year [16], corresponding to ~16 × 10 6 pixels per year. To make the analyses logistically possible, a sampling strategy was employed to sample both burned and unburned pixels in the spatial and temporal strata {Eco, LC, t} defined above as a combination of ecoregions, land cover classes, and time of the year quantized in 16-day periods.
The total sample size in each ecoregion was 100,000 burned and 100,000 unburned pixels (or the entire population if fewer than 100,000 pixels exist in either class). This total sample size was then allocated to the spatial sub-strata (i.e., land cover classes with fire activity in the ecoregion) and temporal strata (i.e., number of 16-day periods of the fire season) proportionally to the area burned, and the sample was extracted randomly with no replacement.
Under stratified random sampling, the probability of inclusion is constant within each stratum and is: where and are, respectively, the sample size and the population size in stratum h. It follows that within each stratum (i.e., each combination of ecoregion Eco, land cover LC, and 16-day period t) the inclusion probability was: where is the extracted sample, and is the total population of the burned area in the stratum. It should be noted that the sample , , was drawn from the total population of burned areas detected in ecoregion Eco and land cover LC, burned in the 14 years of the study period in the same 16 days of the calendar year defined by the interval t. center, bimodal fire season ("Victoria Plains tropical savanna" ecoregion in northern Australia); right, unimodal fire season spanning across calendar years ("Northern Congolian tropical forest-savanna mosaic" ecoregion in northern Africa). The fire season is identified by ranking the 16-day bins in decreasing order by burned area and including in the fire season the bins necessary to reach 90% of the burned area (shown in red). The 16-day periods outside the fire season (shown in blue) were not considered in the analysis.

Sampling Design
The MODIS MCD64A1 global burned area product detects, on average,~4.0 × 10 6 km 2 of burned areas per year [16], corresponding to~16 × 10 6 pixels per year. To make the analyses logistically possible, a sampling strategy was employed to sample both burned and unburned pixels in the spatial and temporal strata {Eco, LC, t} defined above as a combination of ecoregions, land cover classes, and time of the year quantized in 16-day periods.
The total sample size in each ecoregion was 100,000 burned and 100,000 unburned pixels (or the entire population if fewer than 100,000 pixels exist in either class). This total sample size was then allocated to the spatial sub-strata (i.e., land cover classes with fire activity in the ecoregion) and temporal strata (i.e., number of 16-day periods of the fire season) proportionally to the area burned, and the sample was extracted randomly with no replacement.
Under stratified random sampling, the probability of inclusion is constant within each stratum and is: where n h and N h are, respectively, the sample size and the population size in stratum h. It follows that within each stratum (i.e., each combination of ecoregion Eco, land cover LC, and 16-day period t) the inclusion probability was: where ba is the extracted sample, and BA is the total population of the burned area in the stratum. It should be noted that the sample ba Eco,LC,t was drawn from the total population of burned areas detected in ecoregion Eco and land cover LC, burned in the 14 years of the study period in the same 16 days of the calendar year defined by the interval t.

Ecoregion Level Estimation of the Persistence Time
The persistence time is estimated for each stratum, by applying Equations (1)-(8) using the sample of burned pixels ba Eco,LC,t and unburned pixels ub Eco,LC,t . It should be noted that the 16-day grid was used solely for the stratification of the sample, whereas the variation over time of the dNBR values (Equations (2) and (3)) was evaluated based on the daily MCD43A4 time series. As a consequence, for each land cover and ecoregion, up to 23 persistence times (one for each 16-day period t) were obtained, each of them estimated from Equation (8) with daily resolution. This is exemplified in Figure 4 (top), showing the burned area persistence time for the TEOW "Victoria plains tropical savanna" in Northern Australia. The 16-day periods included in the fire season are shaded in gray, and for each of them, the estimated persistence times ∆t * Eco,LC,t in "Shrubland" and "Grassland and Savannah" land cover classes are shown in gray and orange respectively. The persistence time was not computed for the "Forest" land cover, because no time period t had the minimum number of burned area detections defined in Section 3.2.1.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 31 grid was used solely for the stratification of the sample, whereas the variation over time of the dNBR values (Equations (2) and (3)) was evaluated based on the daily MCD43A4 time series. As a consequence, for each land cover and ecoregion, up to 23 persistence times (one for each 16-day period t) were obtained, each of them estimated from Equation (8) with daily resolution. This is exemplified in Figure 4 (top), showing the burned area persistence time for the TEOW "Victoria plains tropical savanna" in Northern Australia. The 16-day periods included in the fire season are shaded in gray, and for each of them, the estimated persistence times Δ * , , in "Shrubland" and "Grassland and Savannah" land cover classes are shown in gray and orange respectively. The persistence time was not computed for the "Forest" land cover, because no time period t had the minimum number of burned area detections defined in Section 3.2.1. In order to obtain meaningful summary metrics of the persistence time Δ * , , , it was necessary to first take into consideration the actual area burned in each time period t (reported in Figure 4, bottom) for the same ecoregion and land cover, and define an appropriate set of weights. This was done by estimating-in each ecoregion and land cover class-which proportion of the total burned area had a given persistence time Δ * , , . The average annual burned area with persistence time Δ * was estimated as follows: In order to obtain meaningful summary metrics of the persistence time ∆t * Eco,LC,t , it was necessary to first take into consideration the actual area burned in each time period t (reported in Figure 4, bottom) for the same ecoregion and land cover, and define an appropriate set of weights. This was done by estimating-in each ecoregion and land cover class-which proportion of the total burned area had a given persistence time ∆t * Eco,LC,t . The average annual burned area with persistence time ∆t * was estimated as follows: where n is the number of years used in the analysis, δ ∆t * Eco,LC,t (∆t * ) is the Dirac measure assuming value 1 if the persistence time for stratum {Eco, LC, t} is equal to ∆t * and 0 elsewhere, ba Eco,LC,t is the area of the sample extracted from each stratum, and π Eco,LC,t is the inclusion probability defined as in Equation (10).
where is the number of years used in the analysis, * , , (Δ * ) is the Dirac measure assuming value 1 if the persistence time for stratum {Eco, LC, t} is equal to Δ * and 0 elsewhere, , , is the area of the sample extracted from each stratum, and , , is the inclusion probability defined as in Equation (10).
Additionally, the interquartile range was computed as the difference between the 25th and 75th percentile of , (Δ * ), and used to measure the variability of the burned area persistence time within each spatial stratum. Figure 5 (bottom) shows the estimated median persistence time (Equation (12)) and the interquartile range for the ecoregion Victoria plains tropical savanna for the two land covers Shrubland and Grassland & Savanna. In this ecoregion, the estimated persistence time ∆ , * was longer for Shrubland land cover (55 days), and 50% of the burned area had an Additionally, the interquartile range was computed as the difference between the 25th and 75th percentile of BA Eco,LC (∆t * ), and used to measure the variability of the burned area persistence time within each spatial stratum. Figure 5 (bottom) shows the estimated median persistence time (Equation (12)) and the interquartile range for the ecoregion Victoria plains tropical savanna for the two land covers Shrubland and Grassland & Savanna. In this ecoregion, the estimated persistence time ∆t * Eco,LC was longer for Shrubland land cover (55 days), and 50% of the burned area had an estimated persistence between 42 and 74 days. In comparison, Grassland & Savanna burned area had a shorter persistence (43 days) and variability (interquartile range of 30-54 days).

Spatial Aggregation by Land Cover and BIOME/Realm
Because of the large number of spatiotemporal strata, it was necessary to summarize the persistence time results meaningfully in order to present and discuss them. The results were first summarized at the ecoregion level by use of the estimated median persistence time ( ∆t * Eco,LC , Equation (12)) and interquartile range, and then spatially aggregated by merging the ecoregions into larger spatial units based on Olson's biomes and realms (Australasia, Antarctic, Afrotropic, Indo-Malay, Neoarctic, Neotropic, Oceania, Paleoarctic) ( Figure 6).

Spatial Aggregation by Land Cover and BIOME/Realm
Because of the large number of spatiotemporal strata, it was necessary to summarize the persistence time results meaningfully in order to present and discuss them. The results were first summarized at the ecoregion level by use of the estimated median persistence time (∆ , * , Equation (12)) and interquartile range, and then spatially aggregated by merging the ecoregions into larger spatial units based on Olson's biomes and realms (Australasia, Antarctic, Afrotropic, Indo-Malay, Neoarctic, Neotropic, Oceania, Paleoarctic) ( Figure 6). Realms provide a subdivision of the main landmasses, and biomes are a convenient stratification unit because of their homogeneity of climate and vegetation and because "they provide a framework for comparisons among units and the identification of representative habitats and species assemblages" [54]. The TEOW biomes have been used to stratify studies for vegetation mapping [90] and animal species distributions [91], climatic models [92,93], anthropogenic urban growth [94], deforestation [95] and, recently, for the validation of remote sensing global burned area products [96,97]. Following the approach from Boschetti, Stehman and Roy [96], the 14 Olson biomes were aggregated into 5 more general biomes: Tropical, Temperate, Boreal, Mediterranean and Desert/Xeric biomes ( Table 2). Realms provide a subdivision of the main landmasses, and biomes are a convenient stratification unit because of their homogeneity of climate and vegetation and because "they provide a framework for comparisons among units and the identification of representative habitats and species assemblages" [54]. The TEOW biomes have been used to stratify studies for vegetation mapping [90] and animal species distributions [91], climatic models [92,93], anthropogenic urban growth [94], deforestation [95] and, recently, for the validation of remote sensing global burned area products [96,97]. Following the approach from Boschetti, Stehman and Roy [96], the 14 Olson biomes were aggregated into 5 more general biomes: Tropical, Temperate, Boreal, Mediterranean and Desert/Xeric biomes (Table 2). The realms of Oceania and Antarctica were excluded from the analysis because of the lack of significant fire activity. Of the 30 possible combinations of the remaining 6 realms and 5 biomes, 25 were valid (the Boreal biome is present only in the Paleoarctic and Neoarctic realms and the Mediterranean biome is not present in the Indo-Malay realm) and were adopted as aggregated spatial units for the summary of the analysis results ( Figure 6).
The distribution of annual burned area as a function of the spectral signal persistence time ∆t * for each spatial unit defined by Biomes (B), Realms (R), and land cover (LC) was simply calculated by summation of all the ecoregions belonging to the spatial unit, and globally by summation over the realms: The median persistence time for each Biome, Realm, and land cover ∆t * B,R,LC is the time corresponding to the 50th percentile of BA B,R,LC (∆t * ), and is computed as in Equation (12).

Fire Activity Distribution
Globally, from 2003 to 2016, the MODIS MCD64A1 product detected, on average, 4.0 × 10 6 km 2 of burned area per year: 0.19 × 10 6 km 2 in Forest, 0.31 × 10 6 km 2 in Shrubland, 3.1 × 10 6 km 2 in Grassland & Savanna, and the remaining 0.40 × 10 6 km 2 occurring in croplands and other land covers not considered in this study. The average annual area burned, reported by land cover and biome and the geographic distribution of area burned by land cover and ecoregion are displayed in Figure 7. The full result tables are presented as supplementary online material (Table S1). Out of the 867 ecoregions defined by the TEOW map, 511 ecoregions had more than 250 km 2 of burned area detected in at least one land cover class over the study period, and were considered in the analysis (Section 3.  The duration of the fire season was calculated in each ecoregion, as described in Section 3.2.1. (Table S2). Figure 8 shows the global distribution of the fire season duration at the ecoregion level. The median duration of the fire season of the ecoregions considered was 112 days (7 × 16-day periods), the shortest 32 days (2 periods) and the longest 304 days (19 periods). The duration of the fire season was calculated in each ecoregion, as described in Section 3.2.1. (Table S2). Figure 8 shows the global distribution of the fire season duration at the ecoregion level. The median duration of the fire season of the ecoregions considered was 112 days (7 × 16-day periods), the shortest 32 days (2 periods) and the longest 304 days (19 periods).

Ecoregion Level Burned Area Persistence Time
The burned area persistence time Δ , , * for each ecoregion (Eco), land cover class (LC) and 16-day period (t) of the fire season (Table S3), as well as the summary metrics , (Δ * ) and ∆ , * (Table S4)

Ecoregion Level Burned Area Persistence Time
The burned area persistence time ∆t * Eco,LC,t for each ecoregion (Eco), land cover class (LC) and 16-day period (t) of the fire season (Table S3), as well as the summary metricsBA Eco,LC (∆t * ) and ∆t * Eco,LC (Table S4) were estimated following the methodology described in Section 3.2.3. The median persistence time∆t * Eco,LC in each ecoregion and land cover, and the corresponding interquantile ranges, are displayed Figure 9.

Ecoregion Level Burned Area Persistence Time
The burned area persistence time Δ , , * for each ecoregion (Eco), land cover class (LC) and 16-day period (t) of the fire season (Table S3), as well as the summary metrics , (Δ * ) and ∆ , * (Table S4)    In Forest, 101 ecoregions (34%), 195 ecoregions (66%), and 230 ecoregions (77%) had a median persistence time shorter than 16, 32, and 48 days, respectively ( Figure 10, orange bars). Figure 9, top row, shows that the ecoregions with the shortest median persistence times are generally in the tropics, and those with the longest persistence times are in North America, the Western Mediterranean, and South Eastern Australia. The ecoregions with the highest interquartile range, which represents the persistence time variability during the fire season, are mostly located in Eurasia. In Forest, 101 ecoregions (34%), 195 ecoregions (66%), and 230 ecoregions (77%) had a median persistence time shorter than 16, 32, and 48 days, respectively ( Figure 10, orange bars). Figure 9, top row, shows that the ecoregions with the shortest median persistence times are generally in the tropics, and those with the longest persistence times are in North America, the Western Mediterranean, and South Eastern Australia. The ecoregions with the highest interquartile range, which represents the persistence time variability during the fire season, are mostly located in Eurasia.

Biome/Realm Aggregation of the Burned Area Persistence Time
The ecoregion-level results were aggregated into summary metrics for each land cover and Realm/Biome spatial unit, applying the methods described in Section 3.2.4 (Table S5). These results are presented in the following section, which is organized by land cover class.

Forest Land Cover
Forests are among the largest carbon pools of the globe [98], and forest fires are one of the main carbon sources for the atmosphere, contributing to 24.8% of the total global carbon emission from

Biome/Realm Aggregation of the Burned Area Persistence Time
The ecoregion-level results were aggregated into summary metrics for each land cover and Realm/Biome spatial unit, applying the methods described in Section 3.2.4 (Table S5). These results are presented in the following section, which is organized by land cover class.

Forest Land Cover
Forests are among the largest carbon pools of the globe [98], and forest fires are one of the main carbon sources for the atmosphere, contributing to 24.8% of the total global carbon emission from fires [24], even though forests contribute to just 4.75% (0.19 × 10 6 km 2 ) of the MODIS global annual area burned. Of this area, the majority was detected in Tropical (65.9%), Temperate (13.8%) and Boreal (19.1%) biomes, with a negligible contribution of the other biomes (Figure 7, top left). Although from an ecological perspective, fire effects persist for decades in most forests, the separability of the spectral signature of burned and unburned areas declines quickly. The global distribution of BA B,R,LC (∆t * ) (including all the 5 aggregated biomes) (Figure 11, bottom right) indicated that 66.4% of the annual burned area (0.09 × 10 6 km 2 per year) had a persistence time of less than 32 days, and less than 48 days for 79.4% of the annual burned area (0.10 × 10 6 km 2 per year); only 7.2% had a persistence time greater than 96 days. The biome-level distribution ofBA B,R,LC (∆t * ) is also displayed in Figure 11 and is complemented by the corresponding box plots reporting the median time∆t * B,R,LC and interquantile range (Figure 12). In Tropical biomes, the distribution ofBA B,R,LC (∆t * ) was unimodal, with short persistence times (80.5% of the area has a persistence under 32 days). In Temperate biomes, the distribution was largely bimodal: the first mode is largely due to burned areas detected in the Paleoarctic realm (where 28.0% of the area has a persistence of less than 16 days); the second mode, at ∆t * ≥ 96 days, was mainly due to burned area detected in Australasia and Neoarctic realm. In boreal biomes, the distribution is trimodal, due to different behavior in the Neoarctic and Paleoarctic realms. The Paleoarctic realm was characterized by a shorter persistence time (∆t * B,R,LC is 44 days) and bimodal distribution, with a first peak at ∆t * < 16 and a second peak at 48 ≤ ∆t * < 64. The Neoarctic had a longer persistence time (∆t * B,R,LC is >96 days) and unimodal distribution. As a result, the global distribution ofBA B,LC (∆t * ) ( Figure 11, bottom right) was a combination of short persistences, mainly from Tropical and Paleoarctic Temperate and Boreal burned areas, and long persistences, mainly from Australasia Temperate and Neoarctic Temperate and Boreal burned areas.  (Figure 7, top left). Although from an ecological perspective, fire effects persist for decades in most forests, the separability of the spectral signature of burned and unburned areas declines quickly. The global distribution of , , (Δ * ) (including all the 5 aggregated biomes) (Figure 11, bottom right) indicated that 66.4% of the annual burned area (0.09 × 10 6 km 2 per year) had a persistence time of less than 32 days, and less than 48 days for 79.4% of the annual burned area (0.10 × 10 6 km 2 per year); only 7.2% had a persistence time greater than 96 days. The biome-level distribution of , , (Δ * ) is also displayed in Figure 11 and is complemented by the corresponding box plots reporting the median time ∆ , , * and interquantile range (Figure 12). In Tropical biomes, the distribution of , , (Δ * ) was unimodal, with short persistence times (80.5% of the area has a persistence under 32 days). In Temperate biomes, the distribution was largely bimodal: the first mode is largely due to burned areas detected in the Paleoarctic realm (where 28.0% of the area has a persistence of less than 16 days); the second mode, at Δ * ≥ 96 days, was mainly due to burned area detected in Australasia and Neoarctic realm. In boreal biomes, the distribution is trimodal, due to different behavior in the Neoarctic and Paleoarctic realms. The Paleoarctic realm was characterized by a shorter persistence time (∆ , , * is 44 days) and bimodal distribution, with a first peak at Δ * < 16 and a second peak at 48 ≤ Δ * < 64. The Neoarctic had a longer persistence time (∆ , , * is >96 days) and unimodal distribution. As a result, the global distribution of , (Δ * ) ( Figure 11, bottom right) was a combination of short persistences, mainly from Tropical and Paleoarctic Temperate and Boreal burned areas, and long persistences, mainly from Australasia Temperate and Neoarctic Temperate and Boreal burned areas.

Shrubland Land Cover
The vast majority of the burned area in the Shrubland aggregated land cover class occurred in the Australasia (73.1%) and Afrotropic (18.7%) realms. In particular, it was a very significant class in the Australasian realm, where it contributed to 46.7% of the average annual burned area. The area burned was almost evenly split between the Tropical (41.6%) and Desert/Xeric (50.0%) biomes, with a negligible contribution of the other biomes (Figure 7, center left). Different persistence times were found in different biomes: shorter in the Tropical biomes, where the mode of the , , (Δ * ) histogram was between 16 and 32 days; and longer in the Desert/Xeric biomes, where the mode was greater than 96 days (Figures 13 and 14). The global histogram is bimodal, and reflects the distinct distribution of these two dominant biomes ( Figure 13).

Shrubland Land Cover
The vast majority of the burned area in the Shrubland aggregated land cover class occurred in the Australasia (73.1%) and Afrotropic (18.7%) realms. In particular, it was a very significant class in the Australasian realm, where it contributed to 46.7% of the average annual burned area. The area burned was almost evenly split between the Tropical (41.6%) and Desert/Xeric (50.0%) biomes, with a negligible contribution of the other biomes (Figure 7, center left). Different persistence times were found in different biomes: shorter in the Tropical biomes, where the mode of theBA B,R,LC (∆t * ) histogram was between 16 and 32 days; and longer in the Desert/Xeric biomes, where the mode was greater than 96 days (Figures 13 and 14). The global histogram is bimodal, and reflects the distinct distribution of these two dominant biomes ( Figure 13).

Grassland & Savanna Land Cover
The vast majority (93.4%) of the burned area for Grassland & Savanna was detected in the Tropical biome. Only 3.9% was detected in the Temperate biome, with a negligible contribution of the other biomes (Figure 7, bottom row). In particular, the majority of the global burned area detections were concentrated in the Tropical biome of the Afrotropic realm (76%). Globally, the mode of the , , (Δ * ) histogram ( Figure 15, bottom right) was between 16 and 32 days (60.5% of the annual burned area). Overall, the persistence time was short: 91% of the area burned had a persistence

Grassland & Savanna Land Cover
The vast majority (93.4%) of the burned area for Grassland & Savanna was detected in the Tropical biome. Only 3.9% was detected in the Temperate biome, with a negligible contribution of the other biomes (Figure 7, bottom row). In particular, the majority of the global burned area detections were concentrated in the Tropical biome of the Afrotropic realm (76%). Globally, the mode of the , , (Δ * ) histogram ( Figure 15, bottom right) was between 16 and 32 days (60.5% of the annual burned area). Overall, the persistence time was short: 91% of the area burned had a persistence  (Figures 15 and  16). At biome level, the histograms of , , (Δ * ) for Tropical and Temperate biomes were both unimodal. The histogram of , , (Δ * ) had its maximum between 16 and 32 days for Tropical biomes and between 32 and 48 days for Temperate biomes. A total of 92.7% of , , (Δ * ) for Tropical biomes and 80.5% for Temperate biomes was associated to Δ * between 16 and 48 days. The global histogram was unimodal, and reflects the distribution of the dominant Tropical biomes and a small contribution from the Temperate biomes ( Figure 15).    (Figures 15 and  16). At biome level, the histograms of , , (Δ * ) for Tropical and Temperate biomes were both unimodal. The histogram of , , (Δ * ) had its maximum between 16 and 32 days for Tropical biomes and between 32 and 48 days for Temperate biomes. A total of 92.7% of , , (Δ * ) for Tropical biomes and 80.5% for Temperate biomes was associated to Δ * between 16 and 48 days. The global histogram was unimodal, and reflects the distribution of the dominant Tropical biomes and a small contribution from the Temperate biomes ( Figure 15).

Ecoregion-Level Burned Area Persistence Time
The results at ecoregion level confirmed that the separability of burned and unburned surfaces quickly decreases after burning, even in ecosystems where the fire effects on the vegetation are persistent for several years. In addition to the dissipation of the charcoal and ashes, vegetation regrowth in burned areas, phenology of unburned vegetation, and several other factors might have contributed to these results.
In this study, the post-fire spectral temporal variations were analyzed by using NBR time series of burned and unburned pixels, as identified by the MCD64A1 global burned area product. Detection errors in the MCD64A1 product influenced the results, by reducing the separability of the burned (Equation (4)) and unburned (Equation (5)) probability distributions. More specifically, omission errors resulted in the inclusion of burned pixels in the unburned population ub, and biased positively the function F UB (Equation (7)). Conversely, commission errors caused the inclusion of unburned pixels in the burned population ba, biasing negatively the function F B (Equation (6)). As a result, both omission and commission errors in the MCD64A1 product resulted in a shorter estimated persistence time (Equation (8)).
Furthermore, the maximum persistence time ∆t * jdoy (Equation (8)) is limited by the length of the post-fire time series of nadir BRDF-adjusted reflectances (MCD43A4). Because at least 7 cloud-free MODIS observations within a 16-day period are needed to perform a valid BRDF inversion [52], in regions with clearly defined dry and wet seasons, or with persistent snow cover in the winter months, no MCD64A4 data were available for part of the year. If the fire season was close to the beginning of the wet or snow season, the estimated persistence time was effectively limited to the duration of the cloud-free period.
Cropland land cover was not considered in our analysis, since the MODIS burned area product is not suited to capturing the size and heterogeneity of cropland burned areas [18,99] because of its geometrical resolution, and because other agricultural practices (e.g., tilling) often immediately follow the fires [100], and have similar spectral signatures [101]. Omission errors in the cropland class of the MCD12Q1 land cover product, therefore, led to the inclusion of MCD64A1 detections over agricultural areas in the burned population ba. Some of these were genuine detections, but had a very short persistence time; some were commission errors and-as detailed above-negatively biased the F B function; in either case, the result was an underestimation of the persistence time. Additional errors in the MCD12 land cover map also influenced the estimated persistence time, by reducing the effectiveness of the stratification by aggregated land cover type.
Finally, caution should be used in considering the estimated burned area persistence time as defined in the present work; it should not be confused with the ecological concept of vegetation recovery to pre-fire level [102][103][104]. The duration of the persistence time can be limited by the recovery of NBR values to pre-fire levels after rapid re-vegetation or other factors, such as snowfall, that occlude the charcoal and bare soil and mask the spectral signature associated with burned areas. Additionally, errors in the MODIS burned area and land cover product could have influenced the stratification of the analysis and the division in the two population of burned and unburned pixels reducing the estimated persistence time.

Forest Land Cover
The majority of the burned area in Tropical biomes was detected in the Afrotropic, Indo-Malay and Neotropic realms (Figure 7, top left). In these three realms, the median persistence time∆t * B,R,LC was 8, 24 and 26 days, respectively (Figures 11 and 12, top left). Compared to the other biomes, tropical forests had shorter persistence times, arguably due to several factors. Many tropical forests have dense canopies that keep surface fuels moist and humidity high, limiting fire spread and duration, and allowing vegetation to recover quickly [105]. New photosynthetically active vegetation can alter the post-fire spectral signature sensed by EO satellites and deteriorate their ability to detect burned areas. For example, after a tropical forest fire in the Congo Basin, a new surface layer of vegetation established (Marantaceae), which covered the burned area in less than 3 months after the burning event [106]. Additionally, the estimated persistence time was influenced by the nature of tropical forest fires. Anthropogenic activities account for the majority of burning events in the tropical Amazon rainforest [107,108], African rainforests [89,109,110] and Indo-Australian rainforests [111]. In central Africa, human-ignited rainforest fires are strongly associated with land cover changes and are influenced by fire activity along the edges of the forest [89]. In Amazonia, fire activity is connected to other disturbances that have spectral similarities with burned areas, such as logging and deforestation, that only complex algorithms can discriminate [112]. In both of these cases, post-fire NBR may recover rapidly in areas planted immediately after clearing, where vegetation regrows within weeks [105]. In Southern China, various forest fire regimes were found [113]; regions characterized by long fire seasons and frequent fires, mainly due to human activity, show lower estimates of persistence time compared to regions characterized by shorter fire season and infrequent fires.
The majority of the burned area in Temperate biomes was detected in Australasia, Neoarctic and Neotropic realms (Figure 7, top left). In these three realms, the median persistence time∆t * B,R,LC was >96, 67 and 29 days, respectively ( Figures 11 and 12, top center). In Australasia, these results were consistent with the occurrence of high-intensity, low-frequency fires in dense, shade-tolerant understory vegetation in many temperate forests in Australia [84,114]. For the Neoarctic realm temperate forests, the median persistence time∆t * B,R,LC was the result of variable∆t * Eco,LC at the ecoregion level. The longest persistence time estimates occurred in mid-elevation, Northern Rockies ecoregions typically affected by large crown fires [115,116]. The shortest persistence times were in the southeast US, which is consistent with frequent surface fires and rapid post-fire grass regeneration observed in that region [117,118].∆t * B,R,LC estimated for the Paleoarctic realm temperate forests was greatly influenced by the presence of anthropogenic fires linked to land conversion for roads and croplands in the western China and Russia-China borders regions [119], in which the MODIS MCD64A1 product detects the majority of burned areas for this combination of realm and realm ( Figure 7, top right).
The Boreal biome is present only in the Neoartic and Paleoartic realms, accounting for 29.9% and 70.1% percent of the annual burned area respectively (Figure 7, top left). The estimated median persistence time∆t * B,R,LC was 44 days in the Paleoartic realm, and >96 days in the Neoarctic (Figures 11 and 12, top right). This large difference was likely due to the different burning conditions in these two realms. Fires in Paleoarctic boreal forests are dominated by low-intensity surface fires, and typically result in smaller-sized fires than Neoarctic boreal forests fires, where high-intensity crown fires are instead predominant [120,121]. Furthermore, the fire temporal distribution in the Paleoarctic realm was consistent with anthropogenic agricultural fires that typically occur in spring and fall and clear extensive areas of surface fuels [100,121]. The size and timing of Paleoarctic fire reduced the spectral separability of burned/unburned pixels in boreal biomes because of limited data quality (e.g., mixed pixels resulting from land-snow patterns and vegetation phenology) and availability (e.g., spatial resolution, cloud-free data availability) [46]. For the Boreal biomes, fire season is often followed closely by the snow season. As discussed in Section 4.2, the estimated persistence time was influenced by missing data due to snow cover, especially when the majority of burned areas were detected at the end of the fire season, such as in the Paleoarctic Boreal realm.

Shrubland Land Cover
The median persistence time∆t * B,R,LC for the Australasia realm was 55 days for the Tropical biome and >96 days for the Desert/Xeric biome. Tropical fires in Shrubland land cover are typically Remote Sens. 2018, 10, 750 23 of 32 characterized by a shorter persistence of char on the ground due to removal by wind and atmospheric agents; the observability is also limited by cloud cover immediately following the fire season [45]. For example, Bowman, Zhang, Walsh and Williams [44] found the temporal persistence of the burned areas in tropical regions of the Australian northern territories to be less than 100 days. Conversely, fire scars in arid zones are typically more persistent, as-once the charcoals and ash have dissipated-there is little vegetation regrowth until the wet season, leading to persistent soil exposure. The NBR value of soil is significantly different from unburned vegetation [122], leading to a longer estimated persistence time.
The median persistence time∆t * B,R,LC in the Afrotropic realm was 48 days and 37 days for the Tropical and the Desert/Xeric biomes, respectively. In both biomes, the persistence time had a large variability, with interquartile ranges of 37-77 days and 35-62 days (Figure 14, grey boxplot of Tropical and Desert/Xeric plots). Tropical shrubland fires are concentrated in the Sahelian zone in North equatorial Africa, including areas in East Africa, and the Kalahari region in South equatorial Africa [123]. In the central Kalahari, the majority of fire events can interest only herbaceous vegetation, while woody fuels remain unburned; only higher-intensity fires are able to burn even in the woody-cover dominated areas [124]. The NBR time series over burned areas showed variable changes in magnitude depending on the amount of unscorched vegetation within a MODIS pixel, and, as a result, the persistence time duration was variable. Finally, the persistence time large variability is potentially due to the different anthropogenic influence of fire activity over managed and protected land [125].

Grassland & Savanna Land Cover
Grasslands and savannas are fire-prone ecosystems, and fire is an essential factor contributing to the maintenance of the ecosystem; where resource availability and climatic control would allow tree seedling establishment, fire decreases tree cover density and maintains grassland and savanna ecosystems [8,126,127]. The short estimate and the small variability of the burned area persistence time for this land cover ( Figure 16, bottom right) could be explained by the common features that grassland and savanna share across the continents. Fine fuels and low fuel loadings result in a shorter persistence of the char residue signal; cloud cover and the frequent presence of extensive smoke aerosol layers result in a limited availability of satellite optical data, and the dominance of surface fires exacerbate the problem of detecting burns in woody savannas [45].
The majority of the burned area in Tropical biomes was detected in the Afrotropic, Australasia and Neotropic realms (Figure 7, bottom left). In these three realms, the median persistence time∆t * B,R,LC was 28, 42 and 27 days, respectively (Figures 15 and 16, top left). In the tropics and subtropics, the total amount of burned area and the interannual variability are controlled by the complex interaction between climate, fuel, and human activity [5,16,128]. The spectral characteristics of vegetation and burned areas change depending on the acquisition date [129] and are highly dynamic; rapid re-vegetation after a fire and the presence of unburned fuel loads can significantly impact the spectral signal [37,45]. The persistence time for Grasslands & Savanna was in accordance with other studies in tropical savannas. Scholes and Walker [130] reported a recovery to pre-fire albedo values 42 days after a fire in a southern African savanna. Likewise, ground-based reflectance was observed to be within 20% of pre-fire values within 14 days of burning events in western African savannas [41,42]. Eva and Lambin [131] analyzed the spectral reflectance temporal dynamics of burned and unburned woodland savannas in Central Africa using ATSR data and concluded that the discrimination between burned and unburned areas was affected by an error of 20% after 35 days from the fire event. The majority of the burned area persistence time for Tropical Australasia realm was longer than the Afrotropic realm ( Figure 16, top left). This difference was likely caused by the different dominant fire regimes of the two realms. In both realms, savanna fire frequency is annual/biannual; Australia has larger, more intense fires, whereas Africa is characterized by smaller, less intense fires [84]. Australian tropical savanna fire intensity rises as the dry season progresses due to more extreme fire weather and fire-prone fuel conditions in the later season [132][133][134]. A variable degree of change in the spectral signal caused by the action of fire was reported in Tropical grasslands and savannas [129], which is consistent with the variability of Grassland & Savanna persistence time results over Afrotropic and Australasian Tropical biomes ( Figure 16, top left).
The majority of the burned area in Temperate biomes was detected in the Afrotropic, Neoarctic, and Paleoarctic realms (Figure 7, bottom left). In these three realms, the median persistence timê ∆t * B,R,LC was respectively 30, 21 and 33 days (Figures 15 and 16, top center). Like other biomes, Grassland & Savanna fire activity in temperate biomes is mainly driven by climate and human activity. In North American temperate grassland ecoregions, antecedent precipitation amount and the drought index significantly influence fire activity as high-precipitation years can increase fuel load and subsequent fire intensity in the drought years that follow [135].
Agricultural fires detected in the MODIS land cover product were discarded in this work. Nevertheless, misclassified crop, managed prairie and pasture fires were likely included in the analysis and could reduce the estimated burned area persistence time, especially in ecoregions of the Paleoarctic realm where the use of fire is widespread as an agricultural management tool for the removal of excess residue and the control of diseases and pests [136]. Typically, temperate grasslands are intensively managed, including fire and grazing practices that alter nutrient cycling and the distribution of organic matter [127,137]. The relatively short burned area persistence time for temperate grassland was also observed in Canadian grassland fires. In these grasslands, Lu and He [138] estimated the post-fire vegetation recovery using multi-temporal Landsat imagery (one, two and three months after the fire) and concluded that grassland has a strong post-fire recovery capacity, especially if there is an adequate water availability. Prescribed fires are commonly used in tallgrass prairie ecosystems and are characterized by a relatively small-sized high degree of combustion heterogeneity, for this reason, MODIS' coarse resolution could be insufficient for an adequate mapping of this type of fire [139].

Conclusions
Coarse-resolution sensors typically have a daily temporal resolution, and it is generally assumed that a sufficient number of observations will be available for mapping burned areas globally. However, moderate-resolution satellites have reduced temporal resolutions (e.g., 16 days for Landsat, 5 days for Sentinel-2A/B), which could potentially lead to large burned area omission errors in ecosystems where the spectral signal associated with burning disappears quickly, especially if more than one post-fire observation is required to detect the burn. This paper attempted to calculate the burned area persistence time to estimate the period of time after the burn date in which burned areas are detectable using remotely sensed data change detection.
The presented methodology estimated the burned area persistence time defined as the duration of the spectral separability of burned/unburned areas, using time series of Normalized Burn Ratio (NBR) spectral index computed from MCD43A4 BRDF-adjusted MODIS reflectances. Burned and unburned areas were defined based on the MCD64A1 MODIS global burned area product. The full global MODIS record for the years 2003-2016 was used.
Among the novel aspects of the proposed methods, a probabilistic method for the analysis of spectral separability over time and space was developed, which can be extended to the analysis of other disturbances and, in general, of other non-permanent land use or land cover changes. The methods were stratified spatially (ecoregion and land cover) and temporally (time of the year) and provided for a rigorous generalization at coarser units.
The results showed that the persistence time was highly variable in time and space. Not only was it different across land cover classes, but within the same class, it depended on biome and realm. The shortest burned area persistence times were found over Tropical Forests, confirming that the separability of burned and unburned surfaces quickly decreases after burning, even in ecosystems where the fire effects on the vegetation are persistent for several years. The persistence times in Tropical