Double-Differenced dNBR: Combining MODIS and Landsat Imagery to Map Fine-Grained Fire MOSAICS in Lowland Eucalyptus Savanna in Kakadu National Park, Northern Australia

: A neglected dimension of the ﬁre regime concept is ﬁre patchiness. Habitat mosaics that emerge from the grain of burned and unburned patches (pyrodiversity) are critical for the persistence of a diverse range of plant and animal species. This issue is of particular importance in frequently burned tropical Eucalyptus savannas, where coarse ﬁre mosaics have been hypothesized to have caused the recent drastic population declines of small mammals. Satellites routinely used for ﬁre mapping in these systems are unable to accurately map ﬁne-grained ﬁre mosaics, frustrating our ability to determine whether declines in biodiversity are associated with local pyrodiversity. To advance this problem, we have developed a novel method (we call ‘double-differenced dNBR’) that combines the infrequent (c. 16 days) detailed spatial resolution Landsat with daily coarse scale coverage of MODIS and VIIRS to map pyrodiversity in the savannas of Kakadu National Park. We used seasonal Landsat mosaics and differenced normalized burn ratio (dNBR) to deﬁne burned areas, with a modiﬁcation to dNBR that subtracts long-term average dNBR to increase contrast. Our results show this approach is effective in mapping ﬁne-scale ﬁre mosaics in the homogenous lowland savannas, although inappropriate for nearby heterogenous landscapes. Comparison of this methods to other ﬁre metrics (e.g., area burned, seasonality) based on Landsat and MODIS imagery suggest this method is likely accurate and better at quantifying ﬁne-scale patchiness of ﬁre, albeit it demands detailed ﬁeld validation.


Introduction
The fire regime is a fundamental concept in fire ecology and organizing principle for fire management. Key elements include the temporal (frequency and seasonality), spatial (area burned, patterns of intensity, severity, and patchiness) metrics [1,2]. The availability of remote sensing products, particularly based on Landsat, MODIS and VIIRS sensors, has revolutionized the quantification and mapping of fire regimes, with the greatest emphasis on frequency, seasonality, area burned and increasingly severity. However, fire mosaics, the pattern of fine-scale fire history and time since fire, comprising patches with distinct severity or burn history, on the order of tens to hundreds of metres, remain understudied and poorly quantified due to the constraints of freely available satellite imagery used for fire mapping. MODIS and VIIRS have the virtue of monitoring fire activity at approximately daily time steps but are limited by their coarse spatial scales (c. 250-375 m) relative to the scale of fire patchiness. By contrast, Landsat has a comparatively high spatial resolution (30 m) but with a semi-monthly return interval, often extended by the prevalence of cloud cover. In tropical savanna environments, this coarse temporal resolution is particularly problematic because the regime of lower severity fire reduces the initial detectability of burn scars, which have little to no change to canopy reflectance, and the rapid post-fire 'regreening' of burned vegetation substantially distorts burned area patterns [3]. Previously, Landsat has been applied to quantification of fire patchiness and regime in Australian savannas at local scales [3][4][5][6][7], but routine broader-scale mapping is known to be relatively coarse compared to the actual grain size of fire mosaics [8].
The population collapse of the diverse assemblage of small mammals, birds and some plant species in north Australian savannas has been hypothesized as being related to the loss of fine-grained fire mosaics that were maintained by millennia of skillful Aboriginal fire management [9,10]. Recent years have seen the development of a fire management strategy in this system based around increased early dry-season burning [11,12]. An important motivation for fire managers to burn early in the dry season is to increase the number of unburned patches due to lower fire intensities and reduced fuel consumption. Early season burning is thought to benefit biodiversity, reduce greenhouse gas emissions, and potentially restore the habitat mosaics sustained by Indigenous fire management [8]. The benefits of large-scale early dry season burning programs, however, are debated and require further evaluation [13].
Here, we introduce a method that can combine the temporal frequency of MODIS with the finer spatial resolution of Landsat to improve detectability of fine-scale patch mosaics. Leveraging a readily available archive of seasonal-composite Landsat data [14], we apply this approach to the Eucalyptus savannas of Kakadu National Park (KNP) in the Northern Territory of Australia. Accurate mapping of fire mosaics in these ecosystems is essential to improve fire management programs across northern Australia as well as to understand the precipitous decline of wildlife species seen in recent decades [15,16] and to improve fire management programs across northern Australia, including KNP.
In this paper we describe an approach to improved fine-scale fire mapping in this system based on easily accessible seasonal composite Landsat data, and compare estimates of burned area totals, seasonality, and patch size statistics based on MODIS/VIIRS fire scars [17] (https://firenorth.org.au/nafi3, accessed on 8 March 2022), Landsat dNBR, and the new product combining both satellites' imagery we call 'double-differenced dNBR'.

Study Area
Kadadu National Park (Figure 1a) is located in the Northern Territory of Australia, and consists of extensive areas of tropical Eucalyptus savanna, with Melaleuca forests along wetlands and waterlines, grading into extensive sedgelands and mangrove systems towards the coast (Figure 1b). A prominent rocky plateau, the Arnhem Plateau, rises sharply from the south-east section of the park (Figure 1c).

Double-Differenced dNBR Methodology that Combines MODIS/VIIRS and Landsat
Our remote sensing approach has three stages: 1. Using seasonal Landsat composites, calculate dNBR for early and late dry season periods, and calculate the median over the shared MODIS/VIIRS fire scars (2000-2021). 2. Subtract these seasonal averages from each year's individual seasonal dNBR. 3. Classify burned and unburned areas using an experimentally determined double dNBR threshold.
Landsat reflectance seasonal composites were accessed from the TERN Seasonal Surface Reflectance data cube [14] which is available Australia wide, for the years 2000-2021, and collates data from a range of NASA Landsat satellites. These mosaics merge imagery in 3-monthly intervals on a standard grid, applying atmospheric and topographic corrections, and filling in cloud-affected pixels using medoid reflectance values where possible.
Three of the four available Landsat composite seasons were employed, March-May (S1), June-August (S2) and September-November (S3), with the December-February season excluded as it comprises the regional wet season. Fire seasons in the region typically run from April through to October, with annual variation depending on the timing of the monsoon, and we expect the S2-S1 difference to capture early season fires, and S3-S2 difference to capture late season fires. Pixels flagged as containing clouds that were unable to be removed by the seasonal compositing process were excluded from analysis. For each season we calculated the median normalized burn ratio, across all years, based on

Double-Differenced dNBR Methodology That Combines MODIS/VIIRS and Landsat
Our remote sensing approach has three stages:

1.
Using seasonal Landsat composites, calculate dNBR for early and late dry season periods, and calculate the median over the shared MODIS/VIIRS fire scars (2000-2021).

2.
Subtract these seasonal averages from each year's individual seasonal dNBR.

3.
Classify burned and unburned areas using an experimentally determined double dNBR threshold.
Landsat reflectance seasonal composites were accessed from the TERN Seasonal Surface Reflectance data cube [14] which is available Australia wide, for the years 2000-2021, and collates data from a range of NASA Landsat satellites. These mosaics merge imagery in 3-monthly intervals on a standard grid, applying atmospheric and topographic corrections, and filling in cloud-affected pixels using medoid reflectance values where possible.
Three of the four available Landsat composite seasons were employed, March-May (S1), June-August (S2) and September-November (S3), with the December-February season excluded as it comprises the regional wet season. Fire seasons in the region typically run from April through to October, with annual variation depending on the timing of the monsoon, and we expect the S2-S1 difference to capture early season fires, and S3-S2 difference to capture late season fires. Pixels flagged as containing clouds that were unable to be removed by the seasonal compositing process were excluded from analysis. For each season we calculated the median normalized burn ratio, across all years, based on composite Landsat pixel value, and the standard differential normalized burn ratio (dNBR [18]) formula using near-infrared (Landsat 8 band 5, 0.85 -0.88 µm) and shortwaveinfrared 2 (Landsat 8 band 7, 2.11-2.19 µm) bands, (B5 − B7)/(B5 + B7). All-year differential normalized burn ratios were then calculated between S1 and S2, and between S2 and S3 to represent a background burn ratio reflectance level for the early (oDNBR early ) and late (oDNBR late ) dry seasons.
For each year, we calculated dNBR between S1 and S2 (dNBR early ), and between S2 and S3 (dNBR late ), and then subtracted the associated seasonal dNBR baseline from each. Subtracting the seasonal composite has the effect of reducing noise from permanent features including topographic effects, vegetation, and soil patterns, thereby increasing the contrast of intermittent features which primarily comprise burn scars in this landscape. Pixels were classified based on a threshold (0.075) determined through exploratory data analysis. In this study, we use a single burnt/not-burnt threshold with the aim of detecting low severity fire in contrast to surrounding unburnt areas, but development of additional thresholds to differentiate severity classes may be possible with ground truthing. In identifying this threshold, we compared the classified burn patches to MODIS/VIIRS-derived burned area patch boundaries obtained from the North Australia and Rangelands Fire Information (NAFI) website [17], which are produced through a semi-automated process employing segmentation of imagery into similar spectral classes based on pre-and post-fire reflectance, with incorporation of active fire hotspots to improve mapping of small fires and cloudobscured areas. and using comparable seasonal definitions. If a pixel was classified as burned in the early dry season, it was excluded from burned area mapping in the late dry season on the assumption it cannot burn twice. A 5 × 5 modal smoothing filter was applied to the Landsat classification to reduce pixel noise, with the acknowledgement this may reduce the ability to detect extremely fine-scale patch dynamics.
Visual inspection of congruence between Landsat and MODIS/VIIRS burned area mapping highlighted that performance was poor in savannas on the rugged, topographically complex Arnhem Plateau likely due to persistent shadows and the presence of extensive areas of bare rock. Subsequent analysis was thus constrained to the lowland Eucalyptus savannas (Figure 2a). Processing was performed in R version 4.2.0 [19] using the raster package [20] for geographic raster processing.
A third combined burned area classification was constructed by clipping the Landsat classification to the burn scars identified in the MODIS/VIIRS mapping, removing burned areas not detected by the coarse classification while retaining unburned holes identified by Landsat within the coarse classification. The intent of this simple, combined method is to reduce noise from false-positive patch detections outside known coarse burn scars, while retaining fine-scale unburnt areas (which are important wildlife habitat) within coarse-scale burn scares that are undetectable with the coarse-scale resolution. A third combined burned area classification was constructed by clipping the Landsat classification to the burn scars identified in the MODIS/VIIRS mapping, removing burned areas not detected by the coarse classification while retaining unburned holes identified by Landsat within the coarse classification. The intent of this simple, combined method is to reduce noise from false-positive patch detections outside known coarse burn scars, while retaining fine-scale unburnt areas (which are important wildlife habitat) within coarse-scale burn scares that are undetectable with the coarse-scale resolution.

Statistical Analysis
For each year within the study area, a confusion matrix with classification statistics was generated to compare Landsat with MODIS/VIIRS burn area pixels. MODIS/VIIRS burned areas were downscaled to the resolution of the Landsat data, and Kappa,

Statistical Analysis
For each year within the study area, a confusion matrix with classification statistics was generated to compare Landsat with MODIS/VIIRS burn area pixels. MODIS/VIIRS burned areas were downscaled to the resolution of the Landsat data, and Kappa, sensitivity, specificity and simple accuracy statistics were calculated using R's caret package [21]. Means and standard deviations of these statistics were calculated across all years.
We calculated the total area of classified burn pixels for each year for the Landsat and MODIS/VIIRS data for early, late, and annual dry seasons. We then compared the total burned area of both Landsat and MODIS satellites.
A patch, in the context of spatial analysis, consists of a contiguous set of cells of the same class. A set of landscape patch metrics, comprising mean patch area, coefficient of variation (C.V.) of patch area, number of patches and total area were calculated for each year for the Landsat, MODIS/VIIRS, and combined burn area classifications operating at the native pixel-scale of burnt and unburnt areas for each classification using the landscapemetrics [22] package in R, and means and standard deviations calculated across years. One-way ANOVAs on log-transformed landscape metric values and Tukey HSD tests were performed using R's multcomp package [23] to detect significant differences between satellites within seasonal subsets. Selected landscape metrics (mean patch size, number of patches, and contagion index) were also calculated on time since fire (TSF) rasters generated from the 2000-2021 fire history for the MODIS/VIIRS, Landsat, and combined methods, and for a subset of long-unburned (>5 years) patches.

Results and Discussion
We have developed a new approach we call 'double-differenced dNBR' to help researchers and fire managers more accurately map fire patchiness in north Australian savannas. This approach should be transferable to comparable savanna landscapes where fire frequency is high, post-fire vegetation recovery is rapid, and which have large areas of homogenous vegetation.
Standard dNBR calculated on the seasonal composite Landsat imagery contained visible, static high-contrast permanent features, generally corresponding to topographically complex areas and waterlines, as demonstrated in Figure 2b for a subset of the study areas. Double-differenced dNBR visibly reduced the contrast of those permanent features, leaving intermittent fire scars visible (Figure 2c). A comparison of MODIS/VIIRS burned area polygons with Landsat classification for the double-difference dNBR shows that Landsat burned areas generally fell within larger MODIS/VIIRS burned areas but revealed reduced boundary area and internal patchiness not evident in MODIS/VIIRS (Figure 2d).
Given the contrasting temporal and spatial resolution of MODIS/VIIRS and Landsat, the burned area product of the former consistently over-estimated the area burned in the lowland Eucalyptus savannas of Kakadu National Park compared to the finer-scale fire scar mapping of the latter (Figure 3). Alignment in total area burned was greatest in the late dry season, where more intense fires and less post-fire recovery produces more easily detectable fire scars, an effect also suggested by [4] in their application of a Landsat time series analysis in Kakadu National Park. The congruence between the two across all years as measured by simple accuracy is maximal in the late dry season (86.7%, SD 5.4%), with an early dry season congruence of 62.8% (SD 5.9%) and an annual congruence of 58.1% (SD 6.1%). The relatively high values of the simple accuracy score may reflect the high prevalence of true negatives in the data, with both methods classifying similar areas as unburned. The Kappa statistic which gives a more conservative measure of agreement suggests classification agreement between the two satellite methods was poorer, with the highest score of 0.41 (SD 0.16) again achieved for the late dry season year, and lower agreement for both the early dry season (0.23 SD 0.09) and the full year (0.22 SD 0.06). This suggests that large areas classified as burned by MODIS/VIIRS were classified as unburned in the Landsat imagery as finer-scale patchiness became detectable. Sensitivity scores were higher for the early (92.8% SD 4.8%) and late (93.2% SD 5.0%) seasons than for analysis of the year as a whole (86.0% SD 7.1%), but all indicate good identification of MODIS burned pixels within the double-dNBR Landsat product. Specificity was lower and variable between years, with scores of 29.6% (SD 30.9%) for the early season, 46.2% (SD 18.4%) for the late season, and 41.2% (SD 12.7%) for the annual comparison, confirming the utility of the combined method in clipping Landsat areas outside known MODIS boundaries to exclude unburnt areas incorrectly classified as burnt. confirming the utility of the combined method in clipping Landsat areas outside known MODIS boundaries to exclude unburnt areas incorrectly classified as burnt. The mismatch in spatial resolution between these MODIS/VIIRS and Landsat satellite products affects estimates of all fire metrics (Figure 3), an expected product of resolution driving accuracy observed in previous comparisons of spatial resolution in fire classification [24] . When we combined both MODIS/VIIRS-and Landsat-classified fire patches, however, the landscape metrics reported results closely aligned more closely with the landscape metrics for Landsat imagery alone, suggesting our double dNBR approach correctly identifies fine-scale fire patch dynamics in Landsat imagery. Mean patch area (Figure 4a) was smaller in Lansdat and combined approach than in MODIS/VIIRS, as would be expected from the significantly larger pixel size in the later. The total burned area (Figure 4b) was lower using the Landsat or combined methods, but results were more strongly aligned to the late dry season MODIS/VIIRS burned area. The coefficient of variation in patch area (Figure 4c) was similar among methods, with variation slightly lower in MODIS/VIIRS, particularly in the late dry season. Further, the combined approach shows positive and negative variation in the magnitude of fire metrics compared to Landsat alone: Larger but fewer patches (Figure 4d) covering a smaller area. These trends hold for both the early and late dry season but are stronger in the late dry season. One-way ANO-VAs performed between satellites were significant at a p < 0.05 level in all cases, except for the coefficient of variation in patch area for the early dry season. Closer concurrence in burned area between the coarse-resolution and fine-resolution methodologies during the late dry season concurs with the findings of [6] for the smaller Landsat subset of the broader landscape they analysed. The mismatch in spatial resolution between these MODIS/VIIRS and Landsat satellite products affects estimates of all fire metrics (Figure 3), an expected product of resolution driving accuracy observed in previous comparisons of spatial resolution in fire classification [24]. When we combined both MODIS/VIIRS-and Landsat-classified fire patches, however, the landscape metrics reported results closely aligned more closely with the landscape metrics for Landsat imagery alone, suggesting our double dNBR approach correctly identifies fine-scale fire patch dynamics in Landsat imagery. Mean patch area (Figure 4a) was smaller in Lansdat and combined approach than in MODIS/VIIRS, as would be expected from the significantly larger pixel size in the later. The total burned area (Figure 4b) was lower using the Landsat or combined methods, but results were more strongly aligned to the late dry season MODIS/VIIRS burned area. The coefficient of variation in patch area (Figure 4c) was similar among methods, with variation slightly lower in MODIS/VIIRS, particularly in the late dry season. Further, the combined approach shows positive and negative variation in the magnitude of fire metrics compared to Landsat alone: Larger but fewer patches (Figure 4d) covering a smaller area. These trends hold for both the early and late dry season but are stronger in the late dry season. One-way ANOVAs performed between satellites were significant at a p < 0.05 level in all cases, except for the coefficient of variation in patch area for the early dry season. Closer concurrence in burned area between the coarse-resolution and fine-resolution methodologies during the late dry season concurs with the findings of [6] for the smaller Landsat subset of the broader landscape they analysed.  Landscape metrics calculated for the final year of data (2021) reveal differences in the spatial pattern and size of patches of different time since fire, including long unburned (>5 years) patches ( Table 1). Patches of a given age were measured as smallest and most numerous in the Landsat classification and significantly larger and fewer in the MODIS/VIIRS classification, with the combined index intermediate but closer to Landsat, as false-positive detections outside confirmed MODIS/VIIRS boundaries were reduced. The contagion index, which represents the tendency for a landscape to be composed of a few large patches or many smaller patches, was lower for the Landsat and combined methods than for MODIS/VIIRS, reflecting the ability of the finer-scale satellite imagery to detect a finer-scale mosaic of fire history. Long-unburned patches presented a similar pattern, with MODIS/VIIRS classifying fewer, larger long-unburned patches, and both Landscape and combined satellites detecting a greater number of patches with a smaller mean area. The mean time since fire, indicative of average return time, was just 2 years when calculated from MODIS/VIIRS imagery, but closer to 5 years using the combined approach. Landscape metrics calculated for the final year of data (2021) reveal differences in the spatial pattern and size of patches of different time since fire, including long unburned (>5 years) patches ( Table 1). Patches of a given age were measured as smallest and most numerous in the Landsat classification and significantly larger and fewer in the MODIS/VIIRS classification, with the combined index intermediate but closer to Landsat, as false-positive detections outside confirmed MODIS/VIIRS boundaries were reduced. The contagion index, which represents the tendency for a landscape to be composed of a few large patches or many smaller patches, was lower for the Landsat and combined methods than for MODIS/VIIRS, reflecting the ability of the finer-scale satellite imagery to detect a finerscale mosaic of fire history. Long-unburned patches presented a similar pattern, with MODIS/VIIRS classifying fewer, larger long-unburned patches, and both Landscape and combined satellites detecting a greater number of patches with a smaller mean area. The mean time since fire, indicative of average return time, was just 2 years when calculated from MODIS/VIIRS imagery, but closer to 5 years using the combined approach. The combined satellites using double-differenced dNBR to remove permanent non-fire features appears to work well in the relatively flat lowland Eucalyptus savannas, refining the broader MODIS/VIIRS burned areas and revealing unburned patchiness that may be overcounted be the coarser MODIS/VIIRS spatial resolution. The topographic complexity of the rocky highlands of the Arnhem Plateau, however, were not able to be overcome by the double-differenced dNBR method, however instead of being applied on a broad landscape scale the method has been successfully used to map individual fire scars in small patches of rock-free savannas that occur on the plateau [25]. It should, therefore, be possible, with pre-masking of rocky areas based on topographic boundaries of spectral signature, to apply this methodology in more topographically complex landscapes. Furthermore, the threshold we identified for classifying burned vs. unburned pixels (0.075) was manually determined by comparing MODIS/VIIRS burn scars with our classifications. Because this is specific to our study region, different local thresholds may need to be identified and applied elsewhere dependent on, e.g., tree density or surface reflectance.
Despite the use of seasonal composites and the low cloud cover present during the northern Australian dry season, cloud cover still negatively impacted the Landsat imagery: MODIS/VIIRS-identified burn scars could not be aligned to composite Landsat imagery in some years due to the prevalence of cloud cover. This has, to some extent, driven the lower burned area total detected with Landsat. Given this, the optimal application of the combined MODIS/VIIRS and Landsat method may be to investigate patchiness of select fires rather than investigating burned area on a coarse landscape scale. The inclusion of additional remote sensing platforms in recent years (e.g., Landsat-9, Sentinel 2), however, should greatly improve detectability of fire and minimize cloud interference. While this study focuses on the application of a double-differencing technique to the dNBR algorithm, as this is most frequently applied to burnt area and severity mapping in Australia, other methods based on pre/post-fire spectral differencing, for example dNDVI (normalized differential vegetation index [26]) or dNDMI (normalized differential moisture index [27]) may benefit equally from the contrast enhancement.
Reliable mapping is essential for monitoring regional fire management practices to ensure the appropriate provision of unburned fire refugia and reduction carbon emissions. High resolution mapping of pyrodiversity is essential for understanding the drivers of precipitous declines in biodiversity in north Australian savannas. Double-differenced dNBR that combined Landsat and MODIS/VIIRS appears to delineate fine-scale patchiness in savanna fires burning more accurately than Landsat alone. However, before double-differenced dNBR can be made routine in northern Australian lowland savannas, or comparable homogenous savanna ecosystems, the approach requires comprehensive ground truthing. Nonetheless, the approach has been effectively used to quantify pyrodiversity in small area of sandsheet savanna that is structurally comparable to lowland Eucalyptus savannas on the Arnhem Plateau [25]. This study found that double-differenced dNBR was ecologically meaningful because it was concordant with ecological monitor- ing that showed the loss of fine-scale mosaics which caused population declines of a fire sensitive conifer (Callitris intratropica).

Conclusions
Our double-differencing approach leverages the strengths of two mainstream satellite products, each with extensive record lengths, to improve mapping of fine-scaled patchiness in a savanna ecosystem. While remote sensing products are continuously improving in their spatial and temporal coverage, resolution, and data availability, the double-dNBR method has utility for retrospective analysis using available data products. All fire mapping approaches, require field validation to confirm the detection of legitimate ecological impacts, and we believe our improved fine-scale mapping method will stimulate field research into mapping fine-scale burn patterns, which is a neglected aspect of savanna ecology. Unmanned aerial vehicles are increasingly employed in mid-scale ecological assessments and we believe they are a prospective technology to achieve this in this system and validate satellite-derived mapping. Improved fire-scar mapping is essential to advance understanding the nexus of pyrodiversity and biodiversity, especially in regard to wildlife habitat condition and long-unburned refugia. However, additional research is required to improve the application of this method to more heterogeneous environments, potentially by masking confounding landscape features based on spectral classification.

Conflicts of Interest:
The authors declare no conflict of interest.