Can Landsat-Derived Variables Related to Energy Balance Improve Understanding of Burn Severity From Current Operational Techniques?

Abstract

Forest managers rely on accurate burn severity estimates to evaluate post-fire damage and to establish revegetation policies. Burn severity estimates based on reflective data acquired from sensors onboard satellites are increasingly complementing field-based ones. However, fire not only induces changes in reflected and emitted radiation measured by the sensor, but also on energy balance. Evapotranspiration (ET), land surface temperature (LST) and land surface albedo (LSA) are greatly affected by wildfires. In this study, we examine the usefulness of these elements of energy balance as indicators of burn severity and compare the accuracy of burn severity estimates based on them to the accuracy of widely used approaches based on spectral indexes. We studied a mega-fire (more than 450 km² burned) in Central Portugal, which occurred from 17 to 24 June 2017. The official burn severity map acted as a ground reference. Variations induced by fire during the first year following the fire event were evaluated through changes in ET, LST and LSA derived from Landsat data and related to burn severity. Fisher’s least significant difference test (ANOVA) revealed that ET and LST images could discriminate three burn severity levels with statistical significance (uni-temporal and multi-temporal approaches). Burn severity was estimated from ET, LST and LSA using thresholding. Accuracy of ET and LST based on burn severity estimates was adequate (κ = 0.63 and 0.57, respectively), similar to the accuracy of the estimate based on dNBR (κ = 0.66). We conclude that Landsat-derived surface energy balance variables, in particular ET and LST, in addition to acting as useful indicators of burn severity for mega-fires in Mediterranean ecosystems, may provide critical information about how energy balance changes due to fire.

1. Introduction

Fire is one of the major disturbance processes in many ecosystems and is particularly prevalent in Mediterranean ecosystems, where it impacts biodiversity [1,2] and other important ecosystem properties. The frequency, severity and size of fires have been increasing in recent decades in Mediterranean ecosystems, due mainly to climate change and modifications in land use [3]. To evaluate fire damage, one of the most commonly used metrics is fire/burn severity [4]. Burn severity groups short- and long-term fire effects on vegetation and soil [5]. Key and Benson [6] distinguished between initial assessment of burn severity, which aims to map burn severity immediately after fire, from an extended assessment, whose goal is to map burn severity when survivorship and mortality are detectable. Increases in burn severity due to changes in climate and land cover may modify vegetation characteristics (composition, resilience, structure) and soil attributes, even in Mediterranean ecosystems that are well adapted to fire [7,8].

Fire reduces evapotranspiration (ET) and increases surface runoff following rainfall [9]. Usually, the increase in surface runoff is lower than the decrease in ET, causing an increase in soil moisture documented by different studies (e.g., [10,11]). However, forested areas affected by high burn severity show a decrease in soil moisture. In this case, the increase in surface runoff exceeds the reduction in ET and leads to drier post-fire soils [9,12]. Post-fire changes in vegetation, soil and water balance corroborate wildfires impact on energy balance [13,14,15]. Fire-induced modifications in vegetation structure and species composition alter latent heat flux and other variables of the energy balance equation [16,17]. Changes in Land Surface Temperature (LST), an important parameter of ET, have been recently identified as a potential indicator of burn severity [18,19,20,21]. Similarly, previous research studies showed the relationship between fire damage and land surface albedo (LSA) [20,22,23,24]. LSA determines the total amount of solar radiation that an ecosystem absorbs [25] and it is subject to seasonality [20]. To our knowledge, there is a lack of studies relating to all three elements of energy balance (ET, LST, and LSA) to burn severity, in particular in Mediterranean forest ecosystems. However, understanding their interactions is critical to evaluating the role of fire on Earth. The loss of vegetation modifies rainfall and water availability (more than 40% of rainfall over Earth’s land is due to ET on average) and has an impact on warming at local scales as well [26]. Forest fires may have an influence on regional climate by modifying energy balance [27,28], and some studies suggest that this may affect climate at global scales [29,30].

There is a variety of approaches to compute LST from thermal satellite data with a satisfactory accuracy (1K or better) and spatial resolution (30 × 30m) [31,32]. Similarly, LSA may be computed from reflective satellite data with a low residual error (less than 0.02) at the same fine spatial resolution [33,34]. However, many previous studies relating to fire effects and ET [13,15,35,36] were based on surface fluxes measured in the field, as it is difficult to estimate ET from satellite data with an adequate accuracy and spatial resolution. Thus, these results are only valid for vegetation similar to the vegetation analyzed by the instrument on a local scale [37]. ET requires a relatively high number of parameters to be computed with reasonable accuracy. Thus, a high temporal resolution (daily frequency or higher) is needed, whereas a fine spatial resolution is advisable. For this reason, different satellite sensors combined with meteorological ground stations are usually needed to estimate ET from satellite data [38]. Among the remote sensing-based ET models, surface energy balance techniques based on LST from thermal infrared data are widely used. Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC) [39] is a well-established surface energy balance-based model [40]. One of the main limitations of METRIC is, however, the need for a well-trained expert to calibrate and run the model [41]. The calibration of METRIC requires the identification of maximum and minimum ET within a satellite image. Some automatic calibration algorithms for the METRIC model have, however, been developed to avoid this limitation (e.g., [42,43]). Among them, the Earth Engine Evapotranspiration Flux (EEFlux) application was designed and implemented in the Google Earth Engine (GEE) platform [44]. The ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) mission has ET as one of the primary output variables [45,46], although at the present it only collects data over the conterminous United States (CONUS) as well as key biomes and agricultural zones around the world and selected FLUXNET [47] validation sites.

Similarly, remote sensing techniques provide a widely used alternative to field-measured burn damage [14,48,49,50,51,52,53]. Operational examples include the European Forest Fire Information System (EFFIS) and the Monitoring Trends in Burn Severity (MTBS) project in the USA. Among the most widely used techniques to estimate burn severity from satellite data, spectral indices stand out due to their simplicity [52,53,54,55,56,57,58]. In particular, using thresholds to classify differenced Normalized Burn Ratio (dNBR) [6] has become a standard to estimate burn severity from optical remotely sensed data [59], specifically from Landsat data [57,60,61]. Additionally, differenced Normalized Difference Vegetation Index (dNDVI) may be used to map burn severity when it is not possible to calculate dNBR due to sensor characteristics [62]. As most of the sensors onboard Unmanned Aerial Vehicles (UAVs) do not record information in shortwave infrared wavelengths, NDVI, and in particular dNDVI, is emerging as an accurate index to estimate burn severity [63]. However, some studies have pointed out some limitations of using spectral indices as a base to estimate burn severity: their interpretation in terms of fire damage is not easy due to the absence of standardized units [64,65] and due to a nonlinear relationship to burn severity [66].

In this context, this research work aims to analyze the relationship between burn severity and variations in ET, LST and LSA after large fires in Mediterranean forests at fine spatial resolution and to verify whether they may act as indicators of burn severity (initial and extended) in these ecosystems. Spectral indexes widely used to estimate burn severity (dNDVI, dNBR) were used as additional variables. In particular, our specific objectives are: (1) to evaluate the changes due to fire in the five variables (ET, LST, LSA, NDVI and NBR) over time: immediately after fire, a month after fire, two months after fire, and one year after fire; (2) to determine whether the five variables can be regarded as indicators of burn severity (initial and extended assessments) from both uni-temporal and multi-temporal perspectives and how many burn severity levels may be distinguished by each variable; (3) to identify pre-fire factors that may also cause changes in the variables: in particular, pre-fire species composition (highly related to fuel model), climate and topography; and 4) to determine the accuracy of the burn severity estimates from these variables. We structured Results and Discussion sections by three main questions that are closely related to these specific objectives: (1) How does fire modify ET, LST, LSA, NDVI and NBR?; (2) How are ET, LST, LSA, NDVI and NBR post-fire images influenced by burn severity and pre-fire factors?; and (3) Is it possible to accurately estimate burn severity from the ET, LST and LSA images?

2. Materials

2.1. Study Area

Our study area is the mega-fire of Pedrogão Grande and environs that burned in 2017 in central Portugal (Figure 1, left). The Joint Research Centre of the European Commission included this fire in its annual report as the fire that caused a total of 66 deaths and important losses in term of forest areas and buildings [67]. In particular, 263 structures were inventoried as being damaged by Ribeiro et al. [68]. The Pedrogão Grande fire burned 458.93 km² between June 17^th 2017 and June 24^th 2017 [69]. After months with low rainfall and high temperatures, 80% of the region was in extreme drought on June 17^th according to the Instituto Português do Mar e da Atmosfera (IPMA). With these conditions, a multi-cell thunderstorm emerged, and lightning and strong winds caused the fire to burn out of control [70,71].

Climate in the area can be described as a mix of Csa (typical Mediterranean) and Csb (humid Mediterranean) Köppen-Geiger classes [72]. Rainy winters and hot and dry summers characterize the typical Mediterranean climate (annual precipitation is between 400–750mm). The humid Mediterranean climate, however, shows a higher annual precipitation (750–2000mm), with summers that are mild and dry.

The highest elevations are to the north and northwest of the affected area. In progression to the south, lower elevations predominate (approximately 400m around Pedrógão Grande). Orography is very heterogeneous. Very steep slopes with values between 45% and 60% or higher predominate in the north and northeast, whereas the area is flatter, with gentle slopes between 0% and 15%, in the south and southwest (basin of Zêzere river). In the extreme Southeast, there are slopes between 30% and 40%.

The Pedrógão Grande area is integrated into the Interior North Pinhal Zone, the largest forest spot in Europe before the fire. It is contained between the basins of the Zêzere and Unhais rivers and the banks of the Pêra and Mega rivers. In this area, plantations of eucalyptus (Eucalyptus globulus Labill.), and maritime pine forest (Pinus pinaster Ait.) are dominant (Figure 1, upper right). However, it is still possible to find cork oak forests (Quercus suber L.), chestnut trees (Castanea sativa Mill.), strawberry trees (Arbutus unedo L.) and shrubs such as gorse (Ulex europaeus L.), red heather (Erica australis L.), white heather (Erica arborea L.), and broom (Cytisus scoparius L.). Acacias and olive trees have great importance in the local economy. Agricultural areas were located close to villages, although both agriculture and grazing activities were low (1–3 livestock units per ha) [71].

2.2. Materials

We downloaded USGS Landsat 8 Operational Land Imager (OLI)/ Thermal InfraRed Sensor (TIRS) surface reflectance images (path/row 204/32, higher-level data product on demand) acquired on 15 June 2017; 1 July 2017; 17 July 2017; 2 August 2017; 18 August 2017; 19 September 2017; 5 August 2018; and 8 October 2018. From them, we computed NDVI, NBR, LST and LSA. METRIC daily ET images were obtained by using the software developed by University of Nebraska-Lincoln [44].

The official fire perimeter and burn severity map (three severity levels) were used as ground reference (estimated geometric accuracy 5m CE90 or better). They were based on Copernicus—Emergency Management Service maps (built using SPOT 6 and 7 imagery) verified in the field and from photointerpretation of Sentinel 2 MSI data by the Portuguese Study Center of Forest Fires [69].

As potential pre-fire factors that may also influence the variables related to energy balance, we included in the study: pre-fire species composition, climate and topography. The official Portuguese Land Cover Map (Carta de Ocupação do Solo, COS15) enabled us to define pre-fire vegetation species composition. COS15 has 48 classes in vector format and has a minimum cartographic unit of one hectare and a minimum distance between lines of 20 meters [73]. Köppen-Geiger climate classes [72] of the study area were extracted from the updated world map of the Köppen-Geiger climate classification [74]. Finally, a 30m Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM V2) provided by USGS (https://search.earthdata.nasa.gov) was employed to derive three topographic factors (elevation, aspect and slope).

3. Methods

3.1. LST Calculation

As in previous studies [75,76], LST was calculated based on the single-channel method proposed by Yu et al. [77]. This method is based on the radiative transfer equation (Equation (1)). See the review work of Li et al. [32] for detailed information on satellite-based methods to obtain LST. Calculations were implemented in the R statistical program (release 3.4.3) [78]

$$\mathrm{B}\left({\mathrm{T}}_{\mathrm{s}}\right)=\left[\frac{\left(\mathrm{B}\left(\mathrm{T}\right)-\mathrm{I}\uparrow \right)}{\mathsf{\tau }\left(\mathsf{\theta }\right)}-\mathrm{I}\downarrow \left(1-\mathsf{\epsilon }\right)\right]\frac{1}{\mathsf{\epsilon }}$$

(1)

where B(T_s) represents ground radiance; B(T), radiance of thermal band 10 with brightness temperature T; τ(θ), atmospheric transmittance for a view zenith angle θ; I↓, and I↑ respectively, denote downwelling and upwelling path radiance; and ε, land surface emissivity calculated from Sobrino et al. [79]. LST is computed from ground radiance according to Planck’s law using Equation (2)

$$\mathrm{LST}=\frac{{\mathrm{K}}_2}{\ln \left(1+{\mathrm{K}}_1/\mathrm{B}\left({\mathrm{T}}_{\mathrm{s}}\right)\right)}$$

(2)

where K1 and K2 are thermal constants defined in metadata of Landsat 8 OLI/TIRS images.

3.2. LSA Calculation

Similar to previous works [23,24], LSA was estimated from surface reflectance based on Liang [80] and Liang et al. [34]. We only considered the total shortwave broadband LSA (LSA) (Equation (3))

$$\mathrm{LSA}=0.356{\mathsf{\alpha }}_1+0.130{\mathsf{\alpha }}_3+0.373{\mathsf{\alpha }}_4+0.085{\mathsf{\alpha }}_5+0.072{\mathsf{\alpha }}_7-0.0018$$

(3)

where α_i indicates the albedo in the ‘i’ Landsat 8 OLI band that is approximated by reflectance.

3.3. ET (METRIC Model)

We include here a summary of the METRIC model. More detailed information about the METRIC-model can be found in Allen et al. [39] and Irmak et al. [81]. In the METRIC model, the energy consumed by the ET process is determined from satellite imagery using the surface energy balance equation (Equation (4)) on a pixel by pixel basis

$$\mathrm{LE}=R_n-G-H$$

(4)

where LE represents latent energy consumed by actual ET (i.e., heat energy used by water in its phase change from liquid to gas during the ET process); R_n, net radiation flux density (W/m²); G, ground heat flux density or sensible heat conducted into the ground) (W/m²); and H, sensible heat flux density convected into the air (W/m²). In METRIC, these three key variables (R_n, G and H) are estimated from satellite data. In particular: R_n from narrow-band reflectance and surface temperature; G, from R_n, surface temperature, and vegetation indices; and H, from surface temperature ranges, surface roughness, and wind speed (see [39]). Next, instantaneous ET (ET_i, expressed in mm/h) is calculated for each pixel using Equation (5)

$${\mathrm{ET}}_i=3600\frac{LE}{\lambda {\rho }_w}$$

(5)

where 3600 is a constant that converts seconds to hours; ρ_w is the density of water (Kg/m³); and λ, is the latent heat of vaporization, expressed in J/Kg, and estimated from surface temperature, expressed in K.

From ET_i, reference ET fraction (ET_rF) is calculated as the ratio of ET_i to the instantaneous tall crop reference ET (ET_r) computed from weather data (Equation (6)). Alfalfa is the tall crop used for METRIC as the tall alfalfa reference resembles maximum energy-limited ET from a well-watered, extensive surface of vegetation.

$$E{T}_{r}F=\frac{E{T}_i}{E{T}_r}$$

(6)

ET_rF is used for extrapolating instantaneous ET to daily ET (ET) expressed in mm/day (Equation (7)). In the METRIC process, it is assumed that the instantaneous ET_rF calculated at image time is constant for the 24-hour period. This assumption is generally valid for agricultural crops (developed to maximize photosynthesis and thus stomatal conductance). For some native vegetation under water stress, however, ET_rF may reduce its value during the afternoon. In this case, ET_rF is modeled from instantaneous ET_rF by developing the appropriate functions.

$$\mathrm{ET}={\mathrm{C}}_{rad}\times E{T}_{r}F\times E{T}_{r-24}$$

(7)

where C_rad represent the correction term used in sloping terrain, in horizontal areas C_rad equals 1.0; and ET_r−24, cumulative 24 h ET_r for the day of the image.

3.4. Database Construction

As a preliminary step to database construction, NBR [82] and NDVI [83] were computed from Landsat surface reflectance images. Additionally, we calculated the difference between pre- and post-fire situations to perform a multi-temporal analysis. In particular, we subtracted post-fire variables from pre-fire ones for ET, LSA, NBR and NDVI, and pre-fire variables from post-fire ones for LST. Variables ending with “_i” represent the initial assessment (difference between pre-fire -06/15/17- and immediately after fire images -07/01/17-), and “_e” an extended assessment (difference between pre-fire -06/15/17- and 1-year after fire images -08/05/18-).

We applied a mean 3 × 3 filter to satellite-derived images before extracting the digital values for the sampling points to minimize positional errors [6,19]. As recommended by Congalton and Green [84], a stratified random sampling was adopted. The number of training samples was related to the surface that each burn severity level occupies on the ground reference map: 2532 samples for high burn severity, 1638 for moderate burn severity and 721 for low burn severity. In addition, the zone outside the fire perimeter was sampled to define the unburned class; we defined 1213 sampling points for the unburned class. A kurtosis normality test enabled us to check its normality.

3.5. Statistical Analysis

After evaluating the statistical characteristics of the variables of interest (ET, LST, LSA NDVI and NBR) over time as a preliminary analysis, we aim to determine how spatial variation in the five variables is related to burn severity. One-way Analysis of Variance (ANOVA) was applied by grouping values into each burn severity level. A Fisher’s least significant difference (LSD) test determined which sample means were significantly different from the others. These analyses were performed at different moments to analyze the temporal evolution of variables: immediately after fire, a month after fire, two months after fire and 1-year after fire (summer and fall). ANOVA analysis enabled us to identify the number of burn severity levels that may be distinguished with statistical significance from the images and to determine whether those variables could be used as indicators of burn severity. In addition, to the uni-temporal post-fire analysis, a similar multi-temporal analysis was performed by subtracting the pre-fire image from the immediately after fire and 1-year after fire images (initial and extended assessment respectively).

Pre-fire conditions also have an influence on the variables of interest. To evaluate the contribution of those determining factors on the variables, a multifactor analysis of variance was performed. It searched for significant interactions amongst the factors. The p-value tests the statistical significance of the contribution of each measured factor, having removed the effects of all other factors. We included as pre-fire factors: dominant vegetation species; Köppen climate category and topographic variables (in particular, elevation, slope and aspect). Vegetation species is usually related to fire damage, although it is not a true measure of fuel state [85]. Specifically, three categories were taken into account in the vegetation species factor—eucalyptus, pine and shrubland—and two Köppen classes in the climate factor—Csa (Mediterranean climate) and Csb (humid Mediterranean climate). Regarding topographic factors, we defined five classes in elevation: 0- 300m, 301-600m, 601-900m and 901m-highest value; five classes in slope: 0°–5°, 5°–10°, 10°–20°, 20°–30°, 30°–highest value; and eight classes in aspect: north (0°–22.5° and 337.5°–360°), northeast (22.5°–67.5°), east (67.5°–112.5°), southeast (112.5°–157.5°), south (157.5°–202.5°), southwest (202.5°–247.5°), west (247.5°–292.5°), and northwest (292.5°–337.5°).

The relationship amongst the five variables was evaluated as well. As in previous studies (e.g., [19,86,87] Pearson’s correlation analysis was utilized to analyze the interactions among the five variables taken into account.

3.6. Burn Severity Mapping

Burn severity (initial assessment) was mapped from both uni-temporal (immediately after fire images) and multi-temporal (pre- and post-fire difference images) perspectives. Classification of the continuous images into four classes—high, moderate and low burn severity levels and unburned—was based on thresholds. The halfway values between the mean values of consecutive classes (computed from the training samples) acted as initial thresholds. These initial thresholds were refined manually to minimize the number of errors of classification [88].

Our study used error matrices and κ statistics [84] to measure the accuracy of each burn severity estimate. In addition, Overall Accuracy (OA), Producer’s Accuracy (PA) (omission error) and User’s Accuracy (UA) (commission error) were computed as well. Prior to error matrix computation, a new stratified random disproportional sampling was applied to the ground reference and to the classified images. In this case we defined 2357 sampling points for high burn severity, 1694 for moderate burn severity, 661 for low burn severity level, and 1140 for the unburned class.

4. Results

4.1. How Does Fire Modify ET, LST, LSA, NDVI and NBR?

Notable differences between burned and unburned areas in every post-fire variable are visually apparent (Figure 2). The images immediately after the fire (second column of Figure 2) showed such differences more clearly. The images one month and two months after the fire (third and fourth column) also enabled us to clearly distinguish burned from unburned areas. Moreover, the two months after the fire image (19 September 2017) enabled us to visually identify a new forest fire in the lower right corner. Even 1-year after fire (right column), burned areas were readily identified in every image except LSA. Moreover, a large forest fire (332 km²) which happened in Ponte das Portelinhas (Sertã) in 2 November 2017 is still evident in the 1-year after fire images. In addition, the capability of ET images to clearly delineate riparian vegetation along the Zêzere River should be noted. Spatial differences inside of the burn perimeter could be noticed as well in the five variables (more clearly in the images closer to the fire date), which suggests sensitivity to burn severity levels.

Table 1. Temporal evolution of statistical characteristics of ET. LST. LSA. NDVI and NBR for burned and unburned samples.

		06/15/17		07/01/17		08/02/17		09/19/17		08/05/18
		μ	σ	μ	σ	μ	σ	μ	σ	μ	σ
ET (mm/day)	Burned	-	-	1.62	1.62	1.62	0.95	1.52	0.98	1.76	1.01
	Unburned	5.58	1.34	6.65	1.16	3.47	0.94	3.26	0.79	3.13	1.21
LST (K)	Burned	-	-	315.79	4.90	319.51	4.72	308.96	3.99	315.25	2.47
	Unburned	305.13	3.80	299.36	3.24	304.10	3.91	298.38	4.05	309.57	4.19
LSA	Burned	-	-	0.08	0.02	0.09	0.03	0.09	0.02	0.16	0.01
	Unburned	0.12	0.03	0.12	0.03	0.11	0.03	0.10	0.03	0.15	0.01
NDVI	Burned	-	-	0.29	0.12	0.30	0.12	0.37	0.11	0.41	0.07
	Unburned	0.72	0.12	0.73	0.13	0.72	0.13	0.68	0.14	0.51	0.11
NBR	Burned	-	-	-0.12	0.24	-0.06	0.22	0.07	0.21	0.29	0.10
	Unburned	0.57	0.16	0.55	0.17	0.55	0.18	0.55	0.22	0.45	0.15

μ: mean value in sampling plots; σ: standard deviation in sampling plots; ET: evapotranspiration; LST: land surface temperature; LSA: land surface albedo; NDVI: normalized difference vegetation index; NBR: normalized burn ratio

displays the statistical characteristics (mean values and standard deviation) of ET, LST, LSA, NDVI and NBR for burned and unburned samples at different times: before the fire (only unburned samples), immediately after fire, one month, two months and 1-year after fire. Knowing the mean values of the burned and unburned samples for each image enabled us to minimize the influence of external and meteorological parameters in the fire-induced changes in the five variables. Fire caused an important initial ET decrease, which eventually reduced its value. ET declined seasonally in unburned samples as well. Patterns of LST variation after fire were similar to the ET trend, but showed a reverse pattern. LST values increased in burned samples. The increase was higher immediately after fire, although it could be observed even a year after the fire. Regarding LSA, burned areas decreased immediately after they burned relative to unburned areas, with the difference decreasing over time. A year after fire, LSA in burned samples increased relative to unburned samples, mainly due to vegetation regrowth. Finally, both NDVI and NBR showed the same tendency: a strong decrease in the burned samples immediately after the fire. These differences declined in the year after fire image. The standard deviation of burned samples immediately after fire was the highest of all images and implies a high spatial heterogeneity that could be related to burn severity.

4.2. How Are ET, LST, LSA, NDVI and NBR Post-Fire Images Influenced by Burn Severity and Pre-Fire Factors?

The results of one-way ANOVA applied to each post-fire image (

Table 2. Fisher’s least significant difference test for post-fire ET, LST, LSA, NDVI and NBR images and burn severity levels.

Evapotranspiration (ET, mm/day)
S	07/01/17		07/17/17		08/02/17		08/18/17		09/19/17		08/05/18		10/08/18
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	0.94	a	1.32	a	1.45	a	0.98	a	1.16	a	1.50	a	1.03	a
M	1.89	b	1.52	b	1.70	b	1.19	b	1.74	b	1.90	b	1.17	b
L	3.36	c	1.73	c	2.04	c	1.51	c	2.26	c	2.31	c	1.45	c
U	6.65	d	2.80	d	3.47	d	2.76	d	3.26	d	3.13	d	1.84	d
Land Surface Temperature (LST, K)
S	07/01/17		07/17/17		08/02/17		08/18/17		09/19/17		08/05/18		10/08/18
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	318.7	a	325.9	a	322.1	a	321.6	a	310.9	a	315.8	a	302.9	a
M	314.5	b	321.8	b	318.4	b	318.2	b	308.1	b	315.1	b	302.5	a
L	309.4	c	317.5	c	314.3	c	314.6	c	305.0	c	313.5	c	300.8	b
U	299.1	d	307.6	d	303.9	d	305.6	d	298.4	d	309.5	d	297.7	c
Land Surface Reflectance (LSA)
S	07/01/17		07/17/17		08/02/17		08/18/17		09/19/17		08/05/18		10/08/18
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	0.07	a	0.08	a	0.08	a	0.08	a	0.08	a	0.16	a	0.10	a
M	0.08	b	0.10	b	0.10	b	0.10	b	0.09	b	0.16	a	0.10	a
L	0.10	c	0.11	c	0.11	c	0.11	c	0.10	c	0.16	a	0.10	a
U	0.12	d	0.12	d	0.11	d	0.11	c	0.10	c	0.15	b	0.08	b
Normalized Difference Vegetation Index (NDVI)
S	07/01/17		07/17/17		08/02/17		08/18/17		09/19/17		08/05/18		10/08/18
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	0.21	a	0.21	a	0.25	a	0.28	a	0.33	a	0.40	a	0.54	a
M	0.30	b	0.29	b	0.32	b	0.34	b	0.38	b	0.41	a	0.55	a
L	0.49	c	0.44	c	0.48	c	0.49	c	0.51	c	0.45	b	0.62	b
U	0.74	d	0.71	d	0.73	d	0.72	d	0.69	d	0.51	c	0.71	c
Normalized Burn Ratio (NBR)
S	07/01/17		07/17/17		08/02/17		08/18/17		09/19/17		08/05/18		10/08/18
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	-0.28	a	-0.21	a	-0.21	a	-0.15	a	1.10	a	0.28	a	0.29	a
M	-0.05	b	0.00	b	0.01	b	0.06	b	1.69	b	0.29	a	0.30	a
L	0.23	c	0.22	c	0.26	c	0.29	c	2.26	c	0.35	b	0.41	b
U	0.57	d	0.55	d	0.56	d	0.56	d	3.25	d	0.46	c	0.53	c

S: burn severity level; H: high burn severity level; M: moderate burn severity level; L: low burn severity level; U: unburned; μ: mean value; HG: homogeneous groups, different letters in each HG column indicate a group of means within which there are no significant differences, showing significant differences among group means; white background: four levels were discriminated; light grey background: three levels were discriminated; dark grey background: two levels were discriminated

) indicated that the five variables tested showed significant differences (p < 0.01) among the mean values of each burn severity level when the 2017 images were taken into account, with the exception of LSA, whose low burn severity and unburned mean values could not be distinguished in the last two 2017 images. Therefore, three burn severity levels could be distinguished in all immediately post-fire and a month after fire images.

Regarding the 2018 images, only ET and LST could discriminate three burn severity levels with statistical significance. ET displayed significant differences among all burn severity levels in the two 2018 images (summer and fall). LST, however, was not distinct between high and moderate burn severity levels in the image acquired in fall 2018 (8^th October 2018). Similarly, NDVI and NBR only discriminated two burn severity levels (high-moderate and low) in the 2018 images. LSA showed the worst performance in the ANOVA analysis as it could only discriminate between the burned and unburned class in both 2018 images. Consequently, ET and LST were the only variables that enabled us to distinguish the three burn severity levels in the 1-year after fire summer image, and ET could also discriminate the three burn severity levels in the 1-year after fire fall image.

Table 3. Fisher’s least significant difference test for the difference between pre-fire and post-fire ET, LST, LSA, NBR and NDVI images and burn severity levels.

Initial Assessment (06/15/17–07/01/17)
S	dET_i		dLST_i		dLSA_i		dNBR_i		dNDVI_i
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	3.76	a	9.52	a	0.05	a	0.78	a	0.46	a
M	3.02	b	5.57	b	0.04	b	0.56	b	0.38	b
L	1.52	c	1.04	c	0.03	c	0.32	c	0.22	c
U	-1.07	d	-5.77	d	0.01	d	0.02	d	-0.01	d
Extended assessment (06/15/17 - 08/05/18)
S	dET_e		dLST_e		dLSA_e		dNBR_e		dNDVI_e
	μ	HG	μ	HG	μ	HG	μ	HG	μ	HG
H	3.20	a	6.84	a	-0.03	a	0.22	a	0.28	a
M	3.00	b	6.27	b	-0.04	a	0.23	a	0.28	a
L	2.57	c	4.98	c	-0.03	b	0.20	b	0.26	b
U	2.46	d	4.44	d	-0.03	b	0.12	c	0.21	c

S: burn severity level; H: high burn severity level; M: moderate burn severity level; L: low burn severity level; U: unburned; μ: mean value; HG: homogeneous groups, different letters in each HG column indicate a group of means within which there are no significant differences, showing significant differences among group means; white background: four levels were discriminated; light grey background: three levels were discriminated; dark grey background: two levels were discriminated

displays the results of the one-way ANOVA from a multi-temporal perspective. Regarding the initial assessment (difference between pre-fire and immediately post-fire image), the five variables (ET, LST, LSA, NBR and NDVI) performed similarly. All of them could discriminate the three burn severity levels with statistical significance. However, regarding the extended assessment (difference between pre-fire and 1-year after fire image), only dET_e and dLST_e could discriminate them. dLSA_e, dNBR_e, and dNDVI_e enabled only the discrimination of two burn severity levels (high-moderate and low). These results are in accordance with the results of the one-way ANOVA from uni-temporal perspective (Table 2).

Table 4. Multifactor analysis of variance for ET, LST, LSA, NDVI and NBR. Initial and extend assessment.

Initial Assessment (07/01/17)
	ET		LST		LSA		NDVI		NBR
Factors	p-value	%	p-value	%	p-value	%	p-value	%	p-value	%
Burn severity	0.000	83.90	0.000	83.87	0.000	40.53	0.000	85.60	0.000	85.66
Climate	0.107	0.00	0.000	0.43	0.197	0.00	0.415	0.00	0.170	0.00
Vegetation	0.000	1.49	0.000	0.77	0.000	17.47	0.000	1.60	0.000	1.22
Elevation	0.085	0.44	0.000	1.65	0.000	1.27	0.000	0.79	0.000	0.94
Slope	0.000	0.58	0.000	0.07	0.000	0.93	0.007	0.63	0.009	1.07
Aspect	0.000	6.20	0.000	4.44	0.000	8.75	0.000	2.04	0.000	0.24
Total explained		92.61		91.23		68.95		90.66		89.13
Total error		7.39		8.76		31.05		9.34		10.87
Extended assessment (08/05/17)
	ET		LST		LSA		NDVI		NBR
Factors	p-value	%	p-value	%	p-value	%	p-value	%	p-value	%
Burn severity	0.000	37.32	0.000	47.24	0.000	2.52	0.000	27.13	0.000	32.47
Climate	0.000	8.01	0.564	0.00	0.567	0.00	0.074	0.00	0.323	0.00
Vegetation	0.000	0.65	0.000	6.02	0.000	20.13	0.000	10.70	0.000	11.08
Elevation	0.000	10.49	0.000	4.38	0.000	6.17	0.002	5.78	0.006	4.17
Slope	0.000	4.23	0.000	1.75	0.000	5.32	0.000	0.94	0.001	1.06
Aspect	0.000	4.98	0.000	8.85	0.000	6.43	0.000	3.05	0.000	4.03
Total explained		65.68		68.24		40.57		47.60		52.81
Total error		34.32		31.75		59.43		52.40		47.19

p-value: p-values test the statistical significance of each of the factors, p-value < 0.01 indicates a significant contribution at the 99% confidence level; %: percentage of variance contributed by each of the factors

displays the results of the multifactor analysis of variance for ET, LST, LSA, NDVI and NBR that used post-fire burn severity, pre-fire vegetation species, climate, elevation, slope and aspect as factors. We observed that burn severity contributes to more than 80% of ET, LST, NDVI and NBR images (see Table 4, initial assessment), with pre-fire vegetation species and aspect being the next two factors with a high contribution. The total variance explained for the factors in those four variables was approximately 90%. With regards to the extended assessment (images acquired 1-year after fire), the total explained variance showed a decrease when compared to the initial assessment. Burn severity contribution also decreased. The topographic contribution increased for all of the variables. The vegetation species factor also increased its contribution to all variables except ET. Climate only had a relevant contribution to ET variance. LSA displayed the lowest values for the total explained variance. In addition, the contribution of burn severity to LSA was relatively low compared to its contribution to the rest of the variables, and the contribution of pre-fire vegetation was relatively high.

Table 5. Mean values of pre-fire (06/15/17) ET, LST, LSA, NDVI and NBR images in different classes taken into account for climate, vegetation species, elevation, slope and aspect.

		ET (mm/day)	LST (K)	LSA	NDVI	NBR
Climate	Csa	4.85	307.44	0.13	0.62	0.44
	Csb	5.31	306.64	0.13	0.63	0.44
Vegetation	Eucalyptus	4.92	307.47	0.12	0.66	0.49
	Pine	4.82	308.03	0.11	0.67	0.49
	Shrub	4.74	307.98	0.12	0.65	0.43
Elevation (m)	0-300	6.34	305.82	0.14	0.64	0.46
	301-600	5.35	307.51	0.13	0.66	0.50
	601-900	4.94	306.61	0.13	0.64	0.47
	>900	3.69	308.21	0.13	0.56	0.33
Slope (º)	0-5	5.63	305.57	0.12	0.65	0.46
	5-10	5.20	306.72	0.13	0.63	0.44
	10-20	4.87	307.59	0.13	0.62	0.43
	20-30	4.88	307.57	0.14	0.62	0.44
	>30	4.82	307.74	0.14	0.61	0.42
Aspect	North	5.42	305.83	0.13	0.66	0.48
	Northeast	4.97	307.33	0.14	0.64	0.45
	East	4.73	308.16	0.14	0.63	0.43
	Southeast	4.70	308.44	0.14	0.60	0.40
	South	4.92	307.68	0.13	0.61	0.42
	Southwest	5.03	307.20	0.13	0.60	0.42
	West	5.38	306.12	0.13	0.62	0.44
	Northwest	5.49	305.54	0.13	0.64	0.46

summarizes the mean values of ET, LST, LSA, NDVI and NBR before the fire for every class, defined in climate, vegetation species and topographic factors. As expected, higher mean values of ET, NDVI and NBR (and lower of LST) were found in the humid Mediterranean climate (Csb class), at lower elevation (0–300m, 300–600m classes) and at lower slope (0°–5°, 5°–10°) ranges because of a high water availability due to the close proximity to the Zêzere River, and on North aspect slopes (more humid). Regarding vegetation species, differences in mean values of the variables were smaller.

4.3. Is it Possible to Estimate Accurately Burn Severity From the ET, LST and LSA Images?

A summary of the accuracy measures of the burn severity estimates from each variable taken into account (uni-temporal and multi-temporal perspective) is displayed in

Table 6. Accuracy summary of the burn severity estimates from ET, LST, LSAT, NDVI and NBR (uni-temporal and multi-temporal perspective).

	Uni-Temporal Perspective (07/01/2017)					Multi-Temporal Perspective (06/15/17–07/01/2017)
	ET	LST	LSA	NDVI	NBR	dET_i	dLST_i	dLSA_i	dNDVI_i	dNBR_i
κ	0.63	0.57	0.44	0.59	0.61	0.55	0.52	0.45	0.65	0.66
PA	0.69	0.65	0.49	0.67	0.70	0.66	0.63	0.58	0.70	0.71
UA	0.68	0.66	0.51	0.70	0.69	0.62	0.60	0.55	0.72	0.70
OA	0.73	0.69	0.61	0.72	0.72	0.67	0.64	0.60	0.75	0.76

. From a uni-temporal perspective, all burn severity estimates except the LSA-based one showed adequate accuracy. The post-fire ET-based estimate achieves the highest κ statistic (0.63), followed by the post-fire the NBR-based estimate (κ = 0.61). When a multi-temporal perspective was adopted, both burn severity estimates based on dET and dLST decreased their accuracy slightly, whereas burn severity estimates based on dNDVI and dNBR increased their accuracy slightly. Consequently, some differences in accuracy were observed in this case: κ statistic was 0.55 in dET-based estimate, and 0.51 in dLST-based, compared to 0.65 and 0.66 in dNDVI- and dNBR-based ones, respectively. The dLSA-based burn severity estimate showed approximately the same accuracy as the LSA-based one in the uni-temporal perspective (too low to be acceptable).

Figure 3 provides a visual comparison of burn severity (initial assessment) maps based on immediately post-fire ET and dNBR_i (used as a reference). From it, we can verify a high level of similarity between the two burn severity maps. A small variation in burn severity level actually associated to aspect changes may be observed in the upper right corner of the ET-based map. We should remember that this aspect contributed 6.20% to the variance in immediately after fire ET image (Table 4).

5. Discussion

5.1. How Does Fire Modify ET, LST, LSA, NDVI and NBR?

Our results agree with previous studies showing that ET decreases after a fire event (e.g. [9,89] in semi-arid; [15,90] in boreal forests; and [13,35,91] in temperate conifer forests). In Mediterranean ecosystem, Sánchez et al. [17] studied the temporal evolution of changes insurface energy balance over an eleven-year period following fire in a maritime pine-shrub mixed forest (Spain). They observed that, after 6 to 7 years, energy fluxes recovered to their pre-fire values in shrubs, while maritime pine forest needed 9 years. Häusler et al. [37], in a forested area near our study area, observed that, similar to maritime pine forest, ET decreased after fire in the Mediterranean Eucalyptus dominated forests, in particular at moderate and high burn severity levels. In eucalypt forests, but in Australia, Nolan et al. [92] found that ET decreased after fire and was higher in forests burned at high severity compared to moderate levels.

Patterns of variation due to fire in LST were similar to ET but inversed. The increase in LST was higher at high burn severity levels than in moderate and low ones. This increase was also observed one year after the fire. These results are consistent with previous studies concluding that LST shows an increase immediately after wildfire that may persist even years after the event [19,21,93,94,95,96,97]. Liu et al. [25] observed that burned forested areas showed a higher LST mean and variability than unburned areas in boreal forest. In our study, variability was higher in burned areas than in unburned ones in only the first year after the fire. In contrast, variance in unburned areas in 2018 exceeded the variability of burned areas. The 2018 decrease in LST variance in burned areas suggests a reduction in the number of burn severity levels that may be differentiated between, as burn severity determines the magnitude of post-fire changes in LST [90].

Immediately after a fire, ash and charcoal accumulation over soil and charred surfaces causes a decrease in LSA [98,99] whose amplitude is influenced by burn severity [22]. However, both 2018 LSA images displayed higher values for burned areas than for unburned ones. The principal reason for this increase is strong regeneration of Mediterranean vegetation following a fire event [100]. Shrubs begin to regrow in a few weeks after the fire event [101]. E. globulus, one of the dominant species in our study area, has a high resprouting capacity, mainly due to a high epicormic bud initiation potential and to the presence of a lignotuber [102]. By contrast, P. pinaster, the other dominant species, is a seeder, with a rapid seed dispersal from the beginning of fire until some months after its extinction [103]. LSA values in the summer 2018 summer image were also observed to be higher than those of the fall image, suggesting that LSA is greatly influenced by seasonality [20]. Ecosystem type and fire regimes, however, also have an important influence on LSA changes in duration, sign and magnitude [25,99].

Mean NDVI and NBR values in burned and unburned surfaces displayed a similar trend: a fall immediately after the fire, whose amplitude followed a decreasing pattern over time. This pattern responds to the changes in radiation recorded by satellite sensors due to both the removal of healthy vegetation and the presence of ash and charcoal [104]. Healthy vegetation absorbs red radiation and reflects NIR radiation. Conversely, burned vegetation reflects more radiation in red and SWIR wavelengths and absorbs NIR radiation, making their spectral signature flatter than the signature of healthy vegetation [6]. However, NDVI is not the most appropriate index to map burned areas in Mediterranean-type regions [105], as it has not been designed to discriminate ash and char, predominant materials in areas recently affected by fire [106]. The main advantage of NDVI is that it can be computed from almost all existing sensors [107]. Although thresholding the differenced form of NBR (dNBR) [6] has become a standard method to estimate burn severity due to its easy computation [60], some limitations of this method have been pointed out (see [66]) and new alternatives have been proposed. Among them, the ones based on char fraction from unmixing reflective bands stand out [66,75,108,109,110].

5.2. How Are ET, LST, LSA, NDVI and NBR Post-Fire Images Influenced by Burn Severity and Pre-Fire Factors?

We clearly observed the dependence of ET, LST, NDVI and NBR on burn severity in our study (Table 4). Fisher’s least significant difference showed that, from an uni-temporal point of view, the three energy balance related variables and the two spectral indices discriminated three burn severity levels with statistical significance the year of the fire event. Conversely, 1-year after fire, only ET (summer and fall images) and LST (summer image) enabled such discrimination (Table 2). From a multi-temporal approach, results similar to the uni-temporal procedure were obtained. Although energy balance variables showed a strong dependence on burn severity, pre-fire vegetation species, and especially topographic factors, also had a relatively important contribution. Climate (Köppen-Geiger classes, [72]) did not contribute to explain the variance in any of the variables used in our study, except ET in the extended assessment. Our study area is a transition zone between Mediterranean (Csa) and humid Mediterranean (Csb), and there is no clear border between these two climatic classes. We believe that the climatic fuzziness resulted in the climate having almost no influence on the five studied variables.

Vegetation species type has an important influence on energy and water exchange between atmosphere and land [25,89,111]. However, in our study vegetation species had a low to moderate contribution to ET, LST, NDVI and NBR in the initial assessment, and a higher contribution to the extended one. We found very similar ET and LSA values in the main pre-fire vegetation species, although broadleaf trees usually show higher values in ET and LSA than coniferous trees [25]. However, Eucalyptus cannot be =considered a broadleaf tree, and thus its ET values are close to the ET values found in pines (Table 5). With regards to topographic factors, they had a reasonable contribution to ET and LSA, and a lower one to the spectral indices. Aspect was the most influential topographic factor in the initial assessment, with elevation together with aspect forming the most relevant topographic factors in the extended one. Previous studies [92,112] have already demonstrated the important effect of topography on ET, through its influence on evaporative demand and forest structure. Furthermore, topography also influences the recovery pattern ET following the fire event, as vegetation regrow relies partly on elevation [91].

5.3. Is it Possible to Accurately Estimate Burn Severity From the ET, LST and LSA Images?

From a uni-temporal perspective, burn severity estimates based on ET achieved the highest accuracy, with (κ = 0.63) being the acceptable accuracy of LST-based estimates (κ = 0.57). From a multi-temporal perspective, both dET_i- and dLST_i-based estimates showed a lower accuracy than dNDVI_i- and dNBR_i-based ones (multi-temporal approach) (Table 6). The obtained accuracy values were in accordance with accuracy values found in previous studies (e.g., [18,60,75,113,114,115,116,117]).

In our study, burn severity estimates based on variables related to energy balance reached similar accuracies to estimates based on commonly used spectral indices. However, according to Parks et al. [118], these variables should be included to model burn severity across broad regions, due to the difficulty in interpreting models based on spectral indices [64,118,119]. Moreover, these variables provide critical information about the alterations in energy balance caused by fire, and in particular, information on water balance. Fire has a clear impact on stand water cycle [120] and hydrology at watershed scale [121]. Forest wildfires, importantly, alter the hydrological processes, which may impact on the growth of the forest and climate at regional scale [122,123].

5.4. Final Considerations and Future Work

Most previous studies on the impact of fire on energy balance have been based on field-measured surface fluxes, and are only valid at very local scales [13,15,91,92,124]. Others rely on variables derived from coarse spatial resolution satellites (e.g., MODIS), and are only valid at global scales [90,125,126]. There are very few studies (e.g., [17,37]) that use satellite-data-derived variables at a fine spatial resolution. The availability of these data makes it possible not only to work at a local-regional scale, but also to study the potential of these energy balance variables as indicators of burn severity. In our study, ANOVA analysis showed that ET, LST and LSA could discriminate three burn severity levels with statistical significance at initial assessment (both from uni- and multi-temporal perspectives). Burn severity estimates from post-fire ET not only reached almost the same accuracy as the widely used method proposed by Key and Benson [6], but also highlighted the important impact of fire on water balance. Since forests absorb and keep rainfall to make it available during dry season [127], wildfires alter the volume stored by forests by removing the vegetation cover [128]. The consequences of that change on water balance may include: modifications of ET from burned soil and vegetation versus unburned ones [37]; alterations in soil moisture [129]; an increase in the volume of precipitation [130]; a post-fire increase in soil erosion (due to higher surface runoff) [131]; lower water quality caused by an increase in suspended sediments [132] and modifications to surface hydrology [128]. Furthermore, unlike spectral indices, ET has a physical meaning, which makes it easier to interpret than spectral indices. Moreover, procedures based on variables related to energy balance such as ET or LST, estimated from fine-scale remotely sensed data, might reduce over-dependence on spectral indices to infer burn severity [65].

Thus, future research studies about the relationship between burn severity and satellite-derived variables related to water balance should be focused on three areas. First, the classification procedure of the variables related to energy balance needs to be revisited, and other classification algorithms should be tested. Second, future studies should evaluate the integration of post-fire energy balance measures with widely used spectral information to increase the accuracy of burn severity estimates. Third, it is necessary to verify the potential use of these variables as indicators of burn severity in other ecosystems and fire regimes.

6. Conclusions

Wildfires not only modify surface reflectance but also induce changes in energy balance. Our study proved the usefulness of Landsat-derived variables related to energy balance to estimate burn severity with a similar accuracy to other well-established satellite-based methods, with the added advantage of providing a better understanding of fire-induced changes in energy balance. We analyzed the post-fire temporal changes in energy fluxes by means of three variables—ET, LST and LSA—within the area affected by a mega-fire in Mediterranean forests, contributing to the evaluation of their important ecological consequences. Although the relationship between these variables and burn severity has been previously recognized, this is the first study that showed their potential as indicators of burn severity and that used them to estimate accurately burned severity.

Our results revealed that fire has a strong impact on energy balance variables at a regional-local scale during the first year after fire: ET decreased and LST and LSA increased. We showed that the fire-induced changes in these variables were related to burn severity. Three burn severity levels could be distinguished with statistical significance from ET and LST images using uni-temporal and multi-temporal approaches. Thus, the proven relationship between ET and fire damage may provide a physical meaning to burn severity in terms of changes in water balance.

We mapped burn severity from each potential indicator using thresholds based on the mean values for each class. Accuracy of burn severity estimated using ET was the highest from an uni-temporal perspective, although the dNBR-based burn severity map had the highest accuracy from a multitemporal one. We conclude that the three Landsat-derived energy balance variables analyzed (especially ET and LST) can be seen as effective indicators of burn severity in mega-fires in Mediterranean ecosystems. In addition, these variables in themselves provide important measures of the impact of fire on energy balance. The proposed method for estimating burn severity may help advance burn severity assessment and our understanding of burn severity patterns. Finally, it takes advantage of new satellites with a similar spatial resolution to Landsat, equipped to provide an accurate estimate of daily ET, such as ECOSTRESS.