Field-Validated Burn-Severity Mapping in North Patagonian Forests

Abstract

Burn severity, which can be reliably estimated by validated spectral indices, is a key element for understanding ecosystem dynamics and informing management strategies. However, in North Patagonian forests, where wildfires are a major disturbance agent, studies aimed at the field validation of spectral indices of burn severity are scarce. The aim of this work was to develop a field validated methodology for burn-severity mapping by studying two large fires that burned in the summer of 2013–2014 in forests of Araucaria araucana and other tree species. We explored the relation between widely used spectral indices and a field burn-severity index, and we evaluated index performance by examining index sensitivity in discriminating burn-severity classes in different vegetation types. For those indices that proved to be suitable, we adjusted the class thresholds and constructed confusion matrices to assess their accuracy. Burn severity maps of the studied fires were generated using the two most accurate methods and were compared to evaluate their level of agreement. Our results confirm that reliable burn severity estimates can be derived from spectral indices for these forests. Two severity indices, the delta normalized burn ratio (dNBR) and delta normalized difference vegetation index (dNDVI), were highly related to the fire-induced changes observed in the field, but the strength of these associations varied across the five different vegetation types defined by tree heights and tree and tall shrub species regeneration strategies. The thresholds proposed in this study for these indices generated classifications with global accuracies of 82% and Kappa indices of 70%. Both the dNBR and dNDVI classification approaches were more accurate in detecting high severity, but to a lesser degree for detecting low severity burns. Moderate severity was poorly classified, with producer and user errors reaching 50%. These constraints, along with detected differences in separability, need to be considered when interpreting burn severity maps generated using these methods.

1. Introduction

Burn severity, defined as the magnitude of ecological change associated with a wildfire [1], is a key element for understanding ecosystem dynamics as it influences processes that shape community structure and composition, nutrient availability, and wildlife habitat [2,3,4]. At broader scales, the spatial distribution of fire effects can determine species replacement and alter landscape dynamics [5,6]. Consequently, the assessment of fire effects is essential for understanding ecosystem resilience and designing management strategies.

Remote sensing allows for the spatial description of burn severity while overcoming cost and logistic limitations of extensive field surveys [7], and numerous mapping methodologies have been proposed and tested (see review by Meng and Zhao [8]). These are based on changes in the spectral signature of fire-affected surfaces, mainly associated with moisture and chlorophyll content of vegetation: red and near-infrared (NIR) reflectance decreases and short-wave-infrared (SWIR) reflectance increases [9]. Several indices combining these bands have been tested for burn severity mapping including the normalized difference vegetation index (NDVI) [10], the normalized burn ratio (NBR) [11], the mid-infrared burned index (MIRBI) [12], and the char soil index (CSI) [13]. As fire effects are dependent on pre-existing vegetation, bi-temporal approaches use the difference between pre- and post-fire indices; however, this requires precautions related to imagery selection to minimize the differences in phenology, solar angle, geometric correction, and sensor calibration [7]. Uni-temporal analyses, based only on post-fire observations, are less affected by these issues but have problems in discriminating burned areas from other spectrally similar areas [14]. Several studies have demonstrated high correlations between field-based reference observations of burn severity and post-fire and differenced versions of spectral indices [15,16,17]. Other protocols use a multitemporal approach to improve the reference data (e.g., the multitemporal dNBR_MT [18] that relies on time series to select reference pixels) or avoid limited image availability due to quality issues (e.g., the relativized burn ratio (RBR) [19], which uses compositing techniques to obtain pre- and post-fire data).

The sensitivity of these methods varies across different vegetation types [20,21] so they need to be field validated to ensure that reliable information is generated. The validation should compare the spectral index to the equivalent surface process or property [1]. Remotely sensed burn severity has proven to be reliable in forests, savannahs, shrublands, and grasslands of different regions [13,22,23]. Some studies have also explored the variation in the sensitivity of indices for detecting changes in different vegetation groups within the same region and have proposed group-specific parameters for burn severity estimation [24,25,26]; grouping criteria have often considered structure, composition, and/or regeneration processes (e.g., contrasting resprouting species and species regenerating mainly from seed). However, field validation still needs to be conducted in many regions with complex vegetation mosaics including North Patagonia.

Fire is a major factor shaping landscapes in North Patagonia and fire-induced vegetation dynamics have been described for the most widespread community types ranging from closed-canopy tall forests of Nothofagus species (either evergreen or deciduous species), dense tall shrublands dominated by resprouting woody species <12 m tall, to grassy open woodlands of the conifer Austrocedrus chilensis [27,28]. Regeneration strategies (e.g., obligate seeding tall trees versus resprouting shrubs) and community flammability have been associated with positive fire feedbacks that maintain a mosaic of pyrophobic forests and pyrophytic shrublands over much of the region [29,30,31]. While forests dominated by fire-sensitive obligate seeders are vulnerable to fire-induced conversion to shrublands, some forests in North Patagonia are dominated by Araucaria araucana, a thick-barked fire-resistant species also capable of resprouting following fire [32,33,34]. Thus, the general model of tall forest conversion to more flammable shrub-dominated vegetation types in northern Patagonia needs to be examined for vegetation types characterized by A. araucana. In addition, fire activity in Araucaria-dominated forests has been shown to be sensitive to both interannual climatic variation and to variability in the frequency of human-set ignitions [35]. Projected shifts in fire regimes caused by climate and land cover changes raise concerns regarding forest resilience in relation to probable increases in the frequency, extent, and severity of fires in North Patagonia [36,37,38] for which improved methods of assessing fire effects related to burn severity are essential.

Araucaria araucana-dominated forests are characterized by a mixed severity fire regime including both high and low severity burns [34,39]. Burn severity varies spatially over relatively short distances, leading to a complex mosaic of differently structured vegetation patches across the landscape. Mixed forests of A. araucana and Nothofagus pumilio or N. antarctica are typically controlled by moderate- to high-severity fires that originate from stands with scattered surviving A. araucana surrounded by small multi-aged groups of A. araucana and a single post-fire cohort of Nothofagus spp. [34]. Post-fire establishment of the obligate seeder N. pumilio decreases under conditions of higher severity [40]; if the accompanying species is N. antarctica, a post-fire cohort is rapidly created by vigorous resprouts that reach 1 m in height within three years [39]. In terms of the understory, higher severity induces decreases in richness and abundance and changes in composition, with exotic species being more abundant in situations of low severity [41]. Regeneration in burned areas is mainly controlled by burn severity and basal area, with influences of topography and biotic interactions, but fire refugia also play a role as seed sources for post-fire dynamics in the longer term [42]. As an important driver of the dynamics of these forests, burn severity is generally assessed in the ecological studies of fire events [43,44,45]. However, field-validation of burn severity mapping methodologies for Patagonian forests is scarce [46]. In the case of Araucaria-dominated forests, the relative version of dNBR (RdNBR) was tested using reference field data gathered 10 years after a fire in Chile [42], but no analogous work is available for the Argentinian forests that are under drier climatic conditions.

The aim of this work was to develop field-validated burn severity maps of two recent (2013–2014) wildfires that affected a complex mosaic of A. araucana-dominated vegetation types in North Patagonia in order to propose a replicable mapping methodology for the region. We address the following questions:

• How do burn-severity estimators derived from satellite images relate to field-based measurements of fire impact in A. araucana vegetation types defined by stand structure and species composition?
• What are the classification thresholds that best discriminate severity classes across the range of A. araucana vegetation types?
• How well do fire-severity maps that are based on the different best performing spectral indices agree?

First, we tested the association between field measurements of burn severity (four years following the fire events) and widely used indices based on remotely-sensed data and conducted a sensitivity analysis to select the best estimators of burn severity across a range of vegetation types. Then, we adjusted the class thresholds for each selected index and assessed their accuracy, contrasting it against other thresholds proposed in previous publications. Finally, we generated burn severity maps for the two studied wildfire events that affected A. araucana forests by applying the most accurate methodologies and evaluated the agreement between them.

2. Materials and Methods

2.1. Study Area

North Patagonian forests grow on the hills and valleys of the Andes in Chile and Argentina, within a steep moisture gradient that decreases from 4000 mm annual precipitation in the west to 500 mm on the eastern side [32,47]. The precipitation is mostly concentrated between May and September and occurs both as rain and snowfall; summer is the dry season. Mean monthly temperatures range from 4 °C to 15 °C, and maximum values reach up to 35 °C in summer. These conditions favor the occurrence of wildfires in summer, which can be more extreme during drought periods. Electric storms are the main natural cause of ignition, but currently, most fires are human-ignited [35,48].

These forests present numerous endemics and are mainly dominated by species of the genus Nothofagus [32,49]. In Argentina, large portions of Patagonian forests are protected under the National Parks Administration. Most of the tree dominant species are relatively sensitive to fire because of their thin to moderately thick bark. Among them, the more fire-sensitive tree species are obligate seeders, whereas the fire-resistant trees and shrubs resprout after fire [50]. All of these species are shade-intolerant and present a regeneration mode that is linked to disturbances. The structure of these forests varies from tall pure or mixed-species stands to <12 m tall stands of multi-stemmed shrubs and small trees including two-tiered stands. Wildfires are a widespread disturbance agent, due partially to the fire-promoting features of some species of the tree canopy and the understory. Consequently, this landscape has been historically shaped by fires, and the mosaic of vegetation types largely reflects the extent and severity of the past events [32,50].

This study focused on two large wildfires that burned between December 2013 and January 2014, near Lakes Ñorquinco and Rucachoroy, in the northernmost area of the Argentinian Patagonian forests (from 39°06′S to 39°13′S and from 71°05′W to 71°27′W; Figure 1). The Ñorquinco fire burned an area of 1626 ha and the Rucachoroy fire extended over 1336 ha, mainly on the valley bottoms and moderate to very steep hillsides with a northern- and southern-aspect, respectively. Both events occurred in summer on days with strong westerly winds within a period of extreme drought, which resulted in widespread intense fires.

2.2. Remote Sensing Data

We used Landsat 8 OLI images downloaded from the USGS repository. As both fires occurred in nearby areas, the same set of images was used in the assessment of both events (Path 232, Row 87; Figure 1). Images with low cloud cover were selected from the Level 1 Landsat Collection, which includes terrain correction, corresponding to the end of the summers that preceded (7 April) and followed the fires up until the year of field survey (23 March 2014; 26 March 2015; 28 March 2016; 15 March 2017; 2 March 2018;

Table 1. Field and remote sensing data: identification and date of wildfires, date of field survey, number of plots, and capture dates of the Landsat 8 OLI images.

Wildfires		Field Survey		Landsat 8 OLI images
ID	Date	No. of Plots	Date	Path–Row	Pre-Fire	Post-Fire
Ñorquinco	29/12/13 to 20/01/14	55	20/02/18 to 02/03/18	232–87	07/04/13	23/03/14; 26/03/15; 28/03/16; 15/03/17; 02/03/18
Rucachoroy		56

). In this way, the first post-fire image corresponds to the end of the growing season when the fires occurred, minimizing the effect of early recovery or invasion, and images of the following years might reflect either of these processes and/or delayed mortality. The selection of close-to-anniversary images allows the comparison of pre- and post-fire images with no concern about the phenology and solar angle effects. All image processing was conducted in QGIS 2.18.20 [51].

Images were converted to surface reflectance with the semi-automatic classification plugin, using dark-object subtraction. Red, near-infrared (NIR), and short-wave infrared (SWIR1 and SWIR2) bands were used to generate indices that are commonly used for burn severity estimation: normalized burn ratio (NBR) [11], normalized difference vegetation index (NDVI) [10], mid-infrared burned index [12] and the char soil index (CSI) [13]. The differenced version of the indices (dNDVI, dNBR, dMIRBI, dCSI) were also included in this work as they measure a magnitude of change, which is closer to the definition of burn severity. These indices have been widely applied in burn severity mapping, with different accuracies in contrasting vegetation types [26,52,53]. To compensate for the effect of the misalignment between the plot and grid centers, index values associated with each mCBI plot were sampled through bilinear interpolation.

2.3. Field-Based Reference Burn Severity Data

Field metrics of burn severity were collected between February and March 2018 in accessible areas of the Ñorquinco and Rucachoroy fires. Prior to the fieldwork, preliminary severity maps were developed using reference dNBR thresholds [11], and were used to identify homogeneous patches. Uniformity is important due to the offset between the circular plot and grid cell centers. As these points are rarely aligned, surveyed information usually combines data registered in more than one pixel. Therefore, this precaution in plot location improves the relationship between the field and remotely sensed data. One hundred and eleven 15 m radius plots were distributed in areas of low, moderate, and high severity across all vegetation types. The number of plots was determined following Chuvieco [54] for the error estimation of classified images with a 95% probability and maximum allowed error of 5%.

Extended assessment, that is sampling the year following the fire, is expected to enable accounting for delayed mortality or survival of vegetation. Due to the timing of the fieldwork (four years after the fires), this analysis constitutes a long-term assessment. One of the dominant species of these forests (A. araucana) is documented to resprout from epicormic and protected terminal buds following fire damage, but many sprouts dry in the following years and only some individuals continue regenerating the crown [39,55]. Consequently, gathering field data a few years after fire events (in this case, in the fourth growing season after the events) should allow for a better characterization of fire impact on these forests. However, the long time lapse from burning to field sampling impedes the evaluation of fire effects on the components of vegetation that recover more quickly (i.e., herbs and low shrubs), limiting the direct application of indices designed for more immediate surveys.

We based our field protocol on the composite burn index (CBI) [11]. In the summer of 2016–2017, we conducted a preliminary assessment where we surveyed 20 plots distributed along the severity gradient and the different vegetation types affected in both fires. Based on this information, we developed a modified CBI index (mCBI) that better fit the specific conditions of this study (

Table 2. Assessment table of the modified composite burn index (mCBI): ranking scale for each factor of the four strata (substrates, tall shrubs, intermediate trees, big trees) corresponding to unburned, low, moderate, and high severity.

Strata	Severity Scale
	Unburned	Low		Moderate		High
	0	0.5	1	1.5	2	2.5	3
Substrates
Litter	intact	charred	absent	absent	absent	absent	absent
Duff	intact	Intact	partially charred	charred	scarce	very scarce	absent
Light Fuels	intact	partially charred	charred/absent	absent	absent	absent	absent
Medium Fuels	intact	Intact	partially charred	partially charred	charred/absent	absent	absent
Heavy Fuels	intact	Intact	intact	intact	partially charred	charred	absent
Tall shrubs (1 to 6 m)
% altered foliage	0%	5%	10%	30%	50%	75%	100%
% resprout	100%	95%	90%	60%	30%	<10%	<1%
Intermediate trees (subcanopy)
% green foliage	100%	90%	80%	60%	40%	<10%	0%
% black foliage	0%	5%	20%	30%	60%	>85%	100%
% charred trunk height	0%	10%	30%	40%	60%	70%	>90%
Big trees (upper canopy)
% green foliage	100%	90%	80%	60%	40%	<10%	0%
% black foliage	0%	5%	20%	30%	60%	>85%	100%
% charred trunk height	0%	10%	30%	40%	60%	70%	>90%

). The modifications included the exclusion of the fast recovering components (herbs and low shrubs) and the elimination of variables that required knowledge of some attributes of the pre-burn vegetation in the plots (e.g., pre-burn presence of small herbs and shrubs and standing dead trees), as a consequence of the lack of pre-fire data and the high structural variability of these forests. Consequently, low severity was characterized by high proportions of green foliage, low char height on trunks, and unburned litter and woody debris; high severity was defined by the opposite conditions; and moderate severity implied an intermediate situation (Figure 2). The proportion of dead fire-killed trees was recorded in every plot.

2.4. Data Analyses

2.4.1. Sensitivity Analysis

We explored the relationship between the field and remotely sensed estimates of burn severity by using scatter plots and correlation analysis. Correlation between mCBI and percentage of dead trees was analyzed to assess the performance of the index. Pearson correlation tests have been used in similar studies [15,56]; we used Spearman tests because the data did not meet the assumption of normality, tested using the ‘fitdistrplus’ function in R. Due to the four-year gap between the fires and the field survey, and considering that both early and delayed resprouting occurs in these forests, we evaluated the most useful acquisition date for severity mapping based on the strength of the association between mCBI and image-derived indices. Indices showing a correlation higher than 0.70 (absolute values) were then subjected to a sensitivity analysis. This was assessed through a separability index [57], calculated as:

$$S_{i-ii}=\frac{{\mu }_i-{\mu }_{ii}}{s{d}_i-s{d}_{ii}}$$

(1)

where S_i-ii is the separability of two consecutive severity classes (low-moderate or moderate-high); μ_i and μ_ii are the mean values of the index for each severity class; and sd_i and sd_ii are the respective standard deviations. Overall separability was obtained as a weighted mean, following Nguyen Tran et al. [26], to assign higher importance to the detection of high severity (0.35 S_low-moderate; 0.75 S_{moderate-high}). To determine constraints on sensitivity related to species-specific reproductive and other traits, this analysis was also conducted separately for each vegetation type (described below) and then summarized as a mean value for each severity index. The higher the value of S, the better the discrimination between two classes; a value of 1 or higher is desired for indices to be used for mapping [58].

Vegetation types were discriminated based on canopy height and regeneration strategy of the dominant species (Figure 3). Tall forests were dominated by trees, with a single closed canopy layer and heights between 30–40 m; multi-stemmed structures of <15 m in height occurred in low forests or two-tiered stands when combined with tall trees. Dominant species were either obligate seeders (S) or resprouters (R). Five vegetation types were identified: Tall/S, dominated by N. pumilio; Tall/R, including stands of N. obliqua or A. araucana; Tall/RS, characterized by mixed stands of N. pumilio-A. araucana; Two-tiered/R, codominated by N. antarctica and A. araucana; and Low/R, dominated by N. antarctica. A. araucana forests varied from open- to closed-canopy stands and, although it can regenerate both by seedlings and resprouts, the crown recovery is the fastest recovery mechanism. The 3 to 6 m tall Chusquea culeou native bamboo occurs in varying quantities in all vegetation types and acts as a highly flammable ladder fuel because it is mostly very fine (grass leaves) or fine (<24 mm diameter) material, and when the culms die following mast events, it leaves massive amounts of fine dead fuels [50].

2.4.2. Thresholding and Accuracy Assessment

To define the class thresholds, we followed two different approaches using mCBI as a reference with class limits of 1.25 (low/moderate) and 2.25 (moderate/high) [59]. First, we fitted linear and polynomial models. Linear and non-linear models have been widely applied for describing relationships between field severity data and spectral indices [17,60,61]. Polynomial models appear to better describe the non-linear relation between these variables [62], but the improvement may be only marginal [17]. Having graphically checked assumptions of normality, homoscedasticity, and independence, best-fit models were selected following the criterion of lowest Akaike score and class thresholds were derived. As an alternative strategy, we determined thresholds by an iterative procedure: increasing values for each threshold (low/moderate and moderate/high) were combined until the highest proportion of samples was correctly classified [42,59].

To select the best performing method for mapping, we compared the accuracies of the classifications derived from the adjusted thresholds and others recommended by Key and Benson [11] and by Assal et al. [42]. The latter study was conducted in Araucaria-dominated forests in Chile, similar to the vegetation in our study area. We applied a bootstrapping method, sampling with reposition, to generate 111 confusion tables and the associated accuracy indices (total accuracy, Kappa, and the producer and user errors for each severity class). Based on this information, we conducted Kruskal–Wallis non-parametric variance analyses to evaluate significance of apparent differences in accuracy.

2.5. Burn Severity Mapping

The two most accurate methods were used to generate burn severity maps for each event. Perimeters of the burned areas were obtained through supervised classification on dNBR and dNDVI layers and manual editing using technical reports (M. Mermoz, pers. com.) and high resolution images as references. Finally, we assessed the agreement between the maps by contrasting the information generated on the proportion of area affected under each severity class.

3. Results

3.1. Index Sensitivity

Field measurements of burn severity were highly and significantly correlated (0.87, p < 0.001;

Table 3. Spearman correlation between the field and remote burn severity estimates for the end of the season when the fires occurred. * Asterisks denote significant correlations (p < 0.001).

	mCBI	%DT	NBR	NDVI	CSI	dNBR	dNDVI
mCBI	1	0.87 *	−0.81 *	−0.85 *	−0.80 *	0.80 *	0.77 *
%DT	0.87 *	1	−0.77 *	−0.83 *	−0.78 *	0.77 *	0.77 *

%DT: Percentage of dead trees.

), indicating that mCBI is a good estimator of post-fire tree mortality. The correlation between mCBI and the tested indices was strongest for data acquired at the end of the summer of 2013–2014 when the fires occurred (>0.70) and declined in the subsequent years for all burn severity indices, except for MIRBI, dMIRBI, and dCSI, which consistently presented low correlations (<0.60; Figure 4). These indices were discarded from further analyses. In the year of the end of the growing season when the fire occurred, the strongest Spearman correlation with mCBI was that of the NDVI (−0.85), followed by NBR, CSI, dNBR, and dNDVI (Table 3). NDVI also had the strongest correlation with %DT (−0.83), whereas the correlations of the other four indices were weaker but similar (0.77 to 0.78).

The highest sensitivity was shown by NDVI and dNBR, but the indices performed differently when vegetation type was considered. In general, separability was around 1 for tall or low stature vegetation types, regardless of the regeneration strategy, but was always under 0.6 for the two-tiered forest (

Table 4. Mean separability of indices calculated for the complete dataset (All), discriminated for each vegetation type: low canopy, dominated by resprouting species (low/R); tall canopy, dominated by resprouting species (tall/R); tall canopy, dominated by seeding species (tall/S); tall canopy, codominated by resprouting and seeding species (tall/R-S); two tiered canopy, dominated by resprouting species (two-tiered/R) and summarized as a general mean (Overall).

		Plots	NBRpost	NDVIpost	CSIpost	dNBR	dNDVI
All		111	0.95	1.06	0.90	1.10	0.93
Canopy	Regen.
Low	R	30	0.92	1.04	0.94	0.94	1.03
Tall	R	25	1.00	1.21	0.92	0.98	1.03
Tall	S	18	1.00	0.73	0.80	1.12	0.89
Tall	RS	27	1.58	1.64	1.40	1.56	1.05
Two tiered	R	11	0.43	0.23	0.29	0.58	0.40
Overall			0.99	0.97	0.87	1.03	0.88

R: resprouter-dominated stand; S: obligate seeder-dominated stand; RS: mixed stand of resprouters and obligate seeders.

). In resprouter-dominated forests, NDVI showed higher separability, whereas NBR was more sensitive for obligate seeders. Severity in stands with both types of species was well detected by all indices (separability > 1). When these partial values were summarized, the most sensitive index was dNBR and the separability of all the others was below 1. dNDVI presented a more balanced sensitivity across vegetation type than NDVI and was particularly more sensitive for the two-tiered stands.

3.2. Thresholds and Accuracy of Selected Indices

Three severity indices were selected for thresholding: dNBR, dNDVI, and CSI. For all indices, the best model was a second order polynomial, with an adjusted R² of 0.77 for dNBR, 0.71 for dNDVI, and 0.70 for CSI (all p < 0.001). The model-derived thresholds were moderately close to those obtained through the iteration process, and for dNBR, the latter approach produced class limits that were the most similar to the field-based reference classes (

Table 5. Accuracy assessment of burn severity classifications. Values in brackets indicate the classification thresholds for low-moderate/moderate-high severity. H and p values correspond to the Kruskal–Wallis variance test for each accuracy index and lowercase letters indicate significantly different groups.

			Low		Moderate		High
Classification	Total Acc.	Kappa	UE	PE	UE	PE	UE	PE
dNBR Ref. 1 (270/600)	77 ± 1	64 ± 2	15 ± 2	11 ± 2	55 ± 2	21 ± 2	1 ± 1	27 ± 2
	B	c	B	c	bc	f	d	a
dNBR Ref. 2 (275/500)	82 ± 1	68 ± 2	15 ±2	11 ±2	46 ± 2	49 ± 2	11 ±1	11 ± 1
	De	d	B	c	d	c	c	e
dNBR It. (300/500)	83 ± 1	70 ± 2	17 ± 2	7 ± 2	41 ± 2	49 ± 2	11 ± 1	11 ± 1
	E	d	Ab	d	e	c	c	e
dNDVI It. (200/300)	83 ± 1	69 ± 2	20±2	3 ±2	33 ±2	59 ±2	13 ±1	9 ± 1
	E	d	A	e	f	b	b	ef
CSI It. (1070/920)	73 ± 1	54 ± 2	12 ± 2	39 ± 2	64 ± 2	44 ± 2	11 ± 1	16 ± 1
	A	a	C	a	a	d	c	c
dNBR Mod. (205/550)	77 ± 1	61 ± 2	6 ± 2	26 ± 2	58 ±2	33 ± 2	11 ± 1	19 ± 2
	B	b	D	b	b	e	c	b
dNDVI Mod. (116/323)	80 ± 1	65 ± 2	6 ± 2	28 ± 2	52 ± 2	32 ± 2	11 ± 1	14 ± 2
	C	c	D	b	c	e	bc	d
CSI Mod. (138/75)	81 ± 1	66 ± 2	17 ± 2	7 ± 2	43 ± 2	67 ± 2	15 ± 2	8 ± 1
	D	c	B	d	de	a	a	f
H	372	310	310	267	353	249	575	539
p	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001

dNBR Ref. 1: thresholds as in Key and Benson (2006); dNBR Ref. 2: thresholds as in Assal et al. (2018). It.: Iterative thresholding; Mod.: Model-derived thresholding. Total Acc.: total accuracy; UE: user’s error; PE: producer’s error.

).

The most accurate classifications were those generated using the iteratively obtained thresholds for the indices dNBR and dNDVI, with total accuracies over 80% and Kappa indices around 70%. Reference class limits taken from the study by Assal et al. [42] in an area in Chile close to our study site produced a similarly accurate categorization (Table 5). These classifications were significantly more accurate than the other tested ones. In terms of class accuracy, the general pattern was shared among all classifications: the lowest errors were found in the high severity class, followed by low severity, while the moderate severity class was the least accurate one (Table 5). For the most accurate protocols, the producer’s error for the high severity class ranged from 9% to 11% and the user’s error varied between 11% and 13%. These errors took values of 3–11% and 15–20% in the case of the low severity class and reached levels of 49–59% and 33–46% for moderate severity, respectively.

3.3. Burn Severity Maps

Burn severity maps were generated using dNDVI and dNBR indices with thresholds obtained by the iteration process (Figure 5). In the case of the Ñorquinco fire, most of the area was burned with high severity (60–70%) while the rest of the area was evenly distributed between moderate and low severity. The agreement between the two maps for this event was 81%, with most of the disagreement occurring in the moderate burn class (42% coincidence). In the case of the Rucachoroy fire, the area was burned in proportions of 50%, 30%, and 20% with high, low and moderate severity, respectively. Again, the two maps coincided in categorizing 81% of the area and the moderate class was the one with the lowest agreement (52% coincidence). Considering dNBR maps as a reference, misclassification was more often due to an underestimation of burn severity (15% of the burned area) rather than an overrating (6%) and mainly occurred in the moderate class (

Table 6. Confusion matrix comparing the severity maps based on dNBR and dNDVI with iteratively defined thresholds including the area burned in the Ñorquinco and Rucachoroy fires. Area expressed in hectares.

		dNBR It.
		Low	Moderate	High	Total (ha)	Agreement
dNDVI It.	Low	658	183	18	859	77%
	Moderate	47	288	267	603	48%
	High	0	136	1703	1839	93%
	Total (ha)	706	608	1988	3302
	Agreement	93%	47%	86%

).

4. Discussion

Our results, based on the collection of field-based reference data, confirm that reliable burn severity estimates can be derived from spectral indices for Araucaria araucana-dominated vegetation types in North Patagonia. We demonstrate that the most commonly employed indices (dNBR and dNDVI) corresponded highly with fire-induced changes observed in the field, although this association was not equally strong across different vegetation types. Time since fire influences this correspondence, with the correlations being higher when spectral indices are obtained from the first post-fire growing season. Based on these indices, we were able to determine class thresholds that produced accurate classifications and identify their constraints as well as finally propose a validated burn severity mapping methodology for the region.

4.1. From Index Suitability to Index Selection

The strong relations of NBR, NDVI, dNBR, dNDVI, and CSI with mCBI suggest the potential use of these spectral indices for burn severity estimation in North Patagonian forests. Previous research in Mediterranean pure and mixed-species forest stands showed that NBR performed better than NDVI for severity estimation and that the association was always weaker in the case of bi-temporal indices [16,52]. A stronger relation to the post-fire uni-temporal indices was also detected in boreal forests and wetlands, where correlation to field severity metrics was higher for NBR than for the other twelve spectral indices including NDVI and CSI [15]. In Australian sclerophyllous forests and woodlands, dNBR outperformed dNDVI in discriminating low severities in a ‘cool’ fire but not in a ‘hot’ fire where classifications where similar [63]. dNDVI was also weaker than dNBR in estimating severity across different vegetation types in the USA, ranging from grasslands and shrublands to conifer and deciduous forests [64,65]. In contrast, CSI was better than NBR for mapping severity in African savannas [13] and MRIBI, CSI, or NDVI were more suitable in wet heath, dry heath, or coniferous stands of European heathlands, respectively [25]. In our study, MIRBI was the only index that showed a weak relation to field severity metrics and was discarded from further analyses.

In all cases, the strength of these relations decreased over time, indicating that burn severity in these forests is better expressed by reflectance in the first growing season following fire events. In Pinus ponderosa dominated forests, the dNDVI correlation with CBI decreased after the second post-fire year, but the dNBR correlation was strong five years later [66]. In North Patagonian forests, the heterogeneity at plot and landscape levels resulted in multiple post-fire scenarios combining various processes. Patches dominated by vigorous resprouters would be expected to mask severity faster than stands of slow resprouters or obligate seeders; a similar situation would occur in open stands where the growth of grasses can be strongly enhanced. In contrast, delayed mortality and slow recovery, as exhibited by A. araucana, would be presumed to be limited or not easily detected in the first period following the fire, leading to an under- or overestimation of severity. Our results suggest that the heterogeneity in species composition and stand structures across this landscape generates a mixture of diverging post-fire reflectance trends, weakening the potential of spectral indices to assess burn severity as time since fire increases. This heterogeneity also generates variability within the growing season, so we compared images at the end of the summer to ensure all species were actively growing. Future research should assess the performance of these indices using a multitemporal approach that could summarize the variability or integrate the information of several post-fire seasons. In the U.S., composite-derived dNBR, RBR, and RdNBR indices outperformed the bi-temporal versions, although improvement in validation statistics was variable in some cases [67].

In our study in North Patagonia, all indices showed an overall suitable sensitivity, but with variable sensitivities across different vegetation types as has been found in studies conducted in other regions. Severity detection by dNBR and dNDVI has been found to perform better in forests than in grasslands and shrublands of temperate and boreal regions in the USA [64], with differences between coniferous- and deciduous-dominated stands [68]. Vegetation class also determined the potential of spectral indices for severity mapping in a heterogeneous landscape in Australia, which was higher for forest covered areas than for grasslands and shrublands [24]. In this same region, index suitability changed along with the regeneration strategy of the dominant species: dNBR was better for resprouter-dominated forests and dNDVI was superior for those composed by obligate seeders [26]. In contrast, our results showed that dNBR sensitivity to burn severity was higher in obligate seeder-dominated stands whereas dNDVI was more sensitive in stands dominated by resprouting species. Differences in tree life forms might explain these opposite results: while the studied Australian forests were all dominated by evergreen angiosperms (Eucalyptus spp.), in our study areas, the dominant trees included deciduous angiosperms and one evergreen conifer, so differences in spectral behavior could be expected. Coinciding with Nguyen Tran et al. [26], CSI sensitivity was lower in all vegetation classes in our study. Regardless of the variation, the separability values were around the expected for indices with good discriminatory power (>1) [58] for nearly all vegetation types.

The effect of vegetation type on indices suitability in our study and in previous similar studies suggests that burn severity estimation could be improved if calibration is conducted separately for each vegetation type. Allen and Sorbel [68] suggest using different scales for coniferous versus deciduous boreal forests. Schepers et al. [25] and Nguyen Tran et al. [26] followed this approach and adjusted predictive models or class thresholds using the best-performing index per vegetation group. This was not feasible in our research because the number of mCBI plots, when stratified by vegetation type, was not sufficient to conduct robust accuracy analyses. However, the investigation of sensitivity differences is relevant for gauging the constraints of the severity estimates.

4.2. From Thresholding to Validated Classifications

In our study, second order polynomials were the best fitting models for all indices, as was previously described for this vegetation type in Chile [42] and other forested regions [16,66,69]. Different models have performed better in other forest types: the sigmoidal model in temperate forests in the USA [70], saturated growth model in Canadian boreal forests [71], and third order polynomial in Russian forests [61], etc. It has also been suggested that using non-linear- instead of linear-models might not refine the severity prediction in spite of the improvement of the fitting parameters [72]. In the cases of dNDVI and dNBR, we found that the model-derived classifications were less accurate than those generated using the iteratively defined thresholds.

All severity estimations in our study performed within the accuracy range reviewed by French et al. [73] with a mean accuracy of 73%, ranging from 59% to 95%. Fire impact has been estimated with similar accuracy in Australian forests, where dNBR performed better than the dNDVI (Kappa index of 77% and 57%, respectively) [63]. In temperate North American forests, dNBR detection of severity was accomplished with accuracy ranging from 62 to 66% [21,70]. The best dNBR and dNDVI classifications presented global accuracies higher than 80% and Kappa indices around 70%, indicating that severity assessment through these methods is reliable. Nevertheless, global indicators should never represent the only evaluating criteria [59]. Unbalanced performance of severity classifications is common, especially on the low and moderate classes [74]. Moderate severity was the least accurately detected class by Soverel et al. [22] in Canadian forests as well as by Kurbanov et al. [61] in Russia. Coinciding with Cocke et al. [21], these studies were best at classifying areas corresponding to high severity. Our results show that, in all cases, high severity was the best discriminated class, with balanced user’s and producer’s errors, and moderate severity was poorly classified (errors around 50%). Low severity class was satisfactorily distinguished, but alternately, one of the two errors was always higher. This is probably due to the high variability of spectral indices in the moderate class, associated with the heterogeneity of fire effects and vegetation types that the field severity indices summarize. This class includes different fire effects on surface and understory modifications that might not be captured by passive sensors [7]. These constraints should be considered when extracting conclusions from the spectral estimation of burn severity.

When subsequent validations are conducted for the same region, the confidence in remote sensing assessments is expected to increase and, ideally, reach a point where the need for field data is minimized [11]. For forests dominated by A. araucana, similar to our study area, the study conducted nearby in Chile by Assal et al. [42] is the only precedent of field validated severity mapping and the RdNBR thresholds they proposed are similar to the dNBR limits we determined. Moreover, the accuracy of the classifications derived from either set of values in our study is higher than the one reported in their work (68% global accuracy and 52% Kappa index), most likely because our study was based on field data that was acquired closer to the date of the fires. Overall, these results indicate that comparable severity assessments can be executed in these forests. The higher accuracy and the different value of our proposed thresholds in comparison to those fixed by Key and Benson [11] sustain the idea of a need to calibrate the methodology to each unique region.

4.3. Application: Burn Severity Mapping

The two most accurate classification protocols (dNBR and dNDVI) were used to generate burn severity maps of the two fire events in the Araucaria araucana-dominated vegetation of North Patagonia. The agreement between these maps reflects the strengths and weaknesses described in the preceding sections: coincidence was high for the high and low severity classes, which were more accurately detected, but declined to 50% for moderately burned areas. Regardless of the spatial distribution, both maps produced the same estimation of proportion of area burned. Similarity in the spatial patterns of NDVI- and NBR-derived severity estimates was detected by Lhermitte et al. [75] in African savannas, who found differences when comparing severity maps generated using the bi-temporal approach and the control pixel method. Comparing maps allows for a demonstration of the impact of methodology in the thematic presentation of burn severity and, thus, can identify the limitations of the generated information. Our results indicate that either approach is similarly reliable for severity mapping in North Patagonia. However, the dissimilar association of NDVI and NBR to different vegetation types could determine the selection of one over the other. We found that greater sensitivity to burn severity was shown by NDVI in areas dominated by resprouting tall shrub species and NBR in obligate seeder stands of tall trees. Further studies specifically comparing the performance of these indices in areas dominated by species of different regeneration strategies would be needed to confirm if one index is preferable, depending on the predominant vegetation type.

5. Conclusions

This work explored the performance of the most frequently used spectral indices for burn severity mapping in North Patagonian forests. Our findings suggest that reliable estimates can be obtained from the dNBR and dNDVI using the adjusted thresholds. As in other forested regions, field-validated maps increase the confidence in these analyses and contribute to a replicable methodology and comparable results. This is particularly important in Patagonian forests where burn severity is usually assessed on a per-fire event basis, with no standardized methods or specific comparison with previously applied thresholds. However, interpretation of the results of our study should consider limitations related to the low-accuracy identification of moderate severity as well as the differential sensitivities of indices for vegetation types of different structure, species composition, and contrasting regeneration strategies of resprouters versus seeders. Similar research including a greater range of vegetation types throughout the Andean-Patagonian region is necessary for generating a burn-severity mapping methodology that performs well across these vegetation types.