Mapping Wetland Burned Area from Sentinel-2 across the Southeastern United States and Its Contributions Relative to Landsat-8 (2016–2019)

: Prescribed ﬁres and wildﬁres are common in wetland ecosystems across the Southeastern United States. However, the wetland burned area has been chronically underestimated across the region due to (1) spectral confusion between open water and burned area, (2) rapid post-ﬁre vegetation regrowth, and (3) high annual precipitation limiting clear-sky satellite observations. We developed a machine learning algorithm speciﬁcally for burned area in wetlands, and applied the algorithm to the Sentinel-2 archive (2016–2019) across the Southeastern US (>290,000 km 2 ). Combining Landsat-8 imagery with Sentinel-2 increased the annual clear-sky observation count from 17 to 46 in 2016 and from 16 to 78 in 2019. When validated with WorldView imagery, the Sentinel-2 burned area had a 29% and 30% omission and commission rates of error for burned area, respectively, compared to the US Geological Survey Landsat-8 Burned Area Product (L8 BA), which had a 47% and 8% omission and commission rate of error, respectively. The Sentinel-2 algorithm and the L8 BA mapped burned area within 78% and 60% of wetland ﬁre perimeters ( n = 555) compiled from state and federal agencies, respectively. This analysis demonstrated the potential of Sentinel-2 to support efforts to track the burned area, especially across challenging ecosystem types, such as wetlands.


Introduction
Wetland ecosystems, characterized by permanent or seasonal saturation with water, hydromorphic soils, and aquatic plants [1], play critical global functions in providing clean water, supporting biodiversity, and regulating climate [2]. Wetlands are essential carbon sinks where saturated soils create anoxic conditions that limit decomposition rates and increase subsurface carbon storage [3][4][5][6]. Thus, although wetlands cover only 6% of the Earth's land cover [7], they store between 20 and 30% of the global soil carbon [8]. Wetlands in the conterminous US, for example, store a total of 11.52 PgC, with palustrine wetlands storing nearly ten-fold more carbon than tidal wetlands [9]. Despite typically high fuel moisture, prescribed fire in wetlands is frequently used to control plant biomass and non-native species, as well as to improve habitat. Wildfires in wetlands are also common, especially during dry periods across the Everglades, Brazilian Pantanal, peat bogs in Indonesia, boreal forests, and other global wetlands [10][11][12]. Fires drive key ecosystem processes (e.g., nutrient cycling, seed scarification, and plant growth) [13], but also consume biomass, emit particulate matter that impacts air quality, and release carbon dioxide to the In this study, we developed an algorithm to map the burned area extent in wetland ecosystems by using the Sentinel-2 image archive across Florida and the coastal portions of Alabama, Georgia, and South Carolina. The USGS Landsat BA [34] is currently the finest resolution burned area product operationally produced across the US. However, Landsat-7 is 22 years into its mission and is likely to be de-commissioned soon. It is critical, therefore, to explore the extent to which Sentinel-2 may supplement the Landsat-8 (L8) BA, particularly in regions with rapid regrowth such as the Southeastern US. To do so, we addressed the following research questions: (1) What is the contribution of clear-sky observations from Sentinel-2, relative to Landsat-8? (2) Can Sentinel-2 effectively map burned areas in wetland ecosystems? (3) What are the individual and combined contributions of Sentinel-2 and Landsat-8 to burned area extent in wetlands, and how do they compare to verified wetland fire perimeters?

Study Area
The study was focused on a subset of the Southeastern US, defined here as Florida, and the coastal portions of Georgia, South Carolina, and Alabama (up to 200 km from the coast, total area = 293,821 km 2 ) ( Figure 2). The region has a humid, subtropical climate. Annual precipitation averages 1315 mm, and annual temperature maximum and minimum average 26 and 14 °C, respectively [54]. Across the region, emergent and forested In this study, we developed an algorithm to map the burned area extent in wetland ecosystems by using the Sentinel-2 image archive across Florida and the coastal portions of Alabama, Georgia, and South Carolina. The USGS Landsat BA [34] is currently the finest resolution burned area product operationally produced across the US. However, Landsat-7 is 22 years into its mission and is likely to be de-commissioned soon. It is critical, therefore, to explore the extent to which Sentinel-2 may supplement the Landsat-8 (L8) BA, particularly in regions with rapid regrowth such as the Southeastern US. To do so, we addressed the following research questions: (1) What is the contribution of clear-sky observations from Sentinel-2, relative to Landsat-8? (2) Can Sentinel-2 effectively map burned areas in wetland ecosystems? (3) What are the individual and combined contributions of Sentinel-2 and Landsat-8 to burned area extent in wetlands, and how do they compare to verified wetland fire perimeters?

Study Area
The study was focused on a subset of the Southeastern US, defined here as Florida, and the coastal portions of Georgia, South Carolina, and Alabama (up to 200 km from the coast, total area = 293,821 km 2 ) ( Figure 2). The region has a humid, subtropical climate. Annual precipitation averages 1315 mm, and annual temperature maximum and minimum average 26 and 14 • C, respectively [54]. Across the region, emergent and forested wetlands are the dominant land-cover type (32%, Figure 2a), followed by evergreen forests (22%) and agriculture (19%, cultivated crops, hay/pasture) [55]. Wetland extent within this study was defined by using the National Wetland Inventory (NWI, version 2.0) dataset. Wetland extent averaged 29 ha km −2 across the region with wetland area dominated by freshwater forested/shrub wetlands (50%), estuarine and marine wetlands (21%), and freshwater emergent wetlands (11%) [56]. wetlands are the dominant land-cover type (32%, Figure 2a), followed by evergreen forests (22%) and agriculture (19%, cultivated crops, hay/pasture) [55]. Wetland extent within this study was defined by using the National Wetland Inventory (NWI, version 2.0) dataset. Wetland extent averaged 29 ha km −2 across the region with wetland area dominated by freshwater forested/shrub wetlands (50%), estuarine and marine wetlands (21%), and freshwater emergent wetlands (11%) [56].

Clear-Sky Contributions of Sentinel-2
Clear-sky observation counts were calculated in Google Earth Engine [57]. All Landsat-7 (non-gap filled) and 8 surface reflectance, and Sentinel-2 top-of-atmosphere images collected in 2016-2019 across the study area were used. Values identified as cloud or cloud shadow in Landsat-7 and Landsat-8 images were masked using cFMask [58]. Cloud and cloud shadow masks were applied to the Sentinel-2 image collection, as described in Section 2.3.1. The Normalized Difference Vegetation Index (NDVI) was calculated for all unmasked pixels across the three image collections. The per-pixel observation count of NDVI was summed annually and for March/April of each year (2016-2019). This two-month period represents the seasonal peak in fire activity [24]. For multisensor comparisons, the per-pixel observation count was summed across sensors. Because Landsat-7 is still active, contributing clear-sky observations from this sensor were also calculated for comparison purposes.

Clear-Sky Contributions of Sentinel-2
Clear-sky observation counts were calculated in Google Earth Engine [57]. All Landsat-7 (non-gap filled) and 8 surface reflectance, and Sentinel-2 top-of-atmosphere images collected in 2016-2019 across the study area were used. Values identified as cloud or cloud shadow in Landsat-7 and Landsat-8 images were masked using cFMask [58]. Cloud and cloud shadow masks were applied to the Sentinel-2 image collection, as described in Section 2.3.1. The Normalized Difference Vegetation Index (NDVI) was calculated for all unmasked pixels across the three image collections. The per-pixel observation count of NDVI was summed annually and for March/April of each year (2016-2019). This two-month period represents the seasonal peak in fire activity [24]. For multisensor comparisons, the per-pixel observation count was summed across sensors. Because Landsat-7 is still active, contributing clear-sky observations from this sensor were also calculated for comparison purposes.

Burned Area Algorithm Development
For burned area, we trained and validated the algorithm entirely on burned pixels within wetlands identified by using the National Wetland Inventory dataset and limited to estuarine and marine wetlands, as well as freshwater emergent, shrub, and forested wetlands [56]. For unburned areas, we trained and validated the algorithm across both wetland and non-wetland land-cover types. Because wetlands often occur in a mosaic with other land cover types, fire events frequently extend beyond a wetland boundary. To avoid mapping only parts of fire events, the burned area algorithm was applied across the entire study area.

Image Collection Preprocessing
All Sentinel-2 top-of-atmosphere (TOA) image collection (Level 1C) images acquired within the study area between 1 January 1 2016, and 31 December 31 2019, were processed. To mask clouds, we used the Sentinel-2 cloud probability layer [59], where pixels with a probability of >50% were classified as clouds and were masked out. Because the Level 1C image collection lacks a cloud shadow mask, cloud shadows were identified and masked by modifying code provided through the Google Earth Engine Community, "Sentinel-2 Cloud Masking with s2cloudless". Cloud shadows were defined as the intersection of the predicted cloud shadow area, derived from the cloud probability layer and the mean solar azimuth angle, and dark areas, identified as low-reflectance, near-infrared (NIR) pixels. To compensate for potential atmospheric contamination, we avoided using single spectral bands as variables and instead relied on normalized spectral indices. Google Earth Engine sets image count limits to control memory use. As Sentinel-2 is delivered as tiles, image stacks can rapidly exceed memory limits. Consequently, all images collected on the same date were mosaicked into a single image to reduce image counts.

Training Data
Training pixels were sampled from Sentinel-2 training images selected to represent four subregions across the Southeast that experience high amounts of wetland fire: (1) the Florida panhandle, (2) the Florida Everglades, (3) the coastal plains of South Carolina, and (4) the Eastern Florida-Georgia border ( Figure 2). In each training image, burned points were randomly sampled within burn perimeters mapped by or created from (1) the Geospatial Multi-Agency Coordination (GeoMAC) [60], now National Incident Fire Services (NIFS), (2) the Wildland Fire Decision Support System (USGS WFDSS, https://wfdss.usgs.gov/wfdss/WFDSS_Data.shtml (accessed on 12 May 2021), or (3) clusters of points identified by the MODIS Active Fire Product (MCD14) [20]. All burned points were restricted to palustrine emergent, palustrine forested/scrub-shrub, or estuarine and marine wetlands, as defined by the NWI dataset, and were checked to ensure that the burned signal was visually observable in the training image. Unburned points were randomly selected outside the burn perimeters to represent all unburned conditions (e.g., unburned wetland, forest, urban, grassland, and agriculture). In each of the four regions, we selected 400 burned training points and 400 unburned training points spread across six images representing seasonal variability (two December-February images, two March-May images, one June-August image, and one September-November image). The one exception was South Carolina, where an extra spring image was added to increase the number of images with burned points present. A total of 3200 training points across 25 images were generated. All training images, as well as the number of burned points per image, are shown in Appendix A Table A1. The variable values for the pixel corresponding to each training point were extracted and used to train the burned area algorithm.

Variables Considered
Spectral indices that have been shown to be helpful in identifying burned area were considered as potential variables and are listed in Table 1 [34]. Because Sentinel-2 lacks thermal bands, indices that included thermal bands were excluded from consideration. Single spectral bands (i.e., blue, red, and NIR) were also excluded from consideration, because we found that random forest models fit to the training data with single bands tended to show higher rates of error when applied to the imagery time series. For each of the 14 spectral indices, we considered the (1) single scene variable, (2) z-score value of the variable, and  were calculated, and change from the monthly mean was calculated by subtracting the specific monthly mean over the time series from the pixel value. To account for potential spectral confusion between increases in wetland water level and burned areas, we also included two Standardized Precipitation Index (SPI) [61] variables, where precipitation was aggregated for the last 30 days (spi30d) and 90 days (spi90d), respectively, relative to the date of the image acquisition. The data were derived from 4 km daily Gridded Surface Meteorological (GRIDMET) data [62]. The precipitation data were standardized relative to 1981-2016, assuming an inverse-normal distribution. The 42 spectral variables combined with the two precipitation variables produced a total of 44 variables that were considered for inclusion in the model. The values of each of the variables for all training points were extracted and exported from Google Earth Engine.

Model Variable Selection and Application
We used the Python package, geemap to access, analyze, and visualize the Sentinel-2 image collection, as well as the GRIDMET data, in the Jupyter environment [77]. Random forest classifiers were fit using the Scikit-learn Python module [78]. We sequentially tested the suite of potential predictors and, in each step, selected the predictor that increased the Random Forest's Area Under the Curve (AUC) the most. In each step, predictors were removed if they had high correlation (≥0.95) with any of the selected predictors. Predictors were selected until the increase in AUC of additional variables was <0.0001. We also evaluated how model accuracy changed with the number of training points to ensure our sample size was sufficient (100 iterations, 20% randomly selected each time as validation set). The selected predictor variables were then used to train a final Random Forest model, using Scikit-learn in a Jupyter Notebook environment (https://jupyter.org/, accessed on 12 May 2021), accessed through Miniconda 3 (https://docs.conda.io/en/ latest/miniconda.html, accessed on 12 May 2021). We ran 300 trees with bootstrapping. The trees were translated into text and compiled into a list of trees readable to Google Earth Engine, using the ee.Classifier.decisionTreeEnsemble method, and applied to all Sentinel-2 images. Variable importance was calculated with Python Scikit-learn as the permutation importance, which reports the mean decrease in accuracy.
Temporary turbidity in open water was a source of commission error and was masked out by identifying permanent open water. Long-term (2016-2019) open water was defined as having a median Normalized Difference Wetness Index (NDWI, normalized difference of green and NIR bands) value of >0.14, where water values tended to be positive and land values tended to be negative. This threshold was selected by iteratively testing threshold performance across the study area to mask out open water bodies but avoid masking out wetlands. During the study period, this region experienced two category 5 hurricanes: Hurricane Irma made landfall on 10 September 2017 and traveled up the western coast of Hurricanes of this magnitude cause substantial flooding and damage; therefore, Sentinel-2 images for the Florida panhandle were excluded between October and December 2018, and Sentinel-2 images within the Everglades and the western coast of Florida were excluded between October and December 2017 to reduce post-hurricane commission error. Both hurricanes showed extended impacts, however, into the following spring, primarily along riparian corridors in the path of the hurricane. To reduce riparian commission error along the path of the hurricane, a pixel was masked out if the normalized difference of the NIR and coastal bands averaged between −0.05 and 0.25 across December, January, and February in the winter following each hurricane. This value range and months maximized distinguishing documented burned area from hurricane-related impacts. In addition, very limited manual editing was performed, but only after the validation analysis was complete. In this step, the direct paths of Hurricane Irma and Hurricane Michael were reviewed for false positives in the spring following each hurricane.
For each year, the classified Sentinel-2 images were consolidated into the burn count for January-May and the burn count for June-December. Because both wildfires and prescribed fires were most common in the January-April period, burn counts of 2 or more were reclassified as burned for this 6-month period. In contrast, burn counts of 3 or more were reclassified as burned for the June-December burn count raster in each year. This enabled us to limit the inclusion of TOA images with a lower than average signal to noise ratio. To reduce commission error, only burned polygons >1 ha were retained. Finally, the spring and fall burned areas were combined to produce an annual raster that was filtered by using a 5 × 5 bit error adaptive filter to reduce noise while preserving edges [79].

Validation
We evaluated the Sentinel-2 burned area classification by using validation points (burned and unburned) derived from WorldView-2 and WorldView-3 images (Figure 2a). Worldview-2 (2 m resolution) and Worldview-3 (1.4 m resolution) images (n = 27) were selected by using disproportionate sampling [80,81] to ensure that enough wetland burned area was considered in the validation (Appendix A Table A2). The image acquisition dates ranged from 26 February 2017 to 6 December 2019. The images were acquired via the NextView License as Level 2A 8-band images. Within each high-resolution image, burned area that occurred within or overlapped wetlands, as defined by the NWI dataset, was manually delineated as (1) visually distinguishable from the high-resolution imagery and (2) verified as a burned area by a secondary source, specifically MODIS active fire points [20] or historical high-resolution imagery available in Google Earth Pro, for burns missed by the MODIS active fire dataset. A total of 319.8 km 2 of burned area across 54 fires, within or overlapping wetlands was delineated, which represented 5.3% of the non-overlapping high-resolution image extent. We then randomly generated validation points with the class proportionate to the amount of burned area cover (5%), so that 1000 points were randomly selected within the burned perimeters and 20,000 points were randomly selected outside the burned perimeters. Each point was visually checked to ensure the accuracy of the assigned class and points were reassigned as needed. The creation of the validation dataset was completed independent of the Landsat BA and the Sentinel-2 burned area produced in this analysis. Points used to train the algorithm did not overlap with points used to validate outputs. For the validation points associated with each high-resolution image, the nearest date, relatively cloud-free Sentinel-2 image was processed to burned/unburned and used to evaluate the Sentinel-2 burned area (Appendix A Table A2). For comparison, we also validated the L8 BA, using the same set of validation points (Appendix A Table A2). Because of variable timing in the fire activity relative to the image collection, the Sentinel-2 and L8 BA images selected for validation were not necessarily the closest date match, but instead, the single image that best mapped the burned area. This decision was made because the nearest date image sometimes had substantial cloud cover. In addition, this approach helped avoid cases in which a burned area was visible in a high-resolution image, but the fire event had not yet occurred in a Sentinel-2 or L8 image collected several days prior. NA values within the selected images were excluded from inclusion in the accuracy statistics.
Omission (OE) and commission (CE) errors were calculated for the "burned" category [82]. The Dice coefficient was calculated as the conditional probability, presented as a percentage, that if one classifier (product or reference data) identifies a pixel as burned, the other one will as well, integrating omission and commission errors [83,84]. The relative bias was also calculated, providing the proportion, presented as a percentage, that burned area is under or overestimated relative to the total burned area of the reference product [41]. To account for the influence of validation points being clustered within each high-resolution image, standard errors were calculated by using the stratified combined ratio estimator [85,86]. Using the numerator and denominator of the accuracy metric for the set of validation points associated with each high-resolution image individually, this approach calculated the variance and standard error of each regional accuracy metric (survey R package) [87].

Analysis
The potential contribution of the Sentinel-2 burned area was evaluated by comparing the annual Sentinel-2 burn area raster to wetland burn perimeters compiled from multiple fire databases (Figure 2b). All burn perimeters were compiled from datasets where they were verified manually. Sources of burn perimeters included (1) the Monitoring Trends in Burn Severity (MTBS) [88], (2) GeoMAC, (3) WFDSS, and (4) perimeter datasets including prescribed fires compiled from the US Forest Service, National Park Service, the Florida Forest Service, and the Nature Conservancy. The burned perimeters were required to overlap an NWI polygon, and burn perimeters < 1 ha were excluded to match the minimum unit size of the Sentinel-2 burned area. A total of 555 wetland fire perimeters were included across the study area, totaling almost 537,000 ha over the four-year period. The Sentinel-2 burned area, the L8 BA, and the combined burned area, all compiled as annual composites, were compared to the wetland burned polygons and evaluated for detection count, defined as the number of burn perimeters that showed an overlap between the burn perimeter and burned area, as well as the total area of the burn perimeters classified as burned. Detection by wetland type and by burned area perimeter size were also evaluated. Wetland type was defined by the majority NWI wetland type within the burned perimeter [56]. For this analysis, wetland types were restricted to (1) freshwater (palustrine) forested and shrub wetland, freshwater (palustrine) emergent wetland, and estuarine and marine wetlands, while riverine, lake, and pond NWI wetland types were excluded.

Sentinel-2 Clear-Sky Observation Contribution
Sentinel-2 can contribute to monitoring burned area by increasing the frequency of clear-sky observations across the Southeast. Adding Sentinel-2 observations to Landsat-8 observations in the peak fire period, March/April, increased the observation count between 43% in 2016 and 313% in 2018 across the study area ( Figure 3). On an annual time-step, adding Sentinel-2 to Landsat-8 in 2016 increased observation count from 17 to 46. With the launch of Sentinel-2B in 2017, the observation count increased from 16 with Landsat-8 alone to 78 by 2019 ( Table 2). The increase in Sentinel-2 observation count in 2018 and 2019 relative to 2016 and 2017 was evident not only in observation count ( Table 2) but also the net contribution of burned area extent, relative to the L8 burned area extent. In Figure 4, for example, we observed a greater "new burned area" identified by Sentinel-2 in 2018 and 2019 across the Everglades, relative to 2016 and 2017.  Percentiles represent increase in the cloudfree observations when including Sentinel-2, relative to using Landsat-8, alone. Observation count for March and April was selected to represent the regional peak in wetland fire frequency. Avg = average. Percentiles represent increase in the cloud-free observations when including Sentinel-2, relative to using Landsat-8, alone. Observation count for March and April was selected to represent the regional peak in wetland fire frequency. Avg = average.

Model Performance and Burned Area Validation
When the model score was related to the number of training points, by 1400 training points, we saw limited further improvements in model score, suggesting that the number of training points used in the development of the burned area random forest model was sufficient (Figure 5a). Out of the 44 potential independent variables tested, nine variables were selected for inclusion in the Sentinel-2 Random Forest model (Figure 5b). The zscore of the NIR-red ratio vegetation index (VI43) and the change from monthly mean of the shortwave infrared and NIR ratio (VI45) showed both high variable importance and the biggest contributions to the AUC. The change from monthly mean of the Enhanced Vegetation Index (EVI) also showed high variable importance but lower contribution to the AUC (Figure 5b). In total, one z-score variable, three change from monthly mean variables, three single scene variables, and both the 30-and 90-day SPI variables were all selected for inclusion in the Random Forest model ( Figure 5). We tested the Spearman correlation values between the selected variables and found a maximum correlation of −0.85 between the change from monthly mean EVI and change from monthly mean NDWI. The average absolute correlation value between all of the selected variables was 0.37. The model's AUC was 0.998, while the model's cross-validation score was 0.92 ( Figure 5).

Model Performance and Burned Area Validation
When the model score was related to the number of training points, by 1400 training points, we saw limited further improvements in model score, suggesting that the number of training points used in the development of the burned area random forest model was sufficient (Figure 5a). Out of the 44 potential independent variables tested, nine variables were selected for inclusion in the Sentinel-2 Random Forest model (Figure 5b). The z-score of the NIR-red ratio vegetation index (VI43) and the change from monthly mean of the shortwave infrared and NIR ratio (VI45) showed both high variable importance and the biggest contributions to the AUC. The change from monthly mean of the Enhanced Vegetation Index (EVI) also showed high variable importance but lower contribution to the AUC (Figure 5b). In total, one z-score variable, three change from monthly mean variables, three single scene variables, and both the 30-and 90-day SPI variables were all selected for inclusion in the Random Forest model ( Figure 5). We tested the Spearman correlation values between the selected variables and found a maximum correlation of −0.85 between the change from monthly mean EVI and change from monthly mean NDWI. The average absolute correlation value between all of the selected variables was 0.37. The model's AUC was 0.998, while the model's cross-validation score was 0.92 ( Figure 5). Strong model fit metrics, however, do not necessarily indicate how well the model will perform when applied to a time series of imagery; additional validation is necessary. The absolute average date gap between Sentinel-2 and the high-resolution validation images was 7.0 days, while the absolute average date gap between the Landsat-8 images and the selected high-resolution images was 11.6 days (Appendix Table A2). The L8 BA showed a 47 ± 5% and 8 ± 3% error of omission and commission for burned area, respectively (Table 3). Sentinel-2 burned area, in comparison, showed a lower error of omission (29 ± 7%) but a higher error of commission (30 ± 4%). However, the Dice coefficient (70 ± 5% relative to 67 ± 5%) and relative bias (2 ± 8% relative to −42 ± 4%) were both stronger in the Sentinel-2 burned area validation relative to the L8 BA validation (Table 3). Strong model fit metrics, however, do not necessarily indicate how well the model will perform when applied to a time series of imagery; additional validation is necessary. The absolute average date gap between Sentinel-2 and the high-resolution validation images was 7.0 days, while the absolute average date gap between the Landsat-8 images and the selected high-resolution images was 11.6 days (Appendix A Table A2). The L8 BA showed a 47 ± 5% and 8 ± 3% error of omission and commission for burned area, respectively (Table 3). Sentinel-2 burned area, in comparison, showed a lower error of omission (29 ± 7%) but a higher error of commission (30 ± 4%). However, the Dice coefficient (70 ± 5% relative to 67 ± 5%) and relative bias (2 ± 8% relative to −42 ± 4%) were both stronger in the Sentinel-2 burned area validation relative to the L8 BA validation (Table 3).

Burned Area Extent and Analysis of Burn Perimeter Detection
Burned area impacts a small proportion of land area annually; therefore, burned area products are often difficult to view at a regional scale. To improve visibility, the proportion of a moving window mapped as burned was generated to indicate the burned area patterns ( Figure 6). Across the study area, the Sentinel-2 burned area showed a similar geographical burn pattern between 2016 and 2019 as the L8 BA Product ( Figure 6). However, some differences were visible, for instance, the Sentinel-2 BA showed a higher proportion of burned area in Southwestern and Southeastern Georgia, while the Landsat-8 BA showed a higher proportion of burned area at the Southern Georgia border and Eastern Florida. Across the study area, averages of 2.1% and 1.1% of the NWI wetlands, defined as estuarine and marine wetlands, as well as freshwater emergent, shrub, and forested wetlands, were mapped as burned each year by Sentinel-2 and the L8 BA, respectively. In comparison, The Landsat Burned Area product mapped 0.4% of the conterminous US as burned each year, with 0.6% of CONUS mapped as burned in 2012, the highest fire year [34].
We compared the Sentinel-2 and L8 burned area to the wetland burn perimeter dataset (n = 555) to quantify the relative contribution of Sentinel-2 burned area. Sentinel-2 exceeded Landsat-8 in all four years in the percent of the wetland burn perimeters detected and in all years except 2016 for the percent of the wetland burn perimeters mapped as burned (Table 4). An example of this finding is shown in Figure 7, where Sentinel-2 identified additional burned area relative to the L8 BA, but showed a similar burn pattern as MODIS active fire points and the L8 BA across the Apalachicola National Forest in the Florida panhandle. Across the four-year period, by count the L8 BA detected 60% of the perimeters, while Sentinel-2 detected 79% of the perimeters and together 79% of the perimeters were detected (Table 4). Although burn perimeters are typically an overestimate of burned area and include unburned areas within the burn perimeter, mapping more area as burned suggests the detection of burned area across a larger range of burn severities. By area, the L8 BA detected 32% of the perimeter area compared to 39% detected by Sentinel-2, and together the two sensors mapped 50% of the perimeter area as burned (Table 4). BA showed a higher proportion of burned area at the Southern Georgia border and Eastern Florida. Across the study area, averages of 2.1% and 1.1% of the NWI wetlands, defined as estuarine and marine wetlands, as well as freshwater emergent, shrub, and forested wetlands, were mapped as burned each year by Sentinel-2 and the L8 BA, respectively. In comparison, The Landsat Burned Area product mapped 0.4% of the conterminous US as burned each year, with 0.6% of CONUS mapped as burned in 2012, the highest fire year [34].   The detection of wetland burn perimeters by Sentinel-2 varied by wetland type and burn perimeter size. Burn perimeters overlapping palustrine forested or scrub-shrub wetlands represented 52% of the perimeters, while the remaining were split between palustrine emergent wetlands (26%) and estuarine and marine wetlands (22%). The Sentinel-2 burned area detected 97% of the palustrine forested or scrub-shrub wetland burn perimeters, compared to 74% of the palustrine emergent wetland burn perimeters, and only 37% of the estuarine and marine wetland burn perimeters (Table 5). An example of burned-area detection in South Carolina across an area dominated by forest/scrub-shrub wetlands shows that Sentinel-2 contributes burned area relative to not only L8 BA but also Landsat-7 BA (Figure 8). In addition to patterns by wetland type, we also observed a pattern by burn perimeter size. While only 40% of the perimeters between 1 and 10 ha were detected, 68% of the burn perimeters 11 to 20 ha were detected, and >97% of the burn perimeters 200-500 ha in size were detected (Table 5). Sentinel-2 burn polygons <1 ha in size were filtered out in postprocessing to reduce speckling, which reduced noise but likely increased omission of smaller (1-10 ha) burn perimeters.  (Table 4). Although burn perimeters are typically an overestimate of burned area and include unburned areas within the burn perimeter, mapping more area as burned suggests the detection of burned area across a larger range of burn severities. By area, the L8 BA detected 32% of the perimeter area compared to 39% detected by Sentinel-2, and together the two sensors mapped 50% of the perimeter area as burned (Table 4).    Although the algorithm was trained using data points from four distinct w regions, the burned area algorithm was applied across the entire study area. F shows examples of how the Sentinel-2 burned area algorithm performed, relative Landsat-8 burned area and the MODIS active fire points, outside of the areas in w was trained, including mapping burned area in wetlands along Florida's easter (Figure 9a-c), and mapping burned area in forest and scrub/shrub land-cover typ ure 9d-i). Although the algorithm was trained using data points from four distinct wetland regions, the burned area algorithm was applied across the entire study area. Figure 9 shows examples of how the Sentinel-2 burned area algorithm performed, relative to the Landsat-8 burned area and the MODIS active fire points, outside of the areas in which it was trained, including mapping burned area in wetlands along Florida's eastern coast (Figure 9a-c), and mapping burned area in forest and scrub/shrub land-cover types (Figure 9d-i).

Discussion
Fires in wetland ecosystems are common across the Southeastern US, as well as globally [11,12,52]; however, existing burned-area products underestimate wetland burned area [24,34], which contributes to error in accurately tracking changes in terrestrial carbon storage and carbon emissions [8,18]. In this effort, we used a machine learning algorithm,

Discussion
Fires in wetland ecosystems are common across the Southeastern US, as well as globally [11,12,52]; however, existing burned-area products underestimate wetland burned area [24,34], which contributes to error in accurately tracking changes in terrestrial carbon storage and carbon emissions [8,18]. In this effort, we used a machine learning algorithm, trained on wetland fires, to map burned area extent (2016-2019), using the Sentinel-2 archive across Florida, as well as parts of Alabama, Georgia, and South Carolina. Burned area in wetlands can go undetected when a burned area is visible but omitted by a burned area algorithm due to low burn severity, burnable wetlands being masked out, or an atypical burn signal [11,38]. Burned area in wetlands can also go undetected when the burned area was never visible during image collection dates because of rapid vegetation recovery and cloud cover [25,26]. The finer spatial resolution and greater frequency of the Sentinel-2 imagery collection can help improve estimates of carbon emissions [89] from wetland fires. With the imminent decommissioning of Landsat-7, and Landsat-9 not yet launched or operational, Sentinel-2 or the harmonized Landsat and Sentinel-2 dataset [90] can potentially be used at a national scale to enhance the USGS BA Product. By training the algorithm using burned area in wetlands only and designing the open water mask to avoid masking burnable wetlands to the extent possible, we demonstrated that Sentinel-2 can be effectively used to support regional-scale efforts to find and track wetland burned area. Additionally, modifying the algorithm to include wetland burned area in other burnable wetland types, such as peat bogs, could improve the model transferability to other regions where wetland fires are common.
Assessment of algorithm accuracy is a critical component of burned area algorithms. The accuracy of our Sentinel-2 burned area (29% OE, 30% CE) was similar to other efforts that used Sentinel-2 to track burned area, for instance, a 33% OE and 31% CE from mapping fires with Sentinel-2 across Tropical and Subtropical African savannah sites [48], and a 27% OE and 19% CE mapping fires from Sentinel-2 across Sub-Saharan Africa [42]. However, in converting from single-image classification, on which our validation was performed, to annual composites of burned area, we observed that a small percentage of images in the time series introduced a disproportionate amount of error. To remove potentially problematic images, we used a seasonally specific burn count threshold. This step, while reducing commission error in the annual mosaic, also contradicted our original argument that observation count is key to wetland fire detection [11,12,25]. An alternative approach is to threshold and segment burn area from predicted burn probability surfaces instead of using burn count [27,34]. Using burn probability instead of burn count potentially retains burned areas that are only visible in a single image. However, when testing this approach in the algorithm development stage, we found that this approach also missed documented burns. This source of error highlights the challenge of minimizing commission error in algorithms applied to dense stacks of imagery, while also minimizing omission error when the burn signal is often subtle in wetland ecosystems. While additional preprocessing, such as conversion to surface reflectance, could potentially reduce noise, other efforts have shown that surface reflectance does not necessarily outperform TOA imagery in wetland and water environments [91,92]. Alternative approaches to identify and remove images with a lower-than-average signal-to-noise ratio may therefore be necessary [93].
Relative to other land-cover types, wetlands show both a high diversity in vegetation type and structure, and a high amount of temporal variability, attributes that make it especially challenging to distinguish burned area extent. This challenge is made more difficult because water and burned area show substantial overlap in their spectral reflectance [21,30,31]. Forested wetlands are abundant across the study area but can be difficult to discern, as they occur under the forest canopy. Water levels in this wetland type respond not only to variability in precipitation, but also typically show an inverse relationship with leaf-cover and leaf-out attributable to tree evapotranspiration [94]. Tidal, emergent wetlands are typically vegetated with salt-tolerant species that change with height above sea level, and experience constant tidally induced changes in water level [1]. Freshwater emergent wetlands, in contrast, primarily respond to variability in precipitation, but wetland vegetation types, such as cattails, bulrushes, saw grass, and maidencane, differ in height and structure from other types of wetland vegetation, such as bladderworts and water lilies, which can float on the water's surface [1]. These differences, which are attributable to wetland type, may be why we found that the Sentinel-2 burned area performed unevenly across different wetland types. For example, while most palustrine forested or scrub-shrub wetland burn perimeters were detected, less than 40% of the estuarine or marine wetland burn perimeters were detected (Table 5). This finding suggests that algorithms focused on wetland areas need to be sensitive to post-burn changes in spectral reflectance that are not typical for other vegetation types (e.g., forests). For example, wetland vegetation can burn off during a fire event, leaving open water, and a fire event can also initiate a rapid burst of vegetation growth leaving the area temporarily greener than nearby unburned areas [32,35]. In the algorithm development stage, we found that these types of post-fire transitions did not represent most wetland fires. However, consideration of non-traditional changes (e.g., to open water or greener) may need to be considered in future efforts to map a more complete record of wetland burned areas. Doing so requires spatially coincident information about land-cover and vegetation types, as well as expected timing in water level and phenology.
Algorithms also need to be sensitive to lower burn severity, which is often associated with prescribed burns [21,28,29]. Explicitly training algorithms on both prescribed fires and wildfires can help improve an algorithm's performance, which we saw, by the Sentinel-2 burned area identifying more area of the perimeter dataset as burned, relative to the Landsat-8 BA. However, mapping prescribed fires as distinct from wildfires will require relating remotely sensed BA products to prescribed fire records and datasets [51], which remains challenging [27].
Wetlands are dynamic systems that experience structural changes from seasonal shifts in water level, episodic fire events, and, in the Southeastern US, tropical storms and hurricanes. In this effort hurricane related flooding and vegetation damage was a substantial source of error. For instance, we observed an increase in commission error in wetlands along the southern and western edge of the Everglades starting in September 2017 with the timing of the error coinciding with the landfall of Hurricane Irma on 10 September 2017 that caused substantial flooding and sediment movement as it traveled up the western coast of the Everglades [95]. When we tracked wetlands misclassified as burned, we found a shift in the spectral signature that persisted from September 2017 until February 2018. In contrast, because the area where Hurricane Michael made landfall was dominated by forest, the region experienced substantial blow-down [96], meaning that the damage persisted for a longer period post-hurricane relative to Hurricane Irma. In both cases, we observed an increase in commission error that was most prominent along riparian corridors. The hurricanes complicated efforts to accurately map burned areas near where the hurricane made land fall; however, we saw a rapid reduction in commission error as the hurricanes moved farther inland. In response to this source of error, we removed imagery at the site where each hurricane made landfall for the October-December post-hurricane period and implemented a temporary riparian mask for December-February to further reduce hurricane-related error. However, sustained impacts from the hurricane remained a source of error for the burned area mapping effort. Specifically training the algorithm on storm damage and impacts may help reduce this type of error.
Increasingly, cubesats and other commercial satellites are collecting multispectral imagery that may have potential applications for mapping burned area [97,98]. However, often these satellites are restricted to visible and NIR bands, which enable algorithms to detect changes in greenness [63,70], but prevent the calculation of some of the most common burn related indices, such as the Normalized Burn Ratio (NBR), which requires shortwave infrared bands [66,67]. In wetlands, using satellites that provide shortwave infrared bands may be especially important. As water levels rise in vegetated wetlands, the spectral reflectance in the NIR declines [38], meaning that saturated wetlands can look very similar to burned areas. Even NBR, however, tends to perform less consistently in unforested areas [73], where the index may respond more to soil wetness than plant coverage [12]. Using Landsat satellites, indices that incorporated the thermal band were found to reduce confusion [34,38]. For sensors that lack the thermal band, such as Sentinel-2, we demon- strated that ancillary sources of data, such as GRIDMET standardized precipitation index can potentially help the algorithm identify when darker-than-normal ground conditions may be attributed to an increase in wetness instead of a burn event. Using precipitation data can also potentially support efforts to identify burned areas where vegetation burned off, revealing open water underneath. However, precipitation datasets, such as GRIDMET, are typically coarse in resolution (4 km), while one of the advantages to using Sentinel-2 is its moderate spatial resolution (20 m) relative to existing global burned area products (≥500 m). Alternatively, it is possible that using a synthetic aperture radar (SAR) satellite, such as Sentinel-1, to estimate changes to water level, or a Sentinel-2 or Landsat-based estimate of inundation extent, could serve a similar role in helping distinguish an increase in wetness from a burn event or when a change from vegetated to non-vegetated water may be attributable to a burn event instead of a rise in water level. As Sentinel-1 now provides regularly collected imagery, opportunities exist to develop burned area algorithms from SAR sensors [48,99]. Although the algorithms, to date, have tended to underperform relative to Sentinel-2 [48], SAR sensors, such as Sentinel-1, have the added benefit of not being sensitive to cloud cover or smoke and therefore would provide many more viable images for burned area mapping across the frequently cloudy Southeast, and can potentially help support efforts to monitor burned area and track disturbance-related changes in wetland biomass [100].

Conclusions
This analysis represented the largest-scale effort to date, of which we are aware, to specifically find and map wetland burned area. We demonstrated that Sentinel-2 can be effectively used to map the wetland burned area, particularly from 2018 onward, as the high count of clear-sky Sentinel-2 observations provide an improved chance of detecting burned areas prior to wetland vegetation recovery. The algorithm performed particularly well in forested and scrub/shrub wetlands. However, greater attention to nontraditional burn trajectories (e.g., conversion to open water and enhanced greenness relative to unburned areas) will be essential to further improve the detection of wetland burned area, particularly in tidally influenced estuarine and marine wetlands, where the water level changes hourly, not just seasonally. While using ancillary datasets, such as precipitation, helped distinguish wet conditions from burn conditions, the challenge of distinguishing these two spectrally similar conditions remains, particularly in tidal wetlands. Wetlands, globally, store a disproportionate amount of soil carbon. This carbon storage is vulnerable to loss during wetland fire events, yet fires in wetlands remain understudied and under-mapped. Algorithms, similar to the one developed here, can be adapted for diverse wetland types across the globe and employed to provide more complete burned-area datasets within wetland ecosystems. An improved understanding of wetland fire occurrence, frequency, and extent is critical for the effective management of wetland fuel loads, carbon storage, and wildlife habitat.   Disclaimer: Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Table A1. Sentinel-2 images selected for training and the number of burned points randomly selected across each image. In each image 67 points were selected to represent unburned, except for the South Carolina training images, where 57 unburned points were selected per image to account for the inclusion of an extra training image. FL = Florida, GA = Georgia, N. = north, S. = south.

Region
Sentinel