Vegetation and Soil Fire Damage Analysis Based on Species Distribution Modeling Trained with Multispectral Satellite Data

Forest managers demand reliable tools to evaluate post-fire vegetation and soil damage. In this study, we quantify wildfire damage to vegetation and soil based on the analysis of burn severity, using multitemporal and multispectral satellite data and species distribution models, particularly maximum entropy (MaxEnt). We studied a mega-wildfire (9000 ha burned) in North-Western Spain, which occurred from 21 to 27 August 2017. Burn severity was measured in the field using the composite burn index (CBI). Burn severity of vegetation and soil layers (CBIveg and CBIsoil) was also differentiated. MaxEnt provided the relative contribution of each pre-fire and post-fire input variable on low, moderate and high burn severity levels, as well as on all severity levels combined (burned area). In addition, it built continuous suitability surfaces from which the burned surface area and burn severity maps were built. The burned area map achieved a high accuracy level (κ = 0.85), but slightly lower accuracy when differentiating the three burn severity classes (κ = 0.81). When the burn severity map was validated using field CBIveg and CBIsoil values we reached lower κ statistic values (0.76 and 0.63, respectively). This study revealed the effectiveness of the proposed multi-temporal MaxEnt based method to map fire damage accurately in Mediterranean ecosystems, providing key information to forest managers.


Introduction
Wildfires have become a major concern in recent years, particularly in the Euro-Mediterranean region [1], with significant ecological consequences on ecosystems [2].Into this region, Spain represents one of the most affected countries, where approximately 180,000 ha were burned during 2017, all in forested areas.There has been an increase over the past decade in catastrophic episodes and in large forest fires (>500 ha).Large fires represented almost 55% of the total burned area in Spain, primarily in the Northwest geographical region (approximately 75% of the surface burned by large fires occurred in this region) [3].The effects of fire on ecosystems are mediated, among other factors, by the fire regime parameters.In this sense, burn severity is one of the most significant and might alter the vegetation recovery capacity [4] and soil properties, issues that must be critical for forest management.Additionally, wildfire effects on vegetation depend on plant functional traits [5].For example, regeneration tends to be more successful for resprouting species [4], but in areas dominated by obligate seeders it could fail, resulting in a reduction of biomass and vegetation cover, modification of landscape structure and shifts in vegetation community composition [6,7].At the same time, burned areas exhibit physical, chemical, mineralogical and biological soil alterations due to high temperatures and charcoal-ash deposition [8,9].Thus, fire modifies soil properties [10], affecting the mineral soil condition of forests [11], and increasing electrical conductivity and pH [12].Both the high temperatures reached during a fire and the increase in soil pH lead to important variations in soil nutrients that are key for post-fire regeneration of vegetation and soil microbiota [13].Loss of vegetation also implies high soil vulnerability to erosion by wind and/or water [6,14].Indeed, these soil changes may influence vegetation structure and/or composition as well as driving vegetation dynamics over time [15,16].
Fire effects may be mitigated by post-fire policies focused on vegetation regeneration and soil preservation.An accurate and rapid monitoring and assessment of fire effects provides key information for designing appropriated management strategies [17,18].Both burned surface area and burn severity are the two most commonly utilized metrics for evaluating fire effects [19].Burn severity or vegetation burn severity refers to the vegetation after fire comprising both short-and long-term impacts (weeks after fire and one year after fire, respectively) [20].With regards to soil burn severity is defined as the loss of organic matter in soil or on its surface [21] and it is a decisive factor for resolving hillslope stabilization and erosion reduction [22].At the same time, vegetation burn severity maps help forest managers to determine the burned area that needs ad hoc management strategies to limit soil and vegetation degeneration [18].This suggests a need for integrated evaluation of the potential damage of fires on ecosystems considering the large surface affected.Remote sensing provides key information for wildfire damage assessment in the short, medium and long-term [6].Examples of its relevance are official programs such as the Monitoring Trends in Burn Severity (MTBS, https://www.mtbs.gov)Project [23] or the Mapping Burn Severity for Burned Area Emergency Response (BAER, https://www.fs.fed.us/eng/rsac/baer/),both in the USA, the European Forest Fire Information System (EFFIS, http://effis.jrc.ec.europa.eu)or the Australian Landgate's Satellite Remote Sensing Services (SRSS, http://srss.landgate.wa.gov.au/fire.php).
Most remote sensing based burn severity/burned area research is based on a classification of a preand post-fire difference image, or of a single post-fire image [24].Different classification algorithms have been used: Threshold [25,26]; support vector machines [27]; object-based image classification [28,29] or neural networks [30,31].We selected the maximum entropy one-class classifier (MaxEnt) [32,33] to model and map burn severity and burned area.The main advantages of MaxEnt include: (1) It has been designed to work with presence-only data; thus, it constitutes an interesting alternative to other machine learning based classifiers [34]; (2) its probabilistic output is easy to interpret and has physical meaning [35] and (3) it is a non-parametric model (the input variables interrelations are not determine a priori) [36].MaxEnt is widely used in ecological studies to model species distributions [37][38][39] and it is increasingly being used in remote sensing based applications like pest potential distribution [40], landslide vulnerability monitoring [41], groundwater potential distribution [42] or fire occurrence modeling [43][44][45][46], so the use of MaxEnt to deal with the effects of burn severity encourages us to check it in the present study.MaxEnt provides the percentage of contribution tables and probability distributions of each target class (three burn severity levels and burned area, in our study), from which burn severity and burned area maps can be built.Thus, we will be able to not only map fire damage but also to know the relative influence of each variable used as input into the algorithm.As independent input variables for MaxEnt, we used fraction images from pre-and post-fire multispectral data and a post-fire land surface temperature (LST) image.
Fraction images were used instead of conventional spectral indices due to their proven suitability in fire damage studies [18,[47][48][49].Their success is related to the characteristics of the short-term post-fire domain that usually produces a mix of remnants of vegetation, ash and burned soil, essentially constituting a sub-pixel question at the spatial resolution of commonly used multispectral sensors [50].In addition, most spectral indices are based on just two spectral bands, neglecting the information included using all bands, whereas fraction images are typically computed using all spectral bands [51].
Another important advantage of fraction images over spectral indices is that they do not need to be calibrated with field data [52].Fraction images have a physical meaning; they estimate the distribution of the different cover classes present in a pixel, which makes them easier to interpret [53].The technique most commonly used to unmix the original image is the spectral mixture analysis (SMA) [54].SMA, however, does not consider within-class spectral variability in the scene as it defines just one endmember spectrum for each class.Conversely, multiple endmember SMA (MESMA) [55] allows multiple endmember spectra in each class, which enables it to overcome the SMA limitation.Fraction images from MESMA have shown their value in a large variety of remote sensing applications [47,[56][57][58][59][60] and in particular in post-fire damage assessment [48,50,52,[61][62][63].Moreover, they have proven to be highly correlated to burn severity field measurements (frequently using the composite burn index, CBI; e.g., [48,49]).
Post-fire LST was also incorporated in our study as an input variable for the MaxEnt classifier.The main reasons for this decision are: (1) Post-fire LST can be regarded as a useful potential indicator of vegetation burn severity [63][64][65][66].Furthermore, although information comprised in thermal infrared bands is usually discarded in fire damage studies, soils after fire exhibit little to no evaporation that may be expressed as high temperatures in the thermal [18]; and (2) a clear relationship between soil temperature and variations in soil characteristics suggests a potential use of LST as the soil burn-severity indicator as shown by previous studies [67][68][69].Additionally, LST has been linked to soil burn severity as LST reveals variations in soil characteristics, principally in organic carbon content that is the main driver of the post-fire physical changes [21].
Within this framework, our study proposes a method to analyze post-fire vegetation and soil damage based on multispectral remotely sensed data in Mediterranean forestry ecosystems.Our main goals are to evaluate the relative contribution of each pre-and post-fire input variable on the three burn severity levels and to map accurately burned area and burn severity in three burn severity scenarios: Total burn severity, soil and vegetation burn severities.From a practical point of view, our aim was to provide some basis for managing forests to reduce fire damage effects.In particular, our specific research questions were: (1) Could burned area, and in particular each of the three considered burn severity levels, be accurately modeled by a multitemporal approach based on MaxEnt trained with preand post-fire fraction images and post-fire LST?What would be the relative influence on the study of fire damage of pre-fire information versus post-fire one?Would pre-and post-fire variables have the same influence on each burn severity level?and (2) how is the MaxEnt-based burn severity map (total burn severity) related to field measured vegetation and soil burn severity?

Study Area
We used as a study case a large wildfire (98.8 km 2 ) that occurred from 21 to 27 August 2017 in the region of La Cabrera, (León, Northwest Spain; Figure 1).It was the largest fire in Spain during 2017 and one of the most significant wildfires in size and damage located in European Union [3].The main environmental conditions that caused the rapid spread of the fire were: Maximum temperatures of 35 • C, low relative air humidity (35%) and most importantly 53 days without any precipitation before the fire, so fuel had a high probability of ignition.
Burned area was limited in the north by "Sierra del Teleno" and "Montes Aquilianos" mountain, to the south by "Sierra de la Cabrera" mountains, and to the west by the mountains of "La Mina".Relief is rugged as can be observed from a topographic map (Figure 1b).The elevation of burned area ranged from 836 m to 1958 m (average altitude = 1326.9m). Figure 1c) shows the slope map of the study area.Slope varied from 0 • to 55.07 • , with an average slope of 19.7 • .A south facing aspect is dominant (Figure 1d).Our study area had a Csa-Mediterranean climate (warm summer) according to the Köppen climate classification [70].The area has cold winters, with frost and precipitation falls as a mixture of water and snow.Summers are warm, exceeding 21 • C average with a maximum temperature of 30 • C and pronounced drought [70].The average annual temperature in the burned area is 8.28 • C (ranging from 4.57 • C to 10.77 • C, Figure 1e).The relief greatly influences the annual rainfall whose average is 1129.12mm/year but oscillates between 639.60 and 1558.2 mm/year (Figure 1f).These characteristics, together with acidic and poor soils have generated a predominance of shrub vegetation, where forests are mainly found in small spots dotted throughout the area and in low riverbank areas.According to the Spanish forest map (Third Spanish National Forest Inventory [71]) the burned area was covered by shrubs (65%), dominated by Erica australis L. heathlands and Genista hystrix broomlands, natural oak woodlands (15%) dominated by Quercus pyrenaica Willd.and Quercus ilex L., forest plantations (6%) of Pinus pinaster Ait. and Pinus sylvestris L., pasturelands (10%) and crops (4%).The fire affected 54.4 km 2 of shrubs, 34.8 km 2 of woodland, 5.8 km 2 of pasture and 3 km 2 of non-forest land.

Field Data
Burn severity was measured in the field during the first three weeks after the wildfire.We established a total of 89 circular field plots (30 m diameter) distributed between burned areas (14 in low burn severity, 19 in moderate and 30 in high) and control (26).These field plots were distributed among the three main ecosystems of the study area: 33 on heathlands (Erica australis); 25 on broomlands (Genista hystrix) and 31 oak woodlands (Quercus pyrenaica), in areas with similar vegetation structure and fire effects.A high-precision Global Positioning System (GPS) receptor (Spectra Precision Mobile Mapper 20) georeferenced the center of each plot (Root-mean-square error -RMSE-< 0.5 m).
The field protocol to quantify CBI was an adaptation of the first proposed CBI procedure [72] (see [2] Table 2 for a detailed description).Basically, different parameters related to burn severity were visually rated on a scale of 0 (unburned) to 3 (highest burn severity) in five different strata (substrate, vegetation <1 m, vegetation 1-5 m, vegetation 5-20 m and vegetation >20 m), and the average value per stratum was computed.Next, burn severity was calculated by averaging the five obtained values (CBI).If we average only the values corresponding to all strata with the exception of substrate, we obtain the vegetation burn severity (CBI veg ), and finally, the value of substrate stratum represents soil burn severity (CBI soil ; for more details see [2]).CBI thresholds for the three considered burn severity classes were defined according to [73]. Figure 2 displays picture examples of field plots for high, moderate and low burn severity level located in the three more representative habitats.

Remotely Sensed Data
Pre-and post-fire MESMA fractions were calculated from Sentinel-2 Multispectral Instrument (MSI) data.The MSI sensor records information in almost the same wavelengths as Landsat 8 Operational Land Imager (OLI) with a finer spatial resolution (10 m for visible and near-infrared bands and 20 m for near-infrared and short wave near-infrared bands).It includes, however, the red-edge domain that has proven useful for detecting and assessing forest disturbances caused by fires [74].The pre-and post-images were downloaded from Copernicus Open access Hub (https://scihub.copernicus.eu)with a processing level 2A: Bottom-of-atmosphere reflectance in cartographic geometry [75].For level 2A, tiles are 100 × 100 km 2 ortho-images in UTM/WGS84 projection.As our area of study was included in two tiles, we needed to mosaic them as a prior step to subsetting to the selected forest fire.The closest images to the fire date without clouds were images acquired 13 August 2017 (pre-fire), and 2 September 2017 (post-fire).An RGB 11:8a:12 color composite from post-fire Sentinel-2 MSI image is displayed in Figure 1.The red-edge spectral bands were rescaled to a spatial resolution of 20 m.
In addition, we downloaded a post-fire Landsat 8 OLI/ Thermal Infrared Sensor (TIRS) scene acquired on 14 October 2017 from the USGS Earth Explorer server (http://earthexplorer.usgs.gov/) to compute post-fire LST.As for the Sentinel-2 MSI data, we searched for imagery with no or very few clouds (no clouds over our study area window), close to the fire date.The level of processing of the scene was L1T: Radiometrically and geometrically corrected image [76].

Additional Data
The Spanish National Center of Geographic Information (CNIG; http://www.cnig.es/)through the Spanish Aerial Orthophotography National Planning (PNOA) agency provided us with a 0.25 m-pixel digital orthophotograph acquired in 2015.Both Sentinel-2A and Landsat 8 data were first coregistered with the orthophotograph.The orthophotograph also assisted us in defining the candidate endmember spectra for the MESMA procedure.Additionally, the 2017 official land use maps of Castilla and León (raster of 10 m of spatial resolution) together with the Spanish forest map (1:50000) were used to digitize polygons over the orthophotograph from which candidate endmember spectra were obtained.A Copernicus Emergency Management Service (EMS) grading map (ID: EMSR227, four burn severity levels) also helped in the selection of candidate endmember spectra for burned areas.Presence-only samples for training the MaxEnt classifier to model each burn severity class were extracted as well from the Copernicus EMS-grading.Training samples to model burned class were located on Copernicus EMS-burnt area map (ID: EMSR227).Finally, a 25 m digital elevation model (DEM) supplied also by the Spanish PNOA was employed to retrieve LST from the Landsat 8 OLI image, together to the Moderate-Resolution Imaging Spectroradiometer (MODIS) water vapor product (MOD05) and global surface summary of the day (GSOD) from National Oceanic and Atmospheric Administration (NOAA), that provided meteorological information regarding visibility.

Methods
The proposed method to evaluate fire damage has four main stages: (1) MESMA procedure.In this first step fraction images were obtained from Sentinel -2 MSI data; (2) LST calculation.A LST image was calculated from Landsat 8 OLI/TIRS data; (3) MaxEnt modeling.We took into account four target-classes (low, moderate and high burn severity levels, plus all burn severity levels combined-burned class).This step includes analysis of different covariates contribution to the model, and model evaluation; and (4) Burn severity and burned area mapping.The optimal thresholds to classify the continuous suitability surfaces produced by MaxEnt were identified and assessed by using CBI values (Figure 3).We organized the selected endmember spectra into three spectral libraries.To unmix the pre-fire image we used green vegetation (GV), non-photosynthetic vegetation (NPV) and soil spectra libraries.Thus, we will be able to analyze their potential contributions to the final MaxEnt models of burn damage.However, to unmix the post-fire image, the char spectral library was included, and we grouped in just one spectral library (NPVS) the endmembers related to NPV and soil as Quintano et al. [62] suggested.Shade was included as the fourth spectral library.By using these libraries we unmixed the data and obtained the fraction images.Constraint values similar to ones suggested by [82] were employed: The maximum and minimum admissible fraction values were set to 1.10 and -0.10 respectively; the maximum allowable shade fraction value set at 0.85 and the maximum allowed RMSE was fixed at 0.025.The unmixing process was repeatedly executed by varying the type and/or number of the spectra incorporated in the spectral libraries until the number of classified pixels was above a threshold, 95% in our case (Figure 3).Finally, we removed the contribution of the shade endmember by shade normalizing the MESMA fraction images, thus emphasizing the contribution of non-shade endmembers.All operations to compute the fraction images were accomplished by the Visualization and Image processing for Environmental Research (VIPER) tools software [81] developed at the Department of Geography at University of California Santa Barbara (https://sites.google.com/site/ucsbviperlab/viper-tools).

LST calculation
LST calculation required the processing of both reflective and thermal Landsat 8 OLI/TIRS bands (Figure 3), as land surface emissivity (; needed to compute LST) was estimated from reflective bands as suggested by Sobrino et al. [83].
Digital numbers of reflective bands were firstly converted to radiance values (L) and atmospherically corrected using the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) algorithm [84].The input parameters required by FLASH (atmosphere models and

MESMA Procedure
We applied MESMA [55] to calculate per-pixel fractional cover of different ground cover classes in the two Sentinel-2 MSI images.As it was done in Quintano et al. [50,62] and Fernández-Manso et al. [56] image endmembers were used to form the candidate spectral library.They are selected from the original data and thus they have the same scale; additionally they can be obtained reasonably easily [77].According to Roth et al. [59], georeferenced polygons were delineated in the fine spatial resolution orthophotographs to provide the candidate endmembers to our multitemporal spectral library.The polygons were defined over just one uniform class with the help of the 2017 official land-use maps of Castilla and León, the Spanish forest map, and the Copernicus EMS grading map.To create the definitive spectral library, the most appropriate endmembers were chosen from the candidate endmembers taken into account three parameters [52]: (1) Minimum average spectral angle (MASA) [78], which includes the endmembers with the lowest average spectral angle; (2) endmember average RMSE (EAR) [79], which chooses the endmembers that yield the lowest RMSE for each class and (3) countbased endmember selection (CoB) [80], which selects the endmembers modeling the largest number of endmembers within their class.Furthermore, as recommended by Roberts [81] our knowledge of the study area as well as of the common spectral signature of the selected spectra helped us in the selection process.These two steps (definition of the candidate endmembers and selection of the definitive endmembers) are key to the success in the following unmixing stage, as they allowed us to identify the most adequate endmembers.
We organized the selected endmember spectra into three spectral libraries.To unmix the pre-fire image we used green vegetation (GV), non-photosynthetic vegetation (NPV) and soil spectra libraries.Thus, we will be able to analyze their potential contributions to the final MaxEnt models of burn damage.However, to unmix the post-fire image, the char spectral library was included, and we grouped in just one spectral library (NPVS) the endmembers related to NPV and soil as Quintano et al. [62] suggested.Shade was included as the fourth spectral library.By using these libraries we unmixed the data and obtained the fraction images.Constraint values similar to ones suggested by [82] were employed: The maximum and minimum admissible fraction values were set to 1.10 and −0.10 respectively; the maximum allowable shade fraction value set at 0.85 and the maximum allowed RMSE was fixed at 0.025.The unmixing process was repeatedly executed by varying the type and/or number of the spectra incorporated in the spectral libraries until the number of classified pixels was above a threshold, 95% in our case (Figure 3).Finally, we removed the contribution of the shade endmember by shade normalizing the MESMA fraction images, thus emphasizing the contribution of non-shade endmembers.All operations to compute the fraction images were accomplished by the Visualization and Image processing for Environmental Research (VIPER) tools software [81] developed at the Department of Geography at University of California Santa Barbara (https://sites.google.com/site/ucsbviperlab/viper-tools).

LST Calculation
LST calculation required the processing of both reflective and thermal Landsat 8 OLI/TIRS bands (Figure 3), as land surface emissivity (ε; needed to compute LST) was estimated from reflective bands as suggested by Sobrino et al. [83].
Digital numbers of reflective bands were firstly converted to radiance values (Lλ) and atmospherically corrected using the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) algorithm [84].The input parameters required by FLASH (atmosphere models and aerosol situations) were defined from MOD05, GSOD from NOAA and DEM from Spanish PNOA.Secondly, the C-correction algorithm [85] allowed us to remove the topographic shadow effects.C-constant was defined by using 10% of pixels [64].Next, land surface reflectance in percentage (ρ) was calculated from the topographic corrected values.Finally, from the reflectance values we computed ε based on NDVI thresholds, red band reflectivity and proportion of vegetation cover (see [83]).
Regarding thermal band processing, LST was calculated using the Landsat 8 OLI /TIRS band 10 (10.5 to 11.2 µm) [2], and the radiative transfer method (a single channel method, see [86].This procedure comprises three stages: (1) Radiometric calibration to convert Digital Numbers (DNs) to radiance values; (2) calculation of ground radiance, including atmospheric correction and ε adjustment and (3) conversion of ground radiance to LST.The radiative transfer equation expresses ground radiance as a function of radiance of the thermal band, atmospheric parameters and ε [86].For thermal band 10 it can be written as: where B(T s ) represents the ground radiance; B(T), the radiance received by the sensor (band 10) with brightness temperature T (from radiometric calibration of DNs); τ(θ) symbolizes the atmospheric transmittance for a view zenith angle θ; I↑ and I↓, respectively, denote the upwelling and downwelling path radiance and ε, the land surface emissivity.In our study, I↑, I↓ and τ(θ) were given by the National Centers for Environmental Prediction (NCEP) profiles [87], ATMCORR tool and, as mentioned above, ε, was calculated following the procedure recommended by Sobrino et al. [83] based on reflective bands.Finally, according to Planck's law, ground radiance was transformed into LST (Kelvin) using Equation ( 2): where K 1 and K 2 are found on Landsat 8 metadata and represent thermal constants.The LST image was rescaled to a spatial resolution of 20 m similar to Sentinel-2 MSI fraction images.

MaxEnt Modeling
Once the covariates for MaxEnt were ready, MaxEnt software (version 3.4.1),developed by [32,33] and updated by [88], was used in our study for modeling independently four target classes: Low burn severity, moderate burn severity, high burn severity and burned area.MaxEnt requires two kinds of input information: (1) Georeferenced locations of presence-only training samples.Training samples for low, moderate and high burn severity target classes were obtained by randomly sampling the considered target class from the Copernicus EMS grading map; training samples for burned class from the Copernicus EMS burnt area map.We used 151 samples of low burn severity, 133 of moderate burn severity, 142 of high burn severity and 126 of burned area; and (2) a set of independent variables (called covariates) that presumably have influence on the target class.We included as covariates the pre-and post-fire shade normalized fraction images and LST image.MaxEnt estimates a probability distribution of the target class by calculating the maximum entropy constrained to the partial knowledge of the target distribution [88].The MaxEnt software (from version 3.4.0onwards) provides by default a complementary log-log (clog-log) output, easy to conceptualize, that estimates the probability of presence of the target class [89].The spatial resolution of the output is the same that the spatial resolution of covariates (20 m in our study).In addition, MaxEnt provides several tools to analyze the contribution of each covariate on the final model.Specifically, the software builds a table reflecting the contribution of each covariate to the model, and three tables indicating the training gain, the test gain and the test AUC, when each covariate is removed and when it is considered alone.
Following [46] the random test percentage was set to 25%.Previous studies (e.g., [90,91]) have recommended adjusting the regularization multiplier value to maximize the fit of the model.We tried six regularization multiplier values (0.2, 0.5, 1, 2, 5 and 10) for each target-class.The regularization multiplier influences how focused the output distribution is.Thus, values lower than 1 (default value) produce output distributions more concentrated on the given presence samples, whereas values higher than 1 produce less localized predictions.The rest of the user-specified parameters were set to their default values as recommended by Dudik et al. [92] and Philips and Dudik [93]: Automatic establishment of the most appropriate complexity level according to the presence-only data sample size; number of background samples = 10000, convergence threshold = 5-10 and maximum iterations = 500.We used 10 replicates with repeated subsampling scheme to obtain more reliable results.
The MaxEnt software provides training and test receiver operating characteristic (ROC) analysis to assess the model, where the area under the curve (AUC) [94] represents the model's capability to predict correctly the presence (sensitivity) and absence (specificity).AUC can be seen as a widely accepted procedure to evaluate the accuracy of a predictive distribution [95].AUC values range from 0 to 1.An excellent model performance is suggested by AUC values higher than 0.9, a moderate good performance by values between 0.7 and 0.9, and a poor model performance by values from 0.5 to 0.7.A diagonal line (AUC = 0.5) represents a random prediction of presence and absence [96].Test AUC has been used to measure the model fit in many previous research works (e.g., [43,97]) because it is not affected by overfitting [98].In our study, we chose the optimal regularization multiplier as the one that maximized the test AUC.

Burned Area and Burn Severity Mapping
For each target class (low, moderate, high burn severities and burned area) MaxEnt produced a continuous output representing its suitability surface.To map the area affected by the fire, we transformed the continuous output for the burned area target class into a binary map by using a threshold.MaxEnt provides a summary of the thresholds that define the area more suitable for the target class maximizing sensitivity (presence) [33], among they: Fixed cumulative value 1, fixed cumulative value 5, fixed cumulative value 10, minimum training presence, 10 percentile training presence, equal training sensitivity and specificity; maximum training sensitivity plus specificity, equal test sensitivity and specificity or maximum test sensitivity plus specificity.In our study, the optimal threshold among the candidate thresholds provided by MaxEnt was chosen by maximizing the κ statistic [99] of the burned area map.An error matrix was built from CBI field-plots by grouping high, moderate and low burn severity classes into the burned class.To balance the number of burned and unburned samples, we added 49 unburned samples by randomly sampling the Copernicus EMS burned area map.
Similarly, to map burn severity from each of the three continuous outputs for the target classes high, moderate and low burn severity were converted into a binary images (unburned/considered burn severity class), and added to obtain a burn severity map.We confirmed empirically that prioritizing the high burn severity level allowed us to solve potential conflicts.As for burned area mapping, we found the optimal thresholds by maximizing the κ statistic of the final burn severity map.Finally, we assessed the relationship of the burn severity map to vegetation and soil burn severity measurements by computing error matrixes based respectively on CBI veg and CBI soil .

Results
The endmember selection process to form the definitive three spectral libraries used by MESMA is outlined in Table 1.We made use of the same set of endmembers covering soil, GV and NPV to unmix both pre-and post-fire Sentinel-2 MSI images.Figure 4 displays as example endmember spectra included in the final post-fire spectral libraries.Sentinel-2 data were successfully unmixed with a low percentage of unclassified pixels for each image (0.13% for post-fire image and 1.17% for pre-fire image; Table 1).The five covariates used as  The GV spectral library was formed by endmember spectra from two categories: Forest and shrub, subdivided into Quercus ilex, Populus sp., Quercus pyrenaica, Pinus sylvestris, Erica australis and Genista hystrix (Table 1).These six subcategories (hierarchical level 3) cover the main and most representative vegetation species of the study area.Regarding the soil spectral library, five soil subcategories were grouped into pervious and impervious surfaces (hierarchical level 2).We only had one category for representing NPV.Char endmember spectra were defined from pixels inside the fire perimeter with the help of the Copernicus EMS grading map.
Sentinel-2 data were successfully unmixed with a low percentage of unclassified pixels for each image (0.13% for post-fire image and 1.17% for pre-fire image; Table 1).The five covariates used as input for MaxEnt: Post-fire shade normalized char fraction image (post_char), post-fire LST (post_LST), pre-fire shade normalized GV fraction image (pre_GV), pre-fire shade normalized NPV fraction image (pre_NPV) and pre-fire shade normalized soil fraction image (pre_soil) are shown in Figure 5. From it, we can observe that pre_GV reveals clearly a previous burned area (blackish tone) that is also observed in pre_NPV (whitish levels); in particular, the wildfire happened in September 2016 and burned 1.208 ha.On the pre_soil image we could differentiate the whitish lines showing roads and firebreaks.Post_char allowed visual discrimination between burned and unburned areas.Moreover there are variations inside the fire perimeter suggesting different burn severity levels.Finally, on the post_LST image we could observe how temperature varies with altitude and aspect.

394
Moreover there are variations inside the fire perimeter suggesting different burn severity levels.

395
Finally, on the post_LST image we could observe how temperature varies with altitude and aspect.Interpretation of Table 2 allowed us to answer our research question 1.It shows the main parameters of the four MaxEnt modeling processes for the average model of the 10 replicates that were run in each process.For each target-class we tried different regularization multipliers from 0.2 to 10 and selected the one that maximized test AUC.We obtained test AUC values higher than 0.8 for the four target classes, indicating a good performance of the models.Test AUC for high burn severity class was the highest, almost reached 0.9, meaning an excellent model performance and low burn severity class presented the lowest one (0.83).Test AUC for burned area class was 0.86.These results enable us to affirm that burned area and burn severity levels can be accurately modeled by MaxEnt trained with post-fire LST and pre-and post-fire fraction images.From the percentage of contribution table (Table 2), post_char is the covariate that contributes the most in the four MaxEnt modeling processes (high, moderate, low burn severity and burned area).Its contribution is especially relevant to the burned area target class (approx.90%).Its contribution was also important when the high and moderate burn severity classes were considered the targets (respectively, approx.80% and 74%).Post_LST had a relatively important contribution for modeling the low burn severity class (approx.25%).Its contribution to the modeling of the moderate burn severity class was also noticeable (10%).Regarding the contribution of pre-fire covariates, pre_GV provided approximately 10% to model the high and low burn severity classes but its contribution to model the moderate burn severity and burned area classes was a bit lower (around 6%).The percentage of contribution of pre_NPV was almost the same for the three burn severity classes (around 8%) and a bit lower for the burned area class where post_char contributed almost exclusively.Finally, the contribution of pre_soil was small for the four target classes (negligible for high burn severity and for the burned area).
Summarizing, post_char presented the maximum contribution percentage to model the high burn severity class (approx.80%), making the contribution of both post_LST and pre_soil almost negligible.Training gain, test gain and test AUC corroborated this fact, and showed that post_char generalized better than the other covariates, making the model more transferable.When the low burn severity class is considered, the joint contribution to the model of both of the post-fire covariates (post_char and post_LST) almost reached 80%.Again, the training gain, test gain and test AUC confirmed this pattern and their generalization capacity.With regard to the moderate burn severity class, the analysis of the contribution to the model of each covariate is similar to the low burn severity class, adding the contribution to the model of post_char and post_LST a little more than 80%.Finally, if considering the burned area class, the post_char contribution to the good fit of the model had the maximum value (89.42%).Again the training gain, test gain and test AUC confirmed this fact.
Figures 6 and 7 show the MaxEnt continuous outputs of the four target classes for the best model (maximum test AUC) of the 10 replicates.They can be interpreted as the suitability surface for each target class occurrence, as they provide an estimation of the target class presence probability.It is notable that high and low burn severity suitability surface (Figure 6a,c) are reciprocal, and they complement well to form the total burned area.Regarding the continuous output for the burned area class (Figure 7a) we visually observed a very good match between the estimated and the real burned area, although there is a slightly overestimation over the previous burned area.Figures 6d and 7b display respectively the burn severity and burned area maps built from these continuous outputs.The accuracy of the burn severity map was tested by using the three CBI metrics (total CBI, CBIveg and CBIsoil; research question 2; Table 3).Regarding the CBI metric, the  statistic reached 0.  The accuracy of the burn severity map was tested by using the three CBI metrics (total CBI, CBI veg and CBI soil ; research question 2; Table 3).Regarding the CBI metric, the κ statistic reached 0.81.A relatively low user's accuracy (UA) value for the low burn severity class (0.64) could be observed, pointing to relatively high commission error.The moderate burn severity class presented modest producer's accuracy (PA) and UA values (0.76), suggesting moderate commission and omission errors.Unburned and high burn severity classes presented the highest PA (approx.0.90) and UA values (approx.0.95).When the CBI veg values were used, the accuracy parameters decreased slightly compared to the CBI metric, but the pattern was very similar.In this case, the κ statistic achieved 0.76 (Table 3).Again, unburned and high burn severity classes showed the highest PA (approx.0.90) and UA values (0.96 for unburned, and 0.82 for high burn severity), though the UA value for the high burn severity class was a bit lower than for the CBI metric, indicating some commission errors.As for the CBI metric, a relatively low UA value for the low burn severity class (0.67) was found.The moderate burn severity class had a relatively low PA value (approx.0.6) and a moderate UA value (0.77), suggesting moderate omission errors.
For the CBI soil metric, global accuracy decreased (κ = 0.63), although the main trend was similar to the total CBI and CBI veg metrics (Table 3).As in the previous cases, unburned and high burn severity classes had the highest PA (0.96 for unburned and 1.00 for high burn severity).Regarding UA values, the burned class had the highest value (0.93) but high burn severity displayed a low UA value (0.54) indicating important commission errors.PA values for the low and moderate classes also were approximately 0.50, pointing to considerable omission errors.Additionally, UA for the moderate burn severity class had a low value as well (0.60), which indicates considerable commission error too.
Finally, regarding the total burned area map, a relatively high κ statistic value was found (0.85).The PA value of the unburned class and UA value of the burned class was relatively high (>0.95).However, the UA value of the unburned and PA value of burned classes were slightly lower (approximately 0.90), pointing to some commission error in the unburned class and some omission error in the burned one.

Discussion
Knowledge of variables that better define fire damage is a decisive point for mapping burn severity accurately and burned area and, thus, developing management policies to preserve vegetation and soil after wildfires.In our analysis of the fire damage effects on Mediterranean forestry systems affected by mega-fire we observed a high relationship between char fraction and burn severity that might mask the contribution of other covariates.This result agrees with previous research [19,28,48,51,62,[100][101][102].Although, to our knowledge there have been no studies relating Sentinel-2 MSI char fraction to burn damage, some studies [29,74,89,[103][104][105] related Sentinel 2-MSI spectral indices to it.In this way, our study confirms the usefulness of Sentinel-2 MSI data to assess post-fire degree of severity.
The contribution of immediately post-fire LST (post_LST) to model burn severity classes increased from high to low burn severity levels, reaching approximately 25%.This fact shows the potential of post_LST as the indicator of burn severity.Different studies [2,62,64] have pointed out the relationship between LST and burn severity although none of them have analyzed the influence of post_LST on each burn severity level independently.Quintano et al. [64] demonstrated that CBI values and LST (immediately after fire) were significantly correlated (R 2 adj = 0.84), and that two burn severity levels could be statistically distinguished from post-fire LST image.In this sense, Fernández-García et al. [2] found a significant relationship between LST images and CBI values in pine fire prone Mediterranean ecosystems.Quintano et al. [62] corroborated this relationship, showing that the inclusion of post-fire LST in a multinomial logistic regression model enabled the increase of final burn severity map accuracy.Similarly, Zheng et al. [66] and Harris et al. [106] showed the increase of accuracy of the burn severity map based on thermally enhanced spectral indices vs. conventional spectral indices.
Our study also showed that the contribution of pre-fire covariates to model burn severity levels and burned area, although lower than the post-fire covariates contribution, was important.Pre_GV and pre_NPV were the covariates that most contributed to the model (approx.20% for low and high burn severity levels), suggesting the influence of pre-fire characteristics of vegetation on burn severity, which agrees with previous research [18,[107][108][109][110][111][112] that showed that, in addition to topography and weather, the most influential factors on burn severity are vegetation composition and structure.Arkle et al. [108] verified correlations between the amount of live vegetation and burn severity.Agee and Skinner [107], Fernandes and Rigolot [109] and García-Llamas et al. [110] demonstrated a significant relationship between pre-fire vegetation structure (high horizontal and vertical fuel continuity) and burn severity.Other pre-fire vegetation parameters such as plant canopy cover, tree density and size, fine fuel accumulations and forest heterogeneity have been also associated to burn severity [18,[111][112][113].
The three continuous outputs from MaxEnt, representing the probability of presence of high, moderate and low burn severity classes, were combined into a burn severity map, which was independently validated using CBI.Its accuracy was relatively high (κ = 0.81).As a reference, [114] in their recapitulation about research works in coniferous forest in North America that estimated burn severity from Normalized Burn Ratio difference (dNBR) and validated the estimation by CBI, found that if three burn severity levels were used, κ statistic varied from 0.62 [115] to 0.37 [116].Working in Mediterranean areas dominated by Pinus pinaster Ait., Quintano et al. [62] reached a κ = 0.66 when using Landsat MESMA fractions and LST from a unitemporal perspective and a κ = 0.45 when using a post-fire Normalized Burn Ratio (NBR) image.Using three burn severities levels Meng et al. [61] obtained a κ = 0.77 at the subcrown scale, and κ = 0.76 at the crown scale from MESMA WorldView-2 fraction images.
Comparing the accuracy of burn severity estimation for vegetation and soil strata (CBI veg and CBI soil metrics), we observed that overall accuracy reached the maximum value for CBI total and the minimum one when the CBI soil metric was used.We found a higher omission error for low and moderate burn severity classes and a higher commission error for high burn severity class when CBI veg was employed instead of CBI total.This pattern was repeated when CBI soil was taken into account, making the differences compared to the CBI total metric more important.Explanation of the overestimation of high burn severity level when CBI veg and particularly CBI soil were considered should be searched by understanding that soil burn severity and vegetation burn severity have different patterns in relation to the vertical structure of the dominant vegetation, in particular, in broomlands and heathlands habitats.From Figure 2 we can observe that in broomlands (left column) although CBI veg values were relatively high for the three burn severity levels, CBI soil values were low, in particular for the low burn severity plot (0.1 vs. 2.25).The same pattern, although less pronounced, was observed in heathlands habitats (central column).In mature broomlands and heathlands the highest aerial phytomass accumulation were not in the lower strata close to the soil, but over one meter above the soil, as occurs in other Mediterranean communities [117], conditioning the differences in burn severity with higher CBI veg than in soil CBI soil .In these mature stages there is a decrease in available light through the shrubland crown, which influences in the mortality of the lowest branches, and may be a determining factor in the vertical distribution of the fuel and in the community structure [118].However, in oak woodlands both CBI metrics (CBI veg and CBI soil ) followed the same trend, because in general this community suffered from lower burn severity than the others because it is associated with a northern aspect and higher humidity, so the severity was lower in strata, vegetation and soil (personal observation).
The combined use of fraction images and distribution models provided the total modeling processes with a physical meaning that facilitated the interpretation of the continuous outputs.These continuous outputs that can be easily interpreted as suitability surfaces for each target class occurrence are one advantage of using MaxEnt [35].Additionally, the use as covariates of fraction images with physical meaning also made understanding of the four modeling processes easier.Previous studies [49,62,119] already pointed out the advantage of fraction images vs. spectral indices for this reason.Thus, the suitability surfaces built by MaxEnt may provide to forest managers additional information that cannot be obtained in just burn severity or burned area maps.Moreover, forest managers could define the thresholds to classify the continuous images into binaries ones that are best adapted to their specific goals and needs [35,120].Although the proposed methodology should be tested in more fire events and in different ecosystems, we believe that it may be a reference to be compared with the results of new studies.In this study, we evaluated this methodology in a large wildfire that represented one of the most significant wildfires in size and damage located in the European Union [97].Our study area was comprised in a Mediterranean ecosystem, however, the methodology could be generalized to other different ecosystems.

Conclusions
Our study analyzed the relative contribution of different pre-and post-fire variables from multispectral data on burned area delimitation and burn severity estimation (low, moderate and high).Both post-fire covariates (char fraction image and LST) added a percentage of contribution of approximately 80% for all burn severity levels (though char fraction had more contribution than LST for high burn severity and LST more contribution for low burn severity level).The rest 20% of contribution percentage was distributed among pre-fire covariates, in particular between GV and NPV.With regards to the burned area target class, the post-fire char fraction image had the highest contribution to the final model, reaching almost 90%, followed by post-fire LST and pre-fire GV and NPV with similar percentages of contribution.The burn severity map reached a κ statistic value equal to 0.81 when total CBI values were used to validate.If CBI veg values were utilized, the accuracy decreased slightly (κ = 0.76), with a greater decrease found for CBI soil (κ = 0.63).Concerning burned area map accuracy, our proposed method achieved a κ statistic equal to 0.85, which indicates a high accuracy.Our results confirmed that accurate burned area/burn severity maps (three burn severity levels) might be obtained using a distribution model trained with pre-and post-fire multispectral data in Mediterranean ecosystems.Additionally, the continuous suitability surfaces provided by MaxEnt for each burn severity level may provide valuable and easy to interpret information that complements the burn severity and burned area maps.Thus, the new proposed multitemporal method may help forest

Figure 2 .
Figure 2. Examples of field plots showing the CBI soil (brown), CBI veg (green) and CBI (orange) values.Upper row: High burn severity level; central row: Moderate burn severity level; lower row: Low burn severity level; left column: Broomlands; central column: Heathlands; left column: Oak woodlands.

24 Figure 3 .
Figure 3. Flowchart of the proposed method to evaluate fire damage.

Figure 3 .
Figure 3. Flowchart of the proposed method to evaluate fire damage.

Figure 4 .
Figure 4. Some spectra from the definitive post-fire spectral libraries.a) Char spectral library, b) green vegetation spectral library; with the most representative woody species; c) non-photosynthetic vegetation (non irrigated lands) and soil spectral library (bare soil, open mine, urban areas, firebreaks and rock).

Figure 4 .
Figure 4. Some spectra from the definitive post-fire spectral libraries.(a) Char spectral library, (b) green vegetation spectral library; with the most representative woody species; (c) non-photosynthetic vegetation (non irrigated lands) and soil spectral library (bare soil, open mine, urban areas, firebreaks and rock).

24 Figure 5 .
Figure 5. From it, we can observe that pre_GV reveals clearly a previous burned area (blackish tone) 391

Figure 6 .
Figure 6.MaxEnt continuous outputs for high, moderate and low burn severity target classes and burn severity map.a) Target class: High burn severity; b) target class: Moderate burn severity; c) target class: Low burn severity; d) total burn severity map.

Figure 6 .
Figure 6.MaxEnt continuous outputs for high, moderate and low burn severity target classes and burn severity map.(a) Target class: High burn severity; (b) target class: Moderate burn severity; (c) target class: Low burn severity; (d) total burn severity map.

Figure 7 .
Figure 7. a) MaxEnt continuous outputs for the burned area target class and b) burned area map.

81 .A
relatively low user's accuracy (UA) value for the low burn severity class (0.64) could be observed, pointing to relatively high commission error.The moderate burn severity class presented modest producer's accuracy (PA) and UA values (0.76), suggesting moderate commission and omission errors.Unburned and high burn severity classes presented the highest PA (approx.0.90) and UA values (approx.0.95).

Figure 7 .
Figure 7. (a) MaxEnt continuous outputs for the burned area target class and (b) burned area map.

Table 1 .
Endmember signatures included in the definitive spectral libraries and summary of unmixing results.

Table 1 .
Endmember signatures included in the definitive spectral libraries and summary of unmixing results.

Table 2 .
Summary of parameters of MaxEnt modeling process.
Note: Values for the average model of the 10 replicas, average values are provided by the MaxEnt software.

Table 3 .
Summary of accuracy parameters.