Detecting Burn Severity across Mediterranean Forest Types by Coupling Medium-Spatial Resolution Satellite Imagery and Field Data

Abstract

In Mediterranean countries, in the year 2017, extensive surfaces of forests were damaged by wildfires. In the Vesuvius National Park, multiple summer wildfires burned 88% of the Mediterranean forest. This unprecedented event in an environmentally vulnerable area suggests conducting spatial assessment of the mixed-severity fire effects for identifying priority areas and support decision-making in post-fire restoration. The main objective of this study was to compare the ability of the delta Normalized Burn Ratio (dNBR) spectral index obtained from Landsat-8 and Sentinel-2A satellites in retrieving burn severity levels. Burn severity levels experienced by the Mediterranean forest communities were defined by using two quali-quantitative field-based composite burn indices (FBIs), namely the Composite Burn Index (CBI), its geometrically modified version CBI (GeoCBI), and the dNBR derived from the two medium-resolution multispectral remote sensors. The accuracy of the burn severity map produced by using the dNBR thresholds developed by Key and Benson (2006) was first evaluated. We found very low agreement (0.15 < K < 0.21) between the burn severity class obtained from field-based indices (CBI and GeoCBI) and satellite-derived metrics (dNBR) from both Landsat-8 and Sentinel-2A. Therefore, the most appropriate dNBR thresholds were rebuilt by analyzing the relationships between two field-based (CBI and GeoCBI) and dNBR from Landsat-8 and Sentinel-2A. By regressing alternatively FBIs and dNBRs, a slightly stronger relationship between GeoCBI and dNBR metrics obtained from the Sentinel-2A remote sensor (R² = 0.69) was found. The regressed dNBR thresholds showed moderately high classification accuracy (K = 0.77, OA = 83%) for Sentinel-2A, suggesting the appropriateness of dNBR-Sentinel 2A in assessing mixed-severity Mediterranean wildfires. Our results suggest that there is no single set of dNBR thresholds that are appropriate for all burnt biomes, especially for the low levels of burn severity, as biotic factors could affect satellite observations.

1. Introduction

Throughout the Mediterranean basin, the year 2017 was particularly catastrophic due to huge wildfire damage in several Mediterranean countries [1]. After a wildfire occurrence, a number of key ecological processes (e.g., tree mortality and regeneration by seeds and resprouts) were affected to different degrees, depending on the extent and the severity of the fire. Thus, post-fire management options, and the following practices to restore the ecological and hydrological functionality of forest-burnt ecosystems, are closely related to both the extent and severity of fire [2,3,4,5,6]. Knowledge of the spatial pattern of wildfire severity is also crucial to estimate post-fire losses of ecosystem services [7,8]. After the fire, detailed and rapid knowledge of the level of damage and its spatial distribution is the first desirable information [9]. At the landscape scale, this information can be summarized by a map of burn severity, where the term ‘severity’ is used to indicate holistically the magnitude of change, caused by aggregate effects of the fire on the vegetation community, with respect to the pre-fire condition [10]. From an ecological point of view, severity is defined as a measure of the extent to which a pre-disturbance vegetation community has been altered and its ecosystem properties have been disrupted [11,12], denoting above- and below-ground organic matter consumption by the fire [13]. Overall effects of fires are classified into first and second order, which correspond respectively to short- and long-term post-fire effects [14]. However, both effects can be quantified in terms of burn severity [13,15]. For the purposes of this work, “burn severity” defines the degree of fire-induced changes, which in turn fits well with ecological post-fire restoration requirements of burnt forest areas [13,14].

Fire severity is related to the amount of heat released by a fire and the duration of heating [14,16] which determine the proportion of individual plants that survive, and thus affect post-fire regeneration processes. In this sense, several attributes have been used to estimate burn severity from field assessment in the post-fire environment [13]. Many such attributes take into account the effects of fires on plant mortality of the dominant stratum, described in terms of basal area removed [12], plant cover reduction [17], branch consumption [18], and the proportion of branches remaining in the canopy [19]. Within this framework, the composite burn index (CBI) was proposed as a quali-quantitative field measure of burn severity experienced by a forest community [10,20]. The CBI field protocol splits the vertical profile of the forest community into five strata comprising understory and overstory vegetation, from the forest floor to the upper canopy, three of which belong to the understory and two to the overstory. Several post-fire plant community attributes, like the degree of litter and deadwood fuel consumption, percentage of foliage alteration, and char height, are visually estimated and numerically scored per stratum. Recently, a new version of the CBI, called the geometrically modified CBI (GeoCBI), was proposed as an alternative index since it is more consistently retrieved from remotely sensed data [15]. Compared to the CBI, the GeoCBI takes into account two plant cover attributes: the fraction of vegetation cover (FCOV) of each stratum and the degree of change in the leaf area index (LAI), exclusively of intermediate and dominant strata.

From an optical perspective, when vegetation is burned, a drastic reduction in visible to near-infrared surface reflectance (i.e., 0.4–2.5 µm) is observed. The depth of the change in question is associated with the degree of severity on vegetation, and is captured as a variation in the value of spectral indices [21]. This opens up opportunities to use satellite optical sensors with resolutions in the 10–30 m range to evaluate burn severity [22]. Indeed, the normalized burn ratio (NBR) was developed for this purpose [14] and is typically used as normalized delta (dNBR) between pre- and post-fire conditions [10,15]. Although other spectral indexes (e.g., dNDVI and dEVI) have been used in burn severity assessment [23], the dNBR represents a widely used reflectance-based index to detect the levels of the severity of fires from medium-resolution remotely sensed data [23,24,25], and is used to produce maps of the burned areas [13,14]. The Landsat satellite family is one of the most widely used [23,26] for the detection of burn severity. Additionally, in recent years, Sentinel-2A has attracted the interest of the remote-sensing community due to an array of improved characteristics, such as wider swaths, lower geometric distortion, higher spatial and spectral resolution, and open data policy [27,28]. Although most recent studies report promising application, both of Landsat-8 OLI and Sentinel-2A MSI, there is an active debate on which satellite is better at detecting burn severity [28,29]. In July 2017, a rapid and partially assessed map of Vesuvius wildfires was produced by Copernicus Emergency Management Service ([EMSR213] Vesuvio: Grading Map). Nevertheless, the remote sensing of burn severity needs field-based calibration of the reflectance-based index dNBR, whose values may differ markedly among the burnt plant community types [30] as a consequence of species composition, stand structure, and season of wildfire occurrence [9,31].

Local studies calibrating the reflectance-based index and remote sensors are necessary to: (i) identify the better relationship between ground-based and remote-sensed burn severity indices, and (ii) determine the most appropriate dNBR thresholds in order to adequately map spatial variability of burn severity. Therefore, in this context, the main objective of this study was to compare the capability of the dNBR reflective spectral index obtained from Landsat OLI and Sentinel 2A MSI sensors in classifying burn severity quantitatively. The case study concerns the 2017 multiple summer wildfire events, which occurred in a Mediterranean spatially heterogeneous forest cover, in the Vesuvius National Park (southern Italy). Specific objectives were (i) to evaluate the accuracy of dNBR from these two satellites as a predictor of burn severity levels assessed by the two field-evaluated composite burn indices (CBI and GeoCBI), and (ii) to verify if the ranges of burn severity developed by Key and Benson [10] fit to our Mediterranean forest. Therefore, knowing the influence of both vegetation cover (of different strata) in a forest stand, and the influence of the pixel spatial resolution on reflectance pixel values, we hypothesized that the combination of GeoCBI and dNBR derived from satellite Sentinel 2A sensors gives better accuracy in capturing spatial variability of mixed-severity fire effects.

2. Materials and Methods

2.1. Study Area, Climate, and Wildfires of 2017

The Vesuvius National Park covers a surface of ~8000 ha, of which 3800 ha is forested, surrounded by a densely urbanized area (Figure 1). The iconic landmark of Vesuvius, with its asymmetric truncated cone-like structure, is the only stratovolcano still active on mainland Europe [32]. On its steep flanks, the soil is depthless and lava rocks frequently appear on the surface, although at lower elevations the soil of the anthropogenic terraces is about 2 m deep [33]. This substrate corresponds pedologically to a sandy, gravel-rich andosol derived from unconsolidated pyroclastic deposits [34,35].

The climate was Mediterranean-type, and in the summer season (JJA) of 2017 the average air temperature was 33 °C (±3 °C) with a total rainfall of about 3 mm. The maximum air temperature of 40 °C was registered in August. In this scenario, several started wildfires occurred from 2 July to 31 August 2017, and the firefighting service extinguished 149 wildfires [8].

2.2. Forest Type

The land use of the Vesuvius National Park is characterized by volcanic geosites (lava fields and crater), forest, agricultural, and high-density urban areas. Forests are dominated by pure and mixed broadleaved stands (1542 ha; 41%), and pure and mixed Mediterranean conifer stands (1504 ha; 40%), followed by shrublands (729.5 ha; 19%) (Figure 1).

Pure and mixed broadleaved coppice stands are dominated by Castanea sativa, followed by stands of Mediterranean oaks, namely Quercus pubescens, Quercus ilex, and the exotic Robinia pseudoacacia. These forests are distributed mainly on the northern and eastern flanks of the Somma-Vesuvius complex. Multistemmed coppice stands exhibit a high density and a low stem size differentiation (diameter at breast height: 26 ± 9 cm, height: 13 ± 4 m). The understory layer covers 80% of the surface area.

Three light-demanding Mediterranean pines (Pinus pinea, followed by Pinus pinaster and, sporadically, Pinus nigra) were planted at the end of the nineteenth century on the southern and southwestern flanks of the Somma-Vesuvius complex to prevent soil erosion and to promote the establishment of late-successional and light-tolerant native Mediterranean broadleaved species. The even-aged pine stands have a low degree of tree size differentiation (diameter at breast height: 35 ± 5 cm; height: 22 ± 4 m). Moreover, the average canopy base height is 16 (±3) m corresponding to a considerable distance of the fuel canopy from the forest floor. Litter layer is a 5–7 cm thick and non-compact undecomposed needles and twigs. Understory is characterized by a low degree of cover (<10%).

The shrubland communities (alien Genista aetnensis, followed by native Spartium junceum and Cytisus scoparius) are distributed around the cone from the ring plain on the northern sector to the southwestern flank of Mt. Vesuvius. Structurally, these shrublands are pluristratified stands with a high density of stem and stool, and a total height of 5–6 m.

2.3. Field Data Collection

To assess field-based burn severity, 90 circular plots, 16 to 20 m in diameter, were selected one year after the wildfires (summer 2018). They were stratified according to the forest type and a preliminary photointerpretation of burn severity levels. To avoid artificial inflation and increase accuracy of classification, plots were preferably located in fuzzy severity categories [36] in order to emphasize fire effects in the low to moderate severity levels.

Here, post-fire conditions were visually evaluated by means of two field-based burn indices (FBIs): the composite burn index (CBI) [10] and its geometrically modified version (GeoCBI) [15,37]. Geographical coordinates of each plot center were recorded with a hand-held decimeter-level differential global positioning system (Trimble GPS Geoexplorer XT).

The FBI protocols (CBI and GeoCBI) applied in this study are both based on field visual assessment and scoring of the fire effects in five strata; three belonging to the understory layer (A—substrate, B—herbs, shrubs, and trees lower than 1 m, C—shrubs and trees with a height ranging from 1 to 5 m); and two belonging to the overstory layer (D—intermediate trees, and E—dominant and co-dominant trees). Hence, in each sampling plot, all attributes describing fire effects in each stratum were visually estimated and scored with a value between 0 and 3 based on ratings of about 30 biophysical variables [10,15]. If a given stratum level or individual attribute is not present, it is not incorporated into the total field-based composite index score. Choosing which FBI values to use as a threshold between severity classes is somewhat subjective [17]. For the purpose of this study, the severity classes reported in Table S1 were used in order to better identify forest stands that experienced low and moderate effects with fire. Unburnt severity identifies the category in which all the biophysical variables scored a value equal to zero. Nearly complete consumption of the surface vegetation, but minimal scorch or mortality of intermediate and dominant canopy strata, and minimal organic soil consumption, rated as low severity. Scorching or mortality of nearly all of the intermediate and dominant strata or litter consumption over most of the area was required for the high severity class, with moderate-low and moderate-high severity ratings falling between these extremes (see numerical thresholds in Table S1).

2.4. Remote Sensing of Burn Severity

This study focused on two mid-resolution satellites, Landsat-8 and Sentinel-2A. The operational land imager (OLI) sensor on board Landsat 8 and the Multi-Spectral Instrument (MSI) on board Sentinel-2A provide spatial resolution of 30 m and 10–20 m, respectively, with global coverage [22]. Two level 2 surface reflectance Landsat-8 scenes (path/row 189/32) acquired on 25 October 2016 and on 28 October 2017, respectively, were retrieved from the USGS Earth Explorer tool. Likewise, two level 2 surface reflectance Sentinel-2A scenes (tile T33TVF) acquired on 30 October 2016 and 15 October 2017 were retrieved from the European Space Agency (ESA) Copernicus Scientific Data Hub. Level 2 surface reflectance data were used for both platforms, ensuring correction for atmospheric effect, terrain, view zenith angle, and sun angle [38,39].

The pairs of acquisition dates were chosen to provide a best estimate of pre-and post-fire conditions, minimizing differences in vegetation phenology, weather conditions and solar zenith, and elevation angles [40], as allowed by available cloud-free datasets.

Normalized burn ratio (NBR) was computed for each scene as [10]:

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

(1)

where NIR and SWIR correspond to the near infrared and short-wave infrared reflectance bands of Landsat-8 OLI (bands 5 and 7, both 30 m resolution) and Sentinel-2A MSI (bands 8A and 12, both 20 m resolution) respectively. Although there are two Sentinel-2A MSI NIR bands (8 and 8A), the Sentinel-2A NBR was derived using the narrower MSI NIR band 8A (855–875 nm) as it is more similar to the Landsat-8 OLI NIR band 5 (851–879 nm) [41]. The pre-fire and post-fire NBR was then used to calculate the differentiated NBR (dNBR):

$${\mathrm{dNBR}=\mathrm{NBR}}_{\mathrm{pre}}-{\mathrm{NBR}}_{\mathrm{post}},$$

(2)

For each set of field-based composite burn index (CBI and GeoCBI) plot coordinates, dNBR values were extracted using bilinear interpolation [42], such that the resulting dNBR is the average of the four nearest pixels. Vector and raster data were managed in R statistical software [43] by using the raster, rgeos, and rgdal packages [44,45,46].

2.5. Burn Severity Classification Approach

Burn severity classes were built starting from 90 field plots, split into randomly selected training (n = 54, 60%) and validation datasets (n = 36, 40%). The training dataset was used for the non-linear regression model analysis, and to evaluate the statistical accuracy of the two satellite-derived dNBR reclassified raster images.

Owing to the non-linearity and the saturation of the relationship between field and satellite burn severity indices [42], the performance of the two dNBR as satellite-metrics of burn severity were tested by applying the non-linear regression model developed by Miller and Thode [47] for Mediterranean wildfires on both field-based indices (FBIs):

$$\mathrm{dNBR}=\mathrm{a}+\mathrm{b}·\exp (\mathrm{FBI}·\mathrm{c}),$$

(3)

where dNBR corresponds to the spectral index computed alternatively from Landsat-8 and Sentinel-2A, FBI is the field-based burn severity index computed using alternatively the CBI and GeoCBI protocols, and a, b, and c represent, respectively, the intercept, slope, and shape parameters of the exponential model. The non-linear regression models were compared using an array of goodness of fit indices: the root mean square error (RMSE), R squared, and Akaike Information Criteria (AIC), choosing the best non-linear regression model with the lower RMSE and AIC.

The four non-linear models of FBIs and dNBRs (built on the training dataset) were used to calculate the ranges of dNBR, corresponding to each CBI burn severity class. For each FBI threshold, the corresponding dNBR thresholds were computed using Landsat 8-OLI and Sentinel 2A MSI raster images. The raster images of both remote sources were then reclassified to assign to each pixel the corresponding specific burn severity class. Analysis was performed in R statistical software [43] by using the nlstool package [48].

2.6. Map Accuracy Assessment

The accuracies of the burn severity map produced by applying USDA-dNBR ranges and the regressed dNBR thresholds produced in the present study were evaluated by a discrete multivariate technique [36,49]. Specifically, in evaluating the accuracy of dNBR-based burn severity classification, a site-specific assessment was proposed using the error matrix approach [36]. The error matrix columns were set to be the reference burn severity classes derived from the field-based composite index, while the rows were set to be the dNBR-classified data derived from the remotely-sensed imageries. Classes were ordered by similarity, along columns and rows of each error matrix, according to the increasing degree of severity. To perform suitable statistical analysis, the error matrices were standardized [36] in order to compare accuracies in classification by the four alternative combinations of the FBIs (CBI and GeoCBI) and dNBRs (Landsat OLI-L8 and Sentinel 2A). For the purpose of this study, the process of error matrix standardization allows the accuracy among burn severity classes to be correctly compared, classified by all combinations of FBIs and dNBRs. In detail, the error matrices were standardized, dividing the diagonal (correct classification) and the marginal (row and column) by the total of each entire matrix.

Accuracy statistics were computed using the complete dataset (90 plots) for USDA-dNBR and the random test dataset consisting of 36 (40% of complete dataset) ground-based plots of burn severity for regression-estimated dNBR threshold ranges. The degree of accuracy of the burn severity raster map generated from the application of the dNBR thresholds was computed by the unweighted Kappa statistic [50], which represents a measure of agreement between the remotely sensed data classification and the ground-based reference data. Therefore, the unweighted K considers the matches between reference and classified data on the main diagonal of the error matrix. Additionally, the K values range from +1 to −1 [49].

Standardized deterministic overall producer’s and user’s accuracy statistics were computed in order to clearly evaluate the individual accuracy of each severity class in terms of both exclusion and inclusion errors. From a map producer’s viewpoint, accuracy determines the proportion of sites correctly mapped, while from the user’s point of view accuracy defines the percentage of classes correctly mapped. The statistical significance of each error matrix was tested by Z-test [36] at significance level α = 0.05. All analyses were performed in R statistical software [43] by using the psych and caret packages [51,52]

3. Results

As a result of the Z test, all analyzed error matrices are statistically significant at α = 0.05 (

Table 1. Error matrix statistics for each combination of ground-based burn index (composite Burn Index (CBI) and geometrically modified CBI (GeoCBI)) and reflectance-based delta Normalized Burn Ratio (dNBR) index derived from Landsat-8 OLI (dNBR_L8) and Sentinel-2A MSI (dNBR_S2). K represents the unweighted kappa value and vâr(K) the corresponding variance, Z is the score of Z test statistics and OA indicates the overall accuracy. Each statistic was computed on the complete field dataset (n = 90 plots). See Table S1 for category and ranges of the ground-based burn severity index, and Key and Benson [10] for categories and ranges of the dNBR index. Significance p-values of Z scores are represented as *, **, and ***, which correspond respectively to p < 0.05, p < 0.01 and p < 0.001.

Error Matrix	K	vâr(K)	Z	OA (%)
CBI - dNBR_L8	0.21	2.90·10−3	3.98 ***	33
GeoCBI - dNBR_L8	0.16	3.24·10−3	2.76 *	30
CBI - dNBR_S2	0.18	2.77·10−3	3.41 ***	30
GeoCBI - dNBR_S2	0.15	3.11·10−3	2.71 *	29

). Meanwhile, both matrices, which have the CBI as reference data, are highly significant, showing a p-value < 0.001 (Table 1). Furthermore, all the kappa values (K) are low, ranging from 0.15 to 0.21, showing that there is very poor agreement between reference data (CBI or GeoCBI) and classified data (dNBR by L8 and S2) (Table 1, Figure 2). Yet the overall accuracies (OA) were low and fell within the narrow range 29–33% (Table 1).

Based on the ranges of burn severity developed by Key and Benson [10] both the dNBR maps were classified. At class level, the producer’s and user’s accuracy steeply decrease from high to low severity classes in both classified dNBR maps (

Table 2. Producer’s and User’s accuracy percentage values (PA and UA respectively) by burn severity categories and for each combination of the ground-based burn index (CBI and GeoCBI) and the reflectance-based dNBR index derived from Landsat-8 OLI (dNBR_L8) and Sentinel-2A MSI (dNBR_S2). See Table S1 for categories and ranges of the ground-based burn severity index, and Key and Benson (2006) for categories and ranges of dNBR index.

Error Matrix	Accuracy	Burn Severity
		High	Moderate High	Moderate Low	Low	Unburnt
CBI - dNBR_L8	Producer’s	75	9	30	0	70
	User’s	92	60	42	0	33
GeoCBI - dNBR_L8	Producer’s	73	8	17	0	75
	User’s	85	40	26	0	43
CBI - dNBR_S2	Producer’s	75	15	15	0	60
	User’s	92	71	27	0	29
GeoCBI - dNBR_S2	Producer’s	73	12	10	14	67
	User’s	85	43	20	3	38

). In the high severity level, the producer’s accuracy ranged between 73–75%, irrespective of the error matrix analyzed. Yet user’s accuracies were higher (92%) in both error matrices, which include CBI as reference data, and the two error matrices, which include GeoCBI as reference data (85%).

For these two error matrices, the producer’s accuracy of the moderate-low burn severity level was in a very narrow range from 17% to 30%. Unlike the producer’s accuracy, user’s accuracies are slightly higher and lie in a wide range of values from 20% to 42%. Yet at moderate-high severity levels, the user’s accuracies were higher for both error matrices that have CBI as reference data than the error matrices that have GeoCBI as the reference data. Low severity levels were detected exclusively by the GeoCBI-dNBR_S2 error matrix with poor statistical accuracy of 14% and 3%, respectively, for the producer’s and user’s accuracy. The unburnt forest surfaces were more consistently retrieved by the Landsat-8 OLI. In both these error matrices, the producer’s accuracies were moderately high, with percentage values respectively of 70% and 75%. These higher producer’s accuracies were also observed for Sentinel-2A MSI, but with lower percentage values than Landsat-8 OLI. Unlike producer’s accuracies, user’s accuracies were lower, lying in the range 33–43%.

Comparison of the goodness of fit statistics shows that the relationship between the two ground-GeoCBI and dNBR metrics derived from Sentinel-2A proved to be the highest when modeled by the Miller and Thode [47] equation model (

Table 3. Goodness of fit statistics and parameters of regression models estimated by combining field-based (CBI and GeoCBI) and satellite-based (dNBR L8 and dNBR S2) burn severity metrics. Each regression model was estimated from the training dataset (n = 54 plots). R², RMSE and AIC correspond respectively to the coefficient of determination, root mean square error, and Akaike’s information criterion. For each regression parameter a, b, and c estimated values, standard error (±SE), student t-value and p-value are reported.

Regression Variables	Statistic	Estimate	Parameter	Value	±SE	t-Value	p-Value
dNBR L8 vs. CBI	R2	0.62	a	60.96	53.87	1.13	0.2631
	RMSE	173.6	b	10.13	11.07	0.92	0.3644
	AIC	718.6	c	1.49	0.37	4.02	0.0002
dNBR L8 vs. GeoCBI	R2	0.57	a	6.08	95.11	0.06	0.9493
	RMSE	185.5	b	44.54	48.30	0.92	0.3608
	AIC	725.8	c	1.01	0.36	2.78	0.0076
dNBR S2 vs. CBI	R2	0.63	a	30.57	58.71	0.52	0.6048
	RMSE	171.8	b	19.33	17.02	1.14	0.2614
	AIC	717.6	c	1.30	0.30	4.40	0.0001
dNBR S2 vs. GeoCBI	R2	0.69	a	−40.23	104.40	−0.39	0.7016
	RMSE	163.6	b	68.76	60.60	1.13	0.2619
	AIC	712.2	c	0.90	0.29	3.12	0.0029

, Figure 3a). However, all the models yielded a moderately high R² ranging between 0.57 and 0.69. The R² decreasing statistic is simultaneously associated with an increase in the RMSE and AIC estimates.

Their smaller values indicate greater parsimony of the non-linear regression model between the GeoCBI and Sentinel-2A based dNBR. Moreover, the two non-linear models in which the regressed variable was the dNBR obtained from Landsat 8 OLI show a slightly worse performance than the models in which the dNBR was obtained from Sentinel-2A MSI. Among the estimated model parameters, the shape/curvature parameter was the only one that was statistically significant for all the combination of ground-based (CBI and GeoCBI) and reflectance-based (dNBRs) burn index metrics (Table 3). Both the parameters intercept and slope are statistically non-significant, following the Student t-test.

The error matrix built from the discretized metrics GeoCBI and Sentinel-2A MSI-based dNBR and based on the Z-test was statistically significant (

Table 4. Error matrix statistics for the single combination of the ground-based burn index (GeoCBI) and reflectance-based dNBR index derived from Sentinel-2A (dNBR_S2). K represents the unweighted kappa value and vâr(K) the corresponding variance, Z is the score of the Z test statistic, and OA indicates the overall accuracy. Each statistic was computed on the validation dataset (n = 36 plots). See Table S1 for category and ranges of the ground-based burn severity index and Figure 3 for categories and regressed thresholds of the dNBR index. Significance p-values of Z scores are represented as *, **, and *** which correspond respectively to p < 0.05, p < 0.01 and p < 0.001.

Error Matrix	K	vâr(K)	Z	OA (%)
GeoCBI - dNBR_S2	0.77	6.90·10−2	9.32 ***	83

), indicating that the produced classification is meaningfully different from a random result.

Therefore, the resulting K value of 0.77, significantly greater than zero, indicates agreement between the remotely sensed classification based on dNBR derived from Sentinel-2A satellite and ground-based GeoCBI reference data. Moreover, the overall percentage accuracy between dNBR remote estimates and GeoCBI field-measured levels of burn severity was moderately high at 83% (Table 4, Figure 3b). Therefore, the resulting burn severity map produced showed a moderate-high overall accuracy (OA) and agreement (K) with reference GeoCBI burn severity levels (Figure 4).

Based on the validation dataset, all the classes ranging from the high to low level of severity showed moderate to high values of producer’s and user’s accuracy statistics (

Table 5. Producer’s and User’s accuracy percentage values by burn severity categories and for the single combination of ground-based burn index (GeoCBI) and reflectance-based dNBR index derived from Sentinel-2A (dNBR_S2). See Table S1 for categories and ranges of the ground-based burn severity index, File S1 for details on error matrix, and Figure 3 for categories and thresholds of the dNBR index.

Error Matrix	Accuracy (%)	Burn Severity
		High	Moderate High	Moderate Low	Low	Unburnt
GeoCBI - dNBR_S2	Producer’s	100	92	88	100	33
	User’s	100	80	70	100	100

). The producer’s accuracy was high in all classes except for the unburnt level, in which it decreased considerably to 33%, while it was consistently high from the user’s point of view, with values ranging from 70% to 100%.

By overlapping map of forest types to the maps of burn severity obtained from Sentinel-2A MSI, the area burnt by summer wildfires in 2017 amounted to 3350.23 ha, corresponding to ~88% of the National Park forest surface area (Figure 5). The Mediterranean conifer forest category was the highest burnt area (1426.58 ha), followed by the broadleaved and shrubland categories (1250.92 ha and 672.58 ha, respectively).

The high fire severity level accounts for 429.14 ha, 11% of the total forested area. About 86% were conifer stands and 14% shrubland and broadleaved categories. The moderate-high severity level accounts for 1170.77 ha, 31% of the total forested area, shared almost equally between conifer stands and broadleaved and shrubland categories.

The moderate-low severity level was 508.69 ha (13.4%) of the total forested area, affecting mainly the broadleaved stands (42%), followed by Mediterranean conifers (35%), and shrubland (23%). The low burn severity level represents the largest burnt surface area with 1241.63 ha, which accounts for 33% of the total forested park area. This severity level was experienced by the broadleaved stand category (60%), followed by Mediterranean conifers (23%) and shrubland stands (16%).

The unburnt surface accounts for 447.71 ha, corresponding to 12% of remaining surfaces of the forest park, represented by broadleaved stands (67%), Mediterranean conifers (20%) and shrubland categories (13%).

4. Discussion

4.1. Burn Severity Map

The dNBR map obtained from Sentinel 2A reflectance data provided moderate-high correlation with GeoCBI field-based estimates of burn severity and likewise provided moderate-high statistical accuracies in classifying burn severity levels, in the heterogeneous Mediterranean forest landscape in question. Combining GeoCBI and Sentinel-2A derived dNBR, the map of burn severity produced is characterized by a moderate agreement [53], as suggested by the K-statistic value (K = 0.77 of Table 4). This accuracy implies that the combination of GeoCBI and Sentinel-2A dNBR may be suitable for the burnt Mediterranean vegetation communities studied. Our finding was in agreement with two recent comparative studies on satellite ability to discriminate levels of burn severity in Mediterranean heterogeneous physiognomy [28,29].

The regressed dNBR-S2A threshold values estimated for discrimination of the high to moderate-high burn severity classes are more similar to the equivalent dNBR threshold values proposed by Key and Benson [10]. This similarity in the dNBR thresholds produces a high producer’s and user’s accuracy in both cases (Table 2; Table 5), suggesting that dNBR values proposed for North American temperate forests [10] could also be used to discriminate the high burn severity classes in Mediterranean forest ecosystems.

Although thresholds and maps of burn severity have been produced in different biomes, from tropical to boreal, most research has been carried out in temperate and Mediterranean ecosystems [37]. Consequently, different dNBR thresholds values for these biomes are reported [54,55,56,57]. Since dNBR thresholds are indeed influenced by an array of factors related to both the satellite used and the ground-based index, we consider our result potentially not generalizable. In contrast, when the degree of burn severity decreases, from moderate-high to low burn severity classes, the ability of the Key and Benson [10] thresholds to retrieve levels of severity becomes inconsistent, as documented by very low values of all accuracy statistics (Table 1; Table 2), irrespective of the satellite-based dNBRs. Indeed, the thresholds estimated by the non-linear regression approach are markedly different from the Key and Benson [10] dNBR thresholds. In detail, we estimated lower dNBR thresholds than those proposed by the authors, showing a decrease of 48.4%, 74.4%, and 25.9% from the moderate-high to lower levels of fire severity. Further, our results suggest that there is no set of dNBR thresholds that is appropriate for all burnt community types, especially for the decreasing level of fire severity.

In some map applications, such as those of land use classification [58], spatial units were aggregated to improve accuracy at the expense of thematic details. In burn severity assessment, it is well recognized that putting together a wide range of burn severity level reduce class numbers, increasing consistency at the expense of map utility [56]. In our study, we define the burn severity levels in order to avoid any misinterpretations of the burn severity map produced and increase their suitability for planning the most appropriate post-fire restoration activities. In the remote sensing literature, different numbers of classes and different ground CBI scores are used to characterize burn severity spatial variability [54,57]. In this context, although some difficulties in defining threshold values to classify burn severity into discrete classes have been documented [59], the number of ground-based categories and the breakpoint values should be made uniform at least for homogeneous physiognomy and biomes to ensure suitable comparison between the burn severity maps produced.

4.2. dNBR Performance

Although the dNBR approach has been questioned because it was initially developed for detecting the burnt vegetation community rather than estimating spatial variability of fire severity levels [14], we found it to be a suitable tool to discriminate the burn severity degree in a heterogeneous Mediterranean forest landscape where mixed-severity effects occur. The sensitivity of the dNBR index to discriminate wildfire effects on vegetation is strongly related to the NIR and SWIR reflectance bands [37,60,61], as they are sensitive to FCOV and LAI [62]. The better performance achieved in retrieving burn severity by Sentinel-2A is consistent with its higher spatial resolution as compared to Landsat-8. Moreover, high efficiency in discriminating fire effects by both the Sentinel-2A MSI narrow NIR band 8, and longer SWIR band 12, was recently reported [28,63]. After a wildfire, NIR reflectance decreases, and SWIR reflectance increases, due to alteration in anatomy and internal structure of fire-stressed leaves, which modify optical properties [64]. Moreover, the reduction in NIR reflectance can be related to a decrease in green biomass [65,66]. In contrast, SWIR band wavelengths are preferentially absorbed by water [21] and thus sensitive to canopy leaf water content when land surfaces are homogeneously covered by vegetation, and to soil moisture when land surfaces are bare or partially covered by vegetation. The contrasting signal in NIR and SWIR reflectance provided by healthy vegetation versus fire-disturbed areas and the different degrees of the resulting signal justify the use of the NBR, and more precisely of the dNBR, in assessing and classifying burn severity.

4.3. Field-Based Burn Index Performance

The more consistent relationship found between satellite-derived dNBR and GeoCBI than between satellite-derived dNBR and CBI, was due to the two spatial extent fractional cover (FCOV) and leaf area index (LAI) attributes employed in the GeoCBI protocol [15]. Indeed, these two visually assessed variables conferred on the GeoCBI a high sensitivity to the reflectance signal response at forests with different degrees of cover in each layer. In contrast, the CBI, originally designed to be operationally retrieved medium-resolution remotely sensed data [10], is most affected by the presence of the overstory canopy and low fractional cover [15]. This was supported by our scores of GeoCBI for moderate-low and low-severity fires, which were lower than the corresponding CBI scores (Figure 2).

Another possible positive attribute of GeoCBI that increases its consistency may be related to linkage with heat pulse of wildfires. It is widely reported that the spatial pattern of burn severity is strongly related to the stand structure and their specific composition [31,61]. In addition, the close spatial proximity of high-severity and low or moderate burnt areas was due to the local topography and weather conditions during the fire. Topographic traits markedly interact with local climate conditions, directly affecting fire propagation [61,67]. These traits, together with the seasonal variation in the moisture content of live [68] and dead fuels [69], and thus their flammability [70,71], had generated a patchwork of mixed-severity wildfires [72,73]. In such situations, the GeoCBI seems able to capture the horizontal and vertical variability of the mixed-severity burnt stands.

4.4. Model Regression and Thresholds

For the wildfires analyzed in this study, the non-linear regression model proposed by Miller et al. [47] provided a realistic and appropriate characterization of burn severity. Overall, non-linear regressions showed moderately high coefficients of determination (Table 3), although the model between GeoCBI ground–based and dNBR Sentinel-2A reflectance-based data exhibited the best performance. These results agree with previous studies, which clearly showed the non-linear relationship between reflectance-based dNBR and ground-based burn severity data [25,74]. On examining the non-linear relationship between GeoCBI and dNBR-S2A (Figure 3a) some defectiveness was found. One of these is related to the insensitivity of the dNBR to the unburnt forest. Indeed, the point cloud around the fitted non-linear regression line in correspondence with GeoCBI scores of zero showed high dispersion, as reported by others in burnt Mediterranean forest ecosystems [28]. This is probably the reason behind the circumstance to apply dNBR only within areas already mapped as burnt [47,55]. This relatively high dispersion of the point cloud strongly affected the intercept parameter of the non-linear regression model, which in turn is statistically non-significant (Table 3). Another defectiveness that attracted our attention concerned the asymptotic behavior of the model from GeoCBI scores lower than ~ 0.5. As a consequence of this behavior, the model potentially underestimated the low and unburnt severity classes, producing in turn a narrow range of dNBR values, in which pixels with low burn severity may be classified. However, the high accuracy in low severity level (Table 5) could be affected by the limited number of field-based samples available for the validation. When the score of ground GeoCBI was higher than ~ 0.5, the ability of the model to discriminate the level of burn severity markedly increased. As a consequence, the empirical model was highly capable of recognizing moderate-low and moderate-high burn severity levels. This ability was related to the curvature/shape parameter of the non-linear regression model that represents the only one that is statistically significant (Table 3). In the literature, a solution for asymptotic behavior was extensively investigated by applying a set of non-linear models based on the saturated growth model [25,47]. The model of Miller et al. [47] is effectively able to solve saturation problems when the value of the field-based burn index exceeds ~2.3 to 2.5 [25] rather than when the index decreases from 0.5 to approximately 0. The increase in SWIR reflectance after fires ceases with field-based index values above ~2 to 2.5, whereas the decrease in NIR reflectance is consistent with increasing burn severity [75], epitomizing the non-linear relationship between field-based burn and reflectance-based dNBR metrics.

An additional consequence of the asymptotic behavior exhibited by the non-linear relationship between GeoCBI and Sentinel 2A-dNBR for decreasing values of GeoCBI (<0.5) is low capability in discriminating and classifying unburnt-to-lightly-burnt patch areas (Figure 3, Table 5). In a mixed-severity fire experienced by Mediterranean vegetation, identifying the size of such patch areas and their patterns provides crucial information about the spatial extent of ecological refugia for fauna and sources of propagules (seeds or buds) for post-fire regeneration [76]. Quantitative analysis of the size, pattern, and proportion of the unburnt patch areas shows a close proximity of the severely burnt area to unburnt or lightly burnt area mainly in a mixed-severity fire [77]. In a protected area, such as the Vesuvius National Park, discriminating the spatial pattern of unburnt to lightly burnt patch areas within fire perimeters is necessary to understand successional ecological processes and their reestablishment following wildfires. Indeed, from a restoration/conservation perspective, these unburnt to lightly burnt areas act as biological legacies [11], providing habitat from which residual species can repopulate the neighboring burned area.

4.5. Factors Influencing the dNBR

One of the main disadvantages of using the bi-temporal spectral index approach, such as the dNBR, is the error increase associated with differences in vegetation phenology or weather conditions due to image data acquisition [55]. In Mediterranean evergreen tree species, phenology is related to the amount of foliage on crowns, which in turn depends on the annual production of new leaves/needles and on the longevity of old leaves. Leaf longevity is 2–3 years for Mediterranean conifers and between 1 and 2 years in evergreen oaks [78,79], which leads to a transient coexistence of old and new leaves cohorts in the crown. In the Mediterranean climate, old leaves become dry and senescent in the summer (leaf turnover approximately occurs from the beginning of July to the end of August), making up about one-third of the canopy. Since this dry matter accumulation in the crown coupled with decreasing live leaf water content [68] can affect optical properties of the tree crown and consequently pre- and post-fire pixel-reflectance response, we advise caution on the timing of pre- and post-fire data acquisition. In our study, we minimized the effect of these two factors by acquiring satellite imageries in two successive years (before and after the wildfire event) with comparable phenological phases and weather conditions (mid to late October). Hence, we believe that the high variability observed in the dNBR-GeoCBI relationships could be attributable to the alteration of reflectance signals caused by the impact on both crown trees and understory vegetation of a massive attack by an alien scale on P. pinea trees.

The pine tortoise scale Toumeyella parvicornis, a Neartic species, first became established in 2014 in southern Italy [80]. Large populations of this harmful pest infest the shoots (current and previous year) of Pinus pinea and P. pinaster, and produce profusely sticky honeydew. This sugar-rich sticky liquid, secreted by the insect, drips and covers all plant surfaces, such as needles, twigs, branches, and stem bark, as well as all the understory strata under pine crowns and the forest floor [81]. Honeydew is a host medium allowing the growth of a sooty mold fungus that blackens pine needles and all surfaces which it covers. Sooty mold is a collective name for several saprophytic fungal species [82]. The black encrustation on the surface of needles, tree crown, and of the understory strata is responsible for changes in the canopy reflectance and causes an appreciable decrease of the wavelength in the NIR range [83,84,85]. Simultaneously the greenness of canopy cover decreases, mimicking, in turn, a spectral signal alteration similar to that of light wildfires on vegetation. The combined effect of black encrustation of pine needles and sooty mold colonization of honeydew is responsible for the high deviation of GeoCBI from the regression model prediction (see RMSE and AIC values in Table 3). Indeed, two types of defectiveness may be observed in dNBR values. The first is a low change in delta between NBRpre and NBRpost wildfire occurrence that leads a low-severity fire to be wrongly classified as a moderate-severity fire. In contrast, the second type of distortion is the high value of dNBR, which leads a low-severity to be ascribed to a moderate-severity fire. As a consequence, in our case these reflectance alterations produce a reduced overall accuracy of the map, mainly in the range of fire severity classes from unburnt up to moderate-low.

The massive attack of the alien pine tortoise scale species can also explain the higher degree of burn severity found in conifers than in broadleaved stands and shrublands (Figure 5). Even though the fuel in the conifer stands was vertically discontinuous, the combination of litter- and retention of three-year-old dead needles in the crown, copiously coated with honeydew, modified fire behavior, favoring the spread of the surface fire to the crown. It cannot be excluded that the flames of the surface fire were borne towards the crown along the honeydew-coated bark-stem (Video S1). Indeed, pine tortoise scale secretions (dropped honeydew sensu Tamaki [86]) chiefly consist of sugars. Usually fructose, glucose, and oligosaccharides are the most abundant chemical component in Mediterranean conifer forest honeydew [87] and, to a lesser extent, terpenes [88], both characterized by high flammability [89]. In addition, the closed canopy cover of the Mediterranean pine stands under study amplified the torching and propagation to crown fire. In recent years, some studies have found that a forest invaded by non-native and invasive insects, or pathogens, is more frequently affected by severe wildfires [90]. From an ecological perspective, this phenomenon can be defined as a linked disturbance [91], which is the interaction between two primary disturbances, in our case the attack of pine tortoise scale and wildfire, which modifies fire behavior, amplifying the high-severity burnt areas.

5. Conclusions

In this study, a burn severity map was produced to spatially identify the forest areas severely damaged by multiple wildfires in 2017, in order to plan the most appropriate post-fire restoration activities efficiently. The dNBR obtained from the Sentinel-2A MSI satellite was found to be closely related to field GeoCBI estimates of burn severity. It outperformed Landsat 8 OLI, demonstrating its potential use in retrieving mixed-severity fire effects experienced by spatially heterogeneous Mediterranean forest ecosystems.

The accuracy of our results is related to both field- and satellite-based data, which, in turn, are both imperfect proxies of burn severity. Indeed, the field indices of burn severity are based on the quali-quantitative judgement procedure that lacks absoluteness [92] while several noise factors hamper satellite image analysis. Additionally, the dNBR thresholds were developed according to the degree of alteration experienced only by the vegetation component of the ecosystem, whereas it is well recognized that the degree of severity can be affected by several other altered ecological components.

Over large burnt forest areas, remotely-sensed severity of fires represents a valid alternative to field samplings and aerial photointerpretation. Knowing the spatial heterogeneity of burn severity could help land-managing institutions to identify action priorities for burnt areas. Indeed, after fire, forest managers need to map the different levels of burn severity in order to schedule their restoration efforts, mostly in Mediterranean areas highly prone to soil erosion affected by autumn rains. From an action priority perspective, after wildfires, the first information as a basis for operations concerns the spatial extent of high-severity burnt vegetation communities. In this context, remotely-sensed burn severity information is capable of detecting rapidly, and accurately, the high-severity burnt forest stands. Therefore, burn severity maps based on satellite-derived spectral indices could represent a useful tool in post-fire management practices, especially in spatial heterogeneous Mediterranean forest communities prone to fire.