How Much of a Pixel Needs to Burn to Be Detected by Satellites? A Spectral Modeling Experiment Based on Ecosystem Data from Yellowstone National Park, USA

Abstract

We present a simple modeling technique based on linear spectral mixture analysis to assess satellite detectability of sub-pixel burned area. Pixel observations are modeled using a linear combination of pure land covers, called endmembers. We executed an experiment using spectral data from Yellowstone National Park, USA. Using endmember samples from spectral libraries, pixel samples were assessed on burn detectability using the widely used differenced Normalized Burn Ratio (dNBR). While individual samples yielded differing results for Landsat 8, Sentinel-2, and the Moderate Resolution Imaging Spectroradiometer (MODIS), the average park-wide detectability of burned area was consistent across satellites. For the commonly used dNBR threshold of 0.15, the results indicated that detectability is reached when around a quarter of a pixel’s area is burned. However, a significant percentage of the modeled burned pixels remained undetectable, especially those with low pre-fire vegetation cover. This has consequences for burned area estimates, as smaller fires in sparsely vegetated terrain may remain undetected in moderate resolution burned area products.

1. Introduction

Wildfires have severe impacts on the natural environment and human lives. They accelerate changes in ecosystems by removing vegetation [1,2], and, as such, can exacerbate soil erosion [3,4]. Furthermore, the release of large amounts of greenhouse gases from fires can contribute to climate warming when ecosystems do not fully recover [5]. In addition, air pollution from the fire’s smoke can be deadly [6,7], and the economic costs of both prevention and damage repair are substantial [8]. For these reasons, monitoring fires is extremely important.

Satellite remote sensing is the only practical means of gathering information on fires over large areas. Two remote sensing techniques are often used to detect fires. Firstly, active fires can be detected by the thermal anomalies from their heat release [9]. Secondly, burned areas can be detected based on surface reflectance changes after the fire [10,11]. Fires remove vegetation, produce charcoal, and can expose soil. Active fires can only be detected during the fire, and detection accuracy thus depends on orbital characteristics of the satellites. In contrast, surface reflectance can be mapped after the fire from cloud-free images. This method is thus more suitable to map the extent of fires [12].

There are, however, also limits to the detection capabilities based on changes in surface reflectance. Image pixels, no matter how small, are always a composite of the constituents within them, and this sub-pixel heterogeneity can impact the analyses [13]. If only a small area within a pixel is burned, this may not be detectable by satellites, which can result in high omission errors [14].

An often-used metric for burned area detection (e.g., [15,16,17,18,19]), the Normalized Burn Ratio (NBR), may also suffer from this limitation. This metric takes advantage of spectral differences between reflectance of healthy vegetation and burned area in the near-infrared (NIR, ∼0.75–0.90 $\mathsf{\mu }$m) and shortwave infrared (SWIR, ∼2.09–2.35 $\mathsf{\mu }$m) domains of the electromagnetic spectrum [20,21]:

$$NBR=\frac{NIR-SWIR}{NIR+SWIR}$$

(1)

Healthy vegetation typically has high reflectance values in the NIR domain, but low values for SWIR. Reflectances for charcoal, soil, and bare rock are typically opposite of that, giving low NIR and high SWIR values in deforested or burned areas. As a result, a significant decrease in NBR is an indication of vegetation loss.

However, if only a small fraction of a pixel has burned, the NBR may not decrease much, leading to the burned area remaining undetectable. A large unknown here is how large the fraction of burned area needs to be before detection is possible. This results in some underestimation of fire impacts, as smaller contributions remain undetected. We aim to improve our understanding of this detection threshold by modeling the response of satellites to increasing burned fractions. To achieve this, the study area of Yellowstone National Park, USA, was selected to model spectral responses of burned area.

2. Materials and Methods

In this project, the linear version of spectral mixture analysis (SMA) was used to calculate satellite spectral response to burned areas. SMA is a widely used technique for modeling the composite nature of pixels [22]. For studies on wildfires, it has been used to estimate fire severity (e.g., [23,24]), intensity (e.g., [25]), and vegetation recovery (e.g., [26]). The technique is relatively simple, has been extensively tested, and has proven to be reliable [22]. The approach in SMA is to model the spectral signatures of a pixel as a mixture of pure land covers, called endmembers. A burned pixel may, for example, be a combination of the endmembers grass, bare soil, and charcoal. Typically, SMA is used to estimate abundances of such endmembers in pixels imaged by satellites. Here, we did not use satellite imagery; instead, we created artificial pixel samples from combining these endmembers and used these to determine satellite response. To our knowledge, no prior study has modeled satellite detectability of burned area in such a way.

Firstly, endmember samples were selected for the study area and obtained from spectral databases (see Section 2.2). Next, the samples were grouped and combined to model sub-pixel environments located within the park. A fire in these environments was subsequently modeled by introducing the spectral contribution of a charcoal sample. Finally, we assessed burned area detectability of our modeled spectral mixtures.

2.1. Study Area

The study area for which endmembers were selected is Yellowstone National Park (shortened here to Yellowstone), a nature preserve in Wyoming, USA [27]. The park covers 8903 km², and is a landscape characterized by high volcanic plateaus eroded by glaciers and rivers and flanked by mountains. Lakes cover around 5% of the park’s area. Elevation ranges between 1637 m and 3512 m, averaging at 2479 m [28]. The area is volcanically active and contains three large calderas [29]. The park was selected as a study area because the available vegetation spectral data are complete and easily accessible. Additionally, wildfires occur often in the park during the summer season. The burned environments for this region were modeled using a linear combination of three endmembers: vegetation, substrate, and charcoal.

The dominant vegetation community in Yellowstone is coniferous forest, which represents around 85% of the surface area. The United States National Park Service classifies them under four communities. Of them, forests dominated by lodgepole pine (Pinus contorta subsp. latifolia) are by far the most common. In addition, there are forests dominated by spruce fir (Picea engelmannii and Abies lasiocarpa), whitebark pine (Pinus albicaulis), and douglas fir (Pseudotsuga menziesii). This last vegetation community is interesting as it has the highest frequency of fire [27]. The vegetation patterns in the region vary greatly with elevation and topography, with douglas fir being located at lower elevations, lodgepole pine between 2000 and 2400 m, and the other communities reaching to the upper tree line [30]. Other vegetation communities in the park are grasslands, meadows, sagebrushes, and hydrothermal areas. These are classified under non-forested vegetation. Finally, scattered throughout the park are smaller communities of aspen forests, wetlands, and streamside vegetation [27].

2.2. Endmember Selection

The reflectance spectra of the environments in Yellowstone were modeled using three categories of endmembers: vegetation, substrate, and charcoal. Vegetation spectral data were obtained from the United States Geological Survey’s (USGS) spectral library, which includes data collected from laboratory, field, and airborne spectrometers. Covering wavelengths between 0.2 and 200 $\mathsf{\mu }$m, the library contains measurements of vegetation, minerals, chemical compounds, and manmade materials. Among this library, averaged spectral reflectance data of the dominant vegetation communities in Yellowstone National Park are available. These are averaged top-of-canopy measurements of vegetation communities, previously collected by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) [31]. The used spectra span wavelengths between 370 and 2500 nm at an interval of 10 nm. The vegetation samples that were used in this study are summarized in

Table 1. Overview of vegetation samples used in the analysis. More detailed descriptions of the samples can be found in the USGS spectral library [31].

Type	Vegetation Community	Dominant Species	No. of Samples
Forest	Douglas fir	Pseudotsuga menziesii	3
	Lodgepole pine	Pinus contorta subsp. latifolia	11
	Spruce fir	Picea engelmannii, Abies lasiocarpa	5
	Whitebark pine	Pinus albicaulis	2
Nonforest	Bacterial mat	Chloroflexus aurantiacus, Synechococcus lividus	1
	Conifer–meadow mix		1
	Grass	Festuca idahoensis	5
	Sagebrush	Artemisia tridentata	4
	Sedge		1
	Willow–sedge mix		1

. It is common in spectral mixture analysis to subdivide the vegetation endmember into two: green and non-photosynthetic parts of the vegetation cover. However, as the vegetation samples used are averages of plant communities, they represent spectra that contain both green and non-photosynthetic parts of the plants. Thus, the vegetation endmember samples from AVIRIS represent the vegetation spectra as they would be measured by a satellite instrument.

The endmembers that were used to simulate the substrates in Yellowstone were obtained from the NASA Jet Propulsion Laboratory’s ECOSTRESS spectral library. This library is a collection of laboratory measurements covering a wavelength range between 350 and 15,400 nm of (green and non-photosynthetic) vegetation, rocks, soils, minerals, and some manmade materials [32]. From this library, both rock and soil samples were obtained to simulate the substrate conditions of the park. Rock spectral samples were selected according to the geology of Yellowstone, as described by Keefer [33]. At the wavelengths used in this work (up to 2350 nm), the spectral samples from this library had a measurement interval of 1 nm.

The inclusion of soil samples is important, as some soil-forming processes influence the spectral signature. For example, increases in organic matter result in lower reflection in the NIR domain [34]. Furthermore, in the presence of charcoal, the Normalized Burn Ratio is highly sensitive to soil type [35]. However, selecting soils is less straightforward compared to rocks, as the characteristics of soils are highly variable. Apart from the lithology in which the soil is formed, it is highly dependent on the local climate characteristics [36]. Soils in Yellowstone are typically in the frigid or cryic temperature regime, although warmer regimes form around hydrothermal areas. The climate of the park is characterized as moderately dry [28], and mostly falls within the udic moisture regime [36]. The most common soil types in the park are mollisols (thick, organic soils) and inceptisols (weakly developed mineral soils); however, entisols, andisols, and alfisols are present as well. For a full overview of the soils in Yellowstone, see Rodman et al. [28].

We only selected soils for the more common lithologies, as these have a higher impact on the park-wide results. Initially, the selected soils matched parent materials with the rocks in the park, could form in the temperature and moisture regimes of Yellowstone [36], and are similar to the soils described by Rodman et al. [28]. The rock and soil samples were grouped in accordance with the geological units, as shown in

Table 2. Overview of substrate samples used in the analysis. More detailed descriptions of the samples can be found in the ECOSTRESS library [32].

Geological Unit	Rock/Soil Name	No. of Samples
Precambrian Gneiss and Schist	Gneiss	19
	Schist	26
	Inceptisol dystrochrept	1
	Inceptisol haplumbrept	1
Paleozoic Formations	Dolomite	3
	Limestone	33
	Sandstone	21
	Shale	25
Mesozoic Formations	Limestone	33
	Sandstone	21
	Shale	25
Tertiary Formations	Conglomerate	3
	Sandstone	21
Diorite Intrusions	Diorite	2
	Granodiorite	4
Absaroka Volcanic Breccias	Andesite	6
	Basalt	35
	Mafic tuff	1
	Inceptisol haplumbrept	1
	Mollisol cryoboroll	1
Yellowstone Tuffs	Felsic tuff	8
	Alfisol fragiboralf	1
	Alfisol haplustalf	2
	Inceptisol xerumbrept	1
Plateau Rhyolite	Rhyolite	24
	Inceptisol cryumbrept	1
	Inceptisol plaggept	1
	Mollisol cryoboroll	1
Basalt Flows	Basalt	35
Quartenary Deposits	Travertine	2
	Alfisol fragiboralf	1
	Alfisol haplustalf	2
	Alfisol paleustalf	1
	Inceptisol cryumbrept	1
	Inceptisol haplumbrept	1
	Inceptisol plaggept	1
	Inceptisol xerumbrept	1
	Mollisol cryoboroll	1
	Mollisol haplustall	1

.

However, some very common parent materials were greatly underrepresented in the soil samples due to a lack of soils that fit all the criteria. In these cases, some soils were chosen that are technically outside of this climatic regime, but fit well with the soil descriptions and parent material. Most notably, there was no soil sample found for volcanic ash that also fit the climatic regimes. Thus, the otherwise well-fitting “inceptisol xerumbrept” (which forms in deserts) was added.

The charcoal endmembers were collected by Veraverbeke et al. [37] from the 2011 Canyon fire scar in Kern County, California. Their spectra were measured at NASA’s Jet Propulsion Laboratory, in the same fashion as the used rock and soil samples from their ECOSTRESS library. Three charcoal endmember samples were available and used.

2.3. Sample Preparation

The first step of the analysis was to calculate the near-infrared (NIR) and shortwave infrared (SWIR) response of the endmember samples (Figure 1). To perform this, we followed the procedure outlined by Barsi et al. [38]. For every wavelength of an endmember sample, the response of a satellite band was calculated:

$${\rho }_{ems,b,\lambda }={\rho }_{ems,\lambda }·{\beta }_{b,\lambda }$$

(2)

where ${\rho }_{ems,b,\lambda }$ is the spectral reflectance of an endmember sample $ems$ as observed by a satellite instrument band b at wavelength $\lambda$. ${\rho }_{ems,\lambda }$ is the spectral reflectance of the endmember sample from library (again at wavelength $\lambda$), and ${\beta }_{b,\lambda }$ is the spectral response of a satellite at that wavelength. Both vary between 0 and 1. The satellites for which detectability was assessed are Landsat 8, Sentinel-2, and the Moderate Resolution Imaging Spectroradiometer (MODIS) (

Table 3. Details of the satellite instrument bands used in the analysis. Used are the near-infrared (NIR) and shortwave-infrared (SWIR) bands for Landsat 8, Sentinel-2 (two satellites; A and B), and the Moderate Resolution Imaging Spectroradiometer (MODIS, two satellites; Terra and Aqua).

		NIR Band			SWIR Band
Satellite Instrument	Band No.	$\mathbf{\lambda }$ Range ($\mathsf{\mu }$m)	Resolution (m)	Band No.	$\mathbf{\lambda }$ Range ($\mathsf{\mu }$m)	Resolution (m)	Data Source
Landsat 8	5	0.85–0.88	30	7	2.11–2.29	30	[39]
Sentinel-2 A & B	8	0.78–0.89	10	12	2.01–2.37	20	[40]
MODIS Terra & Aqua	2	0.84–0.89	250	7	2.11–2.16	500	[41]

). Their response functions are shown in Figure 2. The spectral response functions were linearly interpolated to match the spectral resolution of the endmember sample.

The responses to these reflectances were then averaged over the band’s wavelength range:

$${\rho }_{ems,b}=\frac{{\sum }_{\lambda ={\lambda }_0}^{{\lambda }_1}{\rho }_{ems,b,\lambda }·∆\lambda }{{\sum }_{\lambda ={\lambda }_0}^{{\lambda }_1}{\beta }_{b,\lambda }·∆\lambda }$$

(3)

where ${\rho }_{ems,b}$ is the spectral response for the endmember sample ems of a given instrument band b, ${\lambda }_0$ and ${\lambda }_1$ are the wavelength boundaries of the instrument band, and $∆\lambda$ is the wavelength distance between the measurements of the endmember sample. In this case, multiplying the endmember reflectance with the satellite response functions automatically sets the wavelength range (as the responses are set to 0 outside of the band of the instrument). If $∆\lambda$ is constant over the wavelength domain, Equation (3) can be simplified:

$${\rho }_{ems,b}=\frac{\sum {\rho }_{ems,b,\lambda }}{\sum {\beta }_{b,\lambda }}=\frac{\sum {\rho }_{ems,\lambda }·{\beta }_{b,\lambda }}{\sum {\beta }_{b,\lambda }}$$

(4)

For every endmember sample, Equation (4) was used to calculate the NIR and SWIR band values for the different satellite sensors.

2.4. Modeling Burned Area Response

In linear spectral mixture analysis, spectral endmembers are combined linearly to simulate the spectrum of a pixel [42]:

$${\rho }_p=\sum _{ems=0}^n{\rho }_{ems}·f_{ems}$$

(5)

Here, ${\rho }_p$ is the reflectance of a pixel, which is a weighted average of the reflectances of the endmembers used to simulate the pixel ($0\le {\rho }_{ems}\le 1$), with weights given by their spectral contribution fractions ($0\le f_{ems}\le 1$). In linear SMA, these fractions can also be interpreted to represent the surface area covers of the endmembers, and therefore are non-negative and have to add up to unity.

In this work, pixel spectra were modeled using three endmembers: vegetation, substrate, and charcoal. This yields the following calculation for the pixel reflectances:

$${\rho }_p={\rho }_v·f_v+{\rho }_g·f_g+{\rho }_c·f_c$$

(6)

In this equation, ${\rho }_p$, ${\rho }_v$, ${\rho }_g$, and ${\rho }_c$ are the reflectance values of the pixel and its vegetation (v), substrate (g), and charcoal (c) endmembers, respectively. $f_v$, $f_g$, and $f_c$ are the respective cover contribution fractions. The letter g (from ground) is used for substrate endmembers to disambiguate from the s used for starting values in later equations.

The endmember reflectances in this equation (and further equations) are those obtained by Equation (4). Thus, the obtained pixel reflectance is for a given instrument band and is already scaled with the band’s response function. This means that the NIR and SWIR values obtained from this equation are the reflectances of the pixel as would be obtained from the satellites.

Before a fire, a pixel only contains vegetation and substrate endmembers ($f_c=0$); thus, Equation (6) yields

$${\rho }_{p,s}={\rho }_v·f_{v,s}+{\rho }_g·(1-f_{v,s})$$

(7)

${\rho }_{p,s}$ is the pixel’s (p) pre-fire (starting, s) reflectance. $f_{v,s}$ is the vegetation cover of the pixel sample before the fire. This value was varied between 0.05 and 1 in steps of 0.05, yielding 20 pixel samples per ground-vegetation endmember combination. The substrate endmember takes up the rest of the area of the pixel; therefore, its contribution fraction is equal to $1-f_{v,s}$. From this pre-fire pixel reflectance, the Normalized Burn Ratio before fire ($NB{R}_{p,s}$) could be calculated:

$$NB{R}_{p,s}=\frac{{\rho }_{p,NIR,s}-{\rho }_{p,SWIR,s}}{{\rho }_{p,NIR,s}+{\rho }_{p,SWIR,s}}$$

(8)

The cover contributions for the pixels were then altered by replacing parts of the ground and vegetation with charcoal, simulating a burn in the pixel. The model calculates a decrease in the contribution of vegetation ($f_v$) using a variable fraction of burned vegetation ($f_b$). The fraction of burned vegetation ($0\le f_b\le 1$) gives the fraction of the original vegetation ($f_{v,s}$) that is currently burned. The fraction of vegetation ($f_v$) could then be calculated using

$$f_v=f_{v,s}·(1-f_b)$$

(9)

Along with changes in $f_v$, the contribution of ground ($f_g$) and charcoal ($f_c$) also need to change, as the total fraction of the pixel needs to be 1. The parameter $∆c$ (or $∆$char) gives the change in charcoal cover fraction per unit of vegetation lost (burned). The remaining contributions could thus be calculated using

$$f_c=f_b·f_{v,s}·∆c$$

(10)

$$f_g=1-(f_v+f_c)$$

(11)

$∆c=0$ represents vegetation loss without charcoal input. $∆c=1$ represents a fire in which all vegetation is replaced by charcoal; the contribution of ground stays constant. $0<∆c<1$ represents an increase of both charcoal and bare ground. $∆c>1$ indicates that part of the ground will be covered by charcoal in addition to vegetation loss. The parameter can not be negative, as this would result in an illogical mechanism where forest cover increases as it is burned. In the model, results were calculated for $∆c$’s between 0 and 1 in steps of 0.25.

Vegetation was removed until the vegetation contribution became 0, or the burned area became detectable. The detectability was assessed using the differenced Normalized Burn Ratio (dNBR), which is calculated by taking the difference between the pre- and post-fire NBR. After an increase in the burned vegetation fraction, the pixel’s NBR ($NB{R}_p$) and, subsequently, dNBR were calculated:

$$NB{R}_p=\frac{{\rho }_{p,NIR}-{\rho }_{p,SWIR}}{{\rho }_{p,NIR}+{\rho }_{p,SWIR}}$$

(12)

$$dNBR=NB{R}_{p,s}-NB{R}_p$$

(13)

If the dNBR was greater than a certain threshold, the fire was set to be detectable and the cover fractions of the different endmembers $f_v$, $f_g$, and $f_c$ were saved. A burned pixel was set to be undetectable if the burned fraction reached 1 before the dNBR threshold was exceeded. In the model code, multiple threshold values were used to be able to assess the influence of this parameter. The dNBR thresholds varied between 0.05 and 0.25 in steps of 0.05.

The cover fractions were calculated for every substrate sample combined with every vegetation and charcoal sample. However, we were only interested in the range of these values for a certain combination, not necessarily the outputs for the individual library samples. For example, the combination of gneiss (19 samples) with lodgepole pine forest (11 samples) and charcoal (3 samples) has a total of $19\times 11\times 3=627$ model output values. However, the interest is only in the detectability of the gneiss–lodgepole pine combination. Therefore, of these 627 results, only the minimal, mean, and maximal values were saved. This reduced the number of results for this combination to three per variable.

From the burned fraction $f_b$, the endmember contributions can be calculated. The goal of the model was then to find the value of $f_b$ for which the burn could be detected. This can be achieved using a root-finding algorithm; for example, iteratively, by increasing stepwise $f_b$, testing the resulting dNBR and checking if it is above the threshold. Then $f_b$ could be tweaked to attain closer approximations of the threshold. However, it is possible to directly calculate the burned fraction needed for detection. The methodology is explained in Appendix A. The direct calculation is faster and more precise, allowing more samples to be processed within a reasonable amount of time.

2.5. Data Aggregation

With a high number of endmember sample combinations, aggregation of the results was required to draw conclusions. The results were aggregated at different levels (

Table 4. Result aggregation level details.

Level	Name	Results Count	Notes
0	Endmembers	49,572,000	Results are not saved at this level.
1	Groupings	1,170,000	Contains min, mean, and max values of sample combination groupings.
2	Geologies	300,000	Substrate groupings aggregated to their geological units.
3	Park-wide	6000	Weighted average using abundances of geology–vegetation combinations.

), allowing conclusions to be drawn at different scales.

The direct result outputs are defined at level 0. At this level, the results are given per endmember sample combination; each one is calculated using a combination of one substrate, vegetation, and charcoal spectral sample.

As we are not interested in the variation between the samples of a particular endmember (e.g., between basalt samples), we can reduce the amount of data. We perform this by only saving some statistical properties of a collection of these sample combinations. Here, we save only the mean, minimal, and maximal values of the model results in a collection or grouping, essentially giving one value with uncertainty. This yields results at aggregation level 1. In this case, spectral samples are grouped according to their soil name, rock name (see Table 2), or vegetation community (see Table 1). Every grouping is then a combination of a named substrate with a vegetation community. For example, one of the grouping combinations is basalt with lodgepole pine, containing all model results that are a combination of basalt and lodgepole pine.

Higher aggregation levels than the groupings level require some information on the area covers of the substrate and vegetation combinations. However, the surface area covers of rocks or soils within a geological unit are largely not documented. For aggregation level 2, this is accounted for as follows. We assign every rock and soil to a geological unit (per Table 2). For the mean value of a unit, we assume the rock and soil types to be equally abundant on the surface. The high uncertainty in this assumption is taken into account by obtaining the minimal and maximal value of combinations assigned to a geological unit. These extreme values thus assume that the entire unit is dominated by the most extreme spectral samples within them, providing a constraint on the variability.

The final aggregation level, aggregation level 3, assesses burned area detectability at a park-wide scale. To achieve this, the results for the geological units were weighed with their surface area. To find these weights, a geological map [43] and vegetation habitat map [44] of Yellowstone National Park were reclassified and combined to show the geological unit–vegetation community combinations. Area sizes of these combinations were subsequently set as the weights in the result analysis. The classes and corresponding weights are shown in

Table 5. Estimated area cover % of the lithology–vegetation combinations.

	Douglas Fir	Lodgepole Pine	Spruce Fir	Whitebark Pine	Nonforest	Total
Precambrian Formations	0.675	0.607	0.007	0.130	0.517	1.94
Paleozoic Formations	0.141	0.653	0.038	0.344	0.509	1.69
Mesozoic Formations	0.106	0.919	0.234	0.749	0.694	2.70
Tertiary Formations	0	0.002	0	0.001	0.001	0.003
Diorite Intrusions	0.081	0.079	0.005	0.148	0.236	0.549
Absaroka Volcanics	1.18	5.74	0.842	7.63	5.78	21.2
Yellowstone Tuffs	0.826	10.8	1.26	1.09	0.519	14.5
Plateau Rhyolite	0.073	18.8	1.74	2.14	1.75	24.5
Basalt Flows	0.230	1.39	0.045	0.039	0.145	1.86
Quartenary Deposits	1.11	13.4	6.30	1.97	8.25	31.0
total	4.42	52.5	10.5	14.2	18.4	100

.

At aggregation level 3, the model contains three results (min, mean, max) for the entire park. However, parameterization adds additional dimensions; thus, the park-wide results are dependent on a given dNBR threshold (5 options), $∆$char (5), instrument (4), and starting vegetation fractions $f_{v,s}$ (10), totaling 6000 result values.

3. Results

The outputs of the model vary depending on the parameter settings used. Thus, we analyze the detectability by assessing the influences of the parameters. For most of the results, the starting fraction of vegetation ($f_{v,s}$) is set as the independent variable, as detection characteristics are highly dependent on it.

3.1. Differences between Satellite Sensors

First, we identified the differences introduced by the various satellite instruments. We use the aggregation level 2 for this purpose, as it still contains all the variability of the lower aggregation levels while being attributable to certain areas of the park. The deviation in burned fraction ($f_b$) on detection is shown in Figure 3, illustrating how much vegetation loss is required before detection.

A few observations are of note here. Firstly, deviations are larger with increases in starting vegetation fractions. With higher starting fractions of vegetation, the spectral changes as a result of a burn can be larger. Secondly, the mean values of burned fraction (shown in boxplots) vary only very slightly, deviating up to 0.1 from the average of the satellites. This is in contrast with the extreme sample results (range shown in gray), which can deviate up to 0.5 from the average. We can conclude from this that the satellite sensor influence is relatively unimportant when looking at averaged values of sample combinations, but can vary wildly with individual sample results.

It is important to note that the satellites show a characteristic deviation from the mean. Typically, the Sentinel-2 satellites required higher values of burned fractions before detection, essentially being slightly less sensitive to spectral impacts of fire. The opposite is true for the MODIS instrument, which also had the highest deviations from the other instruments. Finally, Landsat 8 had the least skewed, as well as the lowest absolute, deviation from the mean. However, the sensitivity to spectral changes does not directly translate to burned area sensitivity. The large differences in pixel sizes between the instruments is, for example, another important factor.

For subsequent results, the differences between the satellite instruments are small when compared to the variation introduced by other parameters. Therefore, the data in this section are shown for the Landsat 8 satellite, as it has the most balanced detection characteristics. The results for Sentinel-2A, Sentinel-2B, and MODIS are shown in Appendix B.

3.2. dNBR Threshold Influence

Detectability of a sample pixel was dependent on the detection threshold. When the dNBR threshold increases, a stronger burn spectral signal is required for detection to be reached. As a result, higher fractions of burned area were required, and the number of pixels that were undetectable increased (Figure 4).

In that sense, lower thresholds would be more desirable. However, this also increases the detection of non-burn-related changes in the spectral signal. For example, removal of vegetation changes the spectral signal in a similar way to a burn, as bare soil or rock typically has a lower NBR. We modeled this by setting $∆$char to zero, where the removal of vegetation is compensated only by increases in substrate contribution ($f_g$). The resulting detection for this model run is shown in Figure 5.

Comparing the results of Figure 4 with Figure 5, the conclusion can be taken that on average for the park, a vegetation removal scenario required higher levels of vegetation loss compared to a pure burn scenario. In addition, the percentage of undetectable pixels was higher. However, there seems to be no objectively optimal setting for the dNBR threshold.

3.3. Park-Wide Detectability

Assessing the average number of undetectable pixels as a function of dNBR threshold and $∆$char parameters (Figure 6), some important results become apparent. We can expect a burn’s spectral impact to be somewhere between a pure burn scenario and a vegetation removal scenario; thus, the variability in detectability introduced by the $∆$char parameter indicates a possible range of results. Here, it is notable that a significant fraction of the modeled pixels were undetectable. At the common dNBR threshold of 0.15, we find that even in the best case scenario (a pure burn scenario: $∆$char = 1), more than a quarter of pixels were undetectable.

However, as shown previously, pixels with smaller burnable areas (low starting vegetation fractions) have lower detectability. The results for dNBR = 0.15 are shown in Figure 7. It shows that nearly all pixels with a starting vegetation fraction smaller than 0.2 were undetectable. For the best-case (pure burn) scenario, even at starting vegetation fractions of 0.4, 20% of burned pixels remained undetectable.

Figure 8 shows the required vegetation cover area loss, before either detection or otherwise complete burn of the vegetation. These were calculated by multiplying fraction burned with the starting vegetation fractions, or alternatively by summing charcoal contribution and increases in substrate contribution. The mean values show that even pixels that started completely covered in vegetation ($f_{v,s}=1$) required, on average, a vegetation loss at around 20% of the area of the pixel before detectability was reached. Pixels with a balanced starting contribution of vegetation ($f_{v,s}\approx 0.5$) resulted in the highest vegetation losses before detection, at 30 to 45% of the pixel’s area.

4. Discussion

There are notable uncertainties in our analysis, as denoted by the large range of possibilities in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8. This range resulted from outliers within vegetation–geology groupings. These outliers represent burned pixels in which specific vegetation and geology spectral samples were combined. It is unlikely, yet not entirely impossible, that these combinations are fully representative for real-world combinations. The mean spectral responses of the vegetation–geology groupings showed much less variability, and those likely more realistically represent the real-world variability.

Our results indicate that more than 15% of an average pixel needs to burn before detection is reached, even in pixels with high vegetation covers. This is also illustrated by Randerson et al. [45], who estimated increases in burned area with 25 to 54% when including small fires that were not detected in the MODIS burned area product. The detection results are irrespective of pixel sizes; the estimated fraction of the pixel required to burn does not change with increasing resolution. However, its implication for burned area estimates is dependent on spatial resolution. With smaller pixel sizes, the burned area omission errors should decrease. Small burns cover larger fractions of smaller pixels, and, as such, have increased detectability. This means that burned area detection will be easier when using Landsat 8 and Sentinel-2 compared to using MODIS, which has much larger pixel sizes. For example, Ramo et al. [14] and Glushkov et al. [46] found significant increases in burned area estimates when reassessing MODIS burned area products with the higher-resolution Sentinel-2 imagery. However, the results shown in this work indicate that some omission errors are expected to persist, even for these smaller pixel sizes.

Our results shown in Figure 8 indicate that the highest omission errors are to be expected in sparsely vegetated terrain, as such pixels require nearly complete vegetation loss before detection. This effect may be further exacerbated, as sparsely vegetated areas may be more conducive to smaller fires because of fuel limitations [47]. Additionally, in regions that show high short-term spectral changes unrelated to fire, the small burned areas may be more difficult to discern from non-fire-related spectral shifts. This includes desert spring environments, where, as a result, small burned area detection from satellites is especially limited [48].

As our analysis uses linear spectral mixture analysis, it assumes that the light reflected by the surface interacts with only one component of the pixel and that there is no spectral influence of potential neighboring pixels. While this is a common assumption in spectral mixture analysis, it is also a large simplification of reality, especially in vegetated areas where multiple and nonlinear scattering are common [49].

There are other factors in burned areas that may further influence detectability but were not considered in our analysis. The vegetation community spectral samples were collected from AVIRIS and therefore represent the vegetation at a specific time of the year. In this case, the measurements were taken on 6 August, near the peak of the growing season [31]. The assumption made in this work is that the pre-fire and post-fire satellite measurements occur on the anniversary date, and thus phenological shifts are minimal. However, any changes in spectral reflectance due to phenological shifts, moisture content, or shadows may alter the detectability of the burned areas [50].

We also did not explicitly account for different fire severity levels in our study, and this can be a major influence on the detectability of burned areas [51,52]. In our spectral modeling experiment, the fire effects were modeled by partly or fully replacing vegetation cover by soil and charcoal cover. In low to moderate fire severity plots, vegetation is often scorched. We have not accounted for the possibility of scorched vegetation; however, its inclusion in burned pixels instead of vegetation replacement by charcoal would likely even further decrease detectability.

One method for improving detectability of pixels with lower vegetation fractions is proposed by Miller and Thode [50], in which they scale dNBR by the pre-fire NBR. This results in the so-called Relative differenced Normalized Burn Ratio (RdNBR). This method increases sensitivity of the dNBR when the pre-fire NBR is close to 0. This improves detectability of burned areas in sparsely vegetated terrain. As our results indicate a high dependence on the pre-fire NBR of the pixels on their detectability, the RdNBR metric may yield significant detection improvements and should be further evaluated for burned area mapping in heterogeneous and sparsely vegetated areas.

Burned areas that cover small fractions of a pixel are more likely to be found either as a result of small, isolated fires or near the edges of larger fires. One way of limiting commission errors in the latter is by using region-growing algorithms [53]. In this method, pixels at the core of a burned region are first classified using a strict classification threshold. Next, a more relaxed threshold is used in the vicinity of these core burned pixels. This assumes that changes in NBR in vicinity to core burned pixels are likely resulting from the same fire. It allows more leniency towards smaller burned areas near the fringe of a fire, while changes in NBR unrelated to fire are not included. This technique may be used to improve detectability of burn fringes, but is not suitable for detecting small isolated fires.

We used dNBR thresholds to delineate burned area. While the dNBR has often been used as a burned area discriminator (e.g., [15,54]), we acknowledge that spectral reflectance in individual bands and other spectral indices are also often used (e.g., [11,19]). While small differences may exist depending on which spectral discriminator is used for burned area, we think that there would be commonalities with some of the main findings of our work with regards to, for example, the effect of pre-fire vegetation cover, burn fraction, and threshold values on burned area detectability.

The results presented in this work are representative for Yellowstone National Park, the region where the spectral endmembers were obtained from or selected for. While the park contains varied types of vegetation and substrate, Table 5 shows that the dominant lithologies in the park are various volcanic rocks, in total representing around 60% of the surface cover. In addition, more than 80% of the park’s vegetation is represented by coniferous forests. In other regions, especially those with very different landscape characteristics, burned area detectability may yield different results.

5. Conclusions

We present a viable method for assessing satellite response to burned area. Using spectral data of vegetation, substrates, and charcoal and combining them using linear spectral mixture analysis, it was possible to model pre- and post-fire environments in Yellowstone National Park.

Differences between the spectral sensitivity of the Landsat 8, Sentinel-2, and MODIS to burned area were small. The MODIS instrument is slightly more sensitive towards burn spectral characteristics, while the Sentinel-2 satellites are relatively less sensitive. The Landsat 8 sensor showed a balanced response to burned area. These small influences may lead to slightly higher commission errors for MODIS and relatively higher omission errors for Sentinel-2.

Our results indicated that significant fractions of a pixel’s vegetation need to burn before detection is reached. In addition, a significant percentage of the modeled pixels would remain undetectable. At the widely used dNBR threshold of 0.15, the average park-wide results showed that around a quarter of the pixel needs to be burned before detection becomes possible, and more than a quarter of the burned pixels remained undetectable. These results show that detection of burned area using dNBR may be accompanied by substantial omission errors.

Our spectral sensitivity analysis is independent of spatial scale, yet the implications become larger for coarser resolution sensors such as MODIS. Our results thus further support the continued investments in ongoing large-scale burned area mapping efforts at resolutions around 20 to 30 m from Sentinel-2 and Landsat 8. Using instruments with higher spatial resolutions directly lowers omission errors. Such efforts are especially important to quantify burned area from small isolated fires that often burn small amounts of fuels.