*Article* **Comparison of Physical-Based Models to Measure Forest Resilience to Fire as a Function of Burn Severity**

**José Manuel Fernández-Guisuraga 1,2,\*, Susana Suárez-Seoane 3,4, Carmen Quintano 5,6, Alfonso Fernández-Manso <sup>7</sup> and Leonor Calvo <sup>1</sup>**


**Abstract:** We aimed to compare the potential of physical-based models (radiative transfer and pixel unmixing models) for evaluating the short-term resilience to fire of several shrubland communities as a function of their regenerative strategy and burn severity. The study site was located within the perimeter of a wildfire that occurred in summer 2017 in the northwestern Iberian Peninsula. A preand post-fire time series of Sentinel-2 satellite imagery was acquired to estimate fractional vegetation cover (FVC) from the (i) PROSAIL-D radiative transfer model inversion using the random forest algorithm, and (ii) multiple endmember spectral mixture analysis (MESMA). The FVC retrieval was validated throughout the time series by means of field data stratified by plant community type (i.e., regenerative strategy). The inversion of PROSAIL-D featured the highest overall fit for the entire time series (R<sup>2</sup> > 0.75), followed by MESMA (R2 > 0.64). We estimated the resilience of shrubland communities in terms of FVC recovery using an impact-normalized resilience index and a linear model. High burn severity negatively influenced the short-term resilience of shrublands dominated by facultative seeder species. In contrast, shrublands dominated by resprouters reached pre-fire FVC values regardless of burn severity.

**Keywords:** fractional vegetation cover; MESMA; PROSAIL; recovery; Sentinel-2; wildfire

#### **1. Introduction**

The observed and projected increase in the severity of wildfires in the western Mediterranean Basin may hinder the resilience of plant communities to fire disturbance [1]. Wildfire severity, defined as the fire impact on the ecosystem, and operationally estimated from the amount of above- and belowground plant biomass consumed [2], is one of the key determinants of plant communities' recovery in the first post-fire periods [3]. Although field sampling methods are reliable and accurate for assessing plant community resilience to fire, its wide-scale application is limited in large-scale assessments due to its high time and economic cost [4]. In this sense, synoptic observation of the land surface using remote sensing-based techniques offers an effective way to achieve this objective.

Passive optical sensors of moderate spatial resolution, such as those onboard Landsat or Sentinel-2 satellite missions, offer the potential to accurately detect land cover changes and associated processes over extended time series [5,6]. The Landsat and Sentinel-2

**Citation:** Fernández-Guisuraga, J.M.; Suárez-Seoane, S.; Quintano, C.; Fernández-Manso, A.; Calvo, L. Comparison of Physical-Based Models to Measure Forest Resilience to Fire as a Function of Burn Severity. *Remote Sens.* **2022**, *14*, 5138. https:// doi.org/10.3390/rs14205138

Academic Editor: Nikos Koutsias

Received: 16 September 2022 Accepted: 12 October 2022 Published: 14 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

multispectral imagery archive is potentially useful for assessing spatial patterns of land cover change such as post-fire forest dynamics due to the nominal 30 m/20 m spatial resolution, a consistent scale with that of most landscape-scale land cover changes [7]. In addition, their revisit period provides great coverage for multitemporal studies, especially in areas with frequent cloud cover that can deplete the availability of usable imagery [8].

To date, most studies on ecological resilience to fire have been based on the evaluation of plant community greenness recovery through spectral vegetation indices [9,10], among others. However, this product does not represent physical quantities [11], is affected by the background signal of burned areas [12], and suffers from reflectance saturation at high canopy cover in unburned areas [13]. Conversely, physical-based remote sensing techniques, such as pixel unmixing models and radiative transfer models (RTM), can be used to retrieve vegetation biophysical variables (e.g., fractional vegetation cover (FVC)) as fire resilience indicators [14].

Pixel unmixing models are based on the underlying premise that imagery pixels are constituted by several ground features that contribute to the surface reflectance captured by the optical sensor [15], being the pixel FVC, the physical fraction of vegetation cover in that pixel [16]. In this approach, FVC is directly retrieved from remote sensing data using spectral libraries or image endmembers, without initial field data needs [17]. However, endmember collection can be time-consuming in large, heterogeneous burned landscapes [11]. Among these techniques, multiple endmember spectral mixture analysis (MESMA; [18]) is an advanced method that allows the spectra of several endmembers to characterize each pixel ground component (e.g., the vegetation fraction; [19]). Remarkably, each pixel can be resolved independently with a differing number of endmembers to account for the associated terrain variability [17], for instance different vegetation species within the community. This contrasts with more conventional spectral mixture models, such as linear spectral mixture analysis, in which only one spectrum can be incorporated in each endmember [15].

RTM-based approaches simulate the physical relationships between vegetation biophysical variables, such as the FVC, and the plant community reflectance [20]. These models can be inverted using observed surface reflectance data captured by passive optical sensors for retrieving the FVC, to be implemented as a resilience indicator to fire. RTMs do not need to be parameterized with field data specific to the area of interest, which are only needed for validation of the retrieved biophysical variable [11]. Remarkably, RTM physical relationships are not site-specific, and vegetation recovery can then be monitored in burned landscapes that encompass several communities [21]. RTM inversion is usually performed using machine learning techniques, also known as hybrid inversion [22].

Although physical-based models have been used to monitor post-fire vegetation communities [11,23], to date there are no studies in the literature comparing their effectiveness to assess vegetation resilience to fire in plant communities affected by different severity levels. Therefore, the objective of this work was to compare the potentiality of two physicalbased approaches (pixel unmixing and radiative transfer models) applied to passive optical data for evaluating vegetation resilience as a function of burn severity in shrubland communities with different regenerative traits. FVC was used as a resilience indicator retrieved from a time series of Sentinel-2 imagery, using hybrid RTM inversion and MESMA models.

#### **2. Materials and Methods**

#### *2.1. Study Site Description*

The study site lies within the perimeter of a wildfire that burned 9940 ha of shrub and forest plant communities in the summer of 2017 within Sierra de Cabrera (NW Spain; Figure 1). The site has a rugged topography, with prominent crests and steep slopes, and its altitude ranges from 836–1938 m. The climate is temperate Mediterranean, with a mean annual temperature of 9 ◦C and mean annual precipitation of 850 mm. The wildfire affected, among others, three types of shrub plant communities: shrub communities dominated by facultative seeder species (*Genista hystrix* Lange (gorse) and *Genista florida* L. (broom)) and resprouter species (*Erica australis* L. (heath)). The shrub plant communities of the site

exhibit high spatial variability due to small-scale differences in post-fire recovery patterns and accumulation of burned debris as a consequence of the fire regime variability. Other plant communities affected by the fire include Pyrenean oak forests dominated by *Quercus pyrenaica* Willd, and Scots pine forests dominated by *Pinus sylvestris* L., as well as grasslands in the valley bottoms. The main land use in the site before the wildfire involved extensive livestock farming.

**Figure 1.** Wildfire perimeter and burn severity estimation from the differenced Normalized Burn Ratio (dNBR) thresholds.

Facultative seeder shrub species are able to resprout from belowground organs and to germinate during post-fire conditions; however, in our study site, their resprouting success and vigor is lower than that of obligate resprouters. The latter species rely on resprouting strategy to regenerate after wildfire because they lack a fire-resistant aerial or soil seed bank.

#### *2.2. Sentinel-2 Imagery, Processing and Burn Severity Calculation*

Sentinel-2 multispectral mission comprises two satellites (Sentinel-2A and Sentinel-2B) as part of the Copernicus program. Sentinel-2 provides thirteen bands at different spatial resolutions in the visible, near infrared and short-wave infrared regions: four bands at 10 m, six bands at 20 m and three bands at 60 m [24]. Sentinel-2 bands at a spatial resolution of 10 m were resampled to 20 m using the nearest neighbor interpolation. The bands at 60 m were discarded, as they are used for atmospheric correction and cloud detection processes and are strongly affected by atmospheric effects [20]. Two Sentinel-2 MSI Level 1C scenes were calibrated to surface reflectance (level 2A) with the Atmospheric/Topographic Correction for Satellite Imagery algorithm version 3 (ATCOR-3; [25]) by correcting for topographic and atmospheric effects. Meteorological data from the State Meteorology Agency of Spain (AEMET), the MODIS water vapor product (MOD05), and a digital terrain model (DTM) were used to set appropriate ATCOR-3 input parameters. The scenes were acquired for immediate pre-fire (13 August 2017) and post-fire (2 September 2017) scenarios to estimate burn severity through the differenced Normalized Burn Ratio (dNBR) index [26]. In this study, we selected the dNBR because it was the spectral index most related to field-based burn severity in internal testing compared to relativized indices, as well as in previous research on the site [27]. In addition, the dNBR is the primary spectral index within the European Forest Fire Information System (EFFIS) and a reference approach for initial burn severity assessment [28], which may improve the comparability of the results of our study. The dNBR was validated through the composite burn index (CBI; [29]), measured in 72 field plots of 20 m × 20 m one month after wildfire (initial burn severity assessment). Burn severity was rated in the field between 0 (unburned) and 3 (high severity). Three field burn severity categories were established based on CBI values of each plot: low (CBI < 1.25), moderate (1.25 ≤ CBI ≤ 2.25) and high (CBI > 2.25). These widely accepted CBI thresholds in the literature correspond to those proposed by [30]. These CBI thresholds were used to define three burn severity categories based on a dNBR thresholding approach using a linear regression model: low (dNBR < 384), moderate (384 ≤ dNBR ≤ 659) and high (dNBR > 659) (Figure 1). The R2 of the linear model was equal to 0.84.

#### *2.3. Physical-Based Models*

Pixel unmixing and radiative transfer models were used to retrieve the FVC from three Sentinel-2 MSI Level 2A scenes (processed from level 1C following the methodology described in Section 2.2) acquired during the biomass peak of the study site between 2017 and 2020: (i) 1-week pre-fire (13 August 2017 (equivalent to pre-fire dNBR calculation)), (ii) 2-weeks post-fire (2 September 2017 (equivalent to post-fire dNBR calculation)), and (iii) 3-years post-fire (18 July 2020).

#### 2.3.1. Multiple Endmember Spectral Mixture Analysis (MESMA)

Candidate endmember spectra for MESMA models were extracted from Sentinel-2 scenes (image endmembers) rather than using spectral libraries (reference endmembers) because image endmembers are acquired at the same resolution of the imagery and are influenced by the same atmospheric imagery corrections [19,31]. The first post-fire Sentinel-2 image (2-weeks post-fire) was used to collect endmembers to ensure the presence of enough non-photosynthetic material corresponding to burned vegetation. We used 500 training polygons consisting of uniform ground patches encompassing a single vegetation, soil or charred material type. Polygon size was set to encompass at least four Sentinel-2 pixels. We ensured a separation between polygons of 100 m and a uniform distribution among plant community types [32]. Training areas for soil, green vegetation (*Genista hystrix*, *Genista florida* and *Erica australis*) and non-photosynthetic vegetation (charred debris) were delineated using very high spatial resolution orthophotographs at a spatial resolution of 0.5 m. The Iterative endmember selection (IES) algorithm [33] was used in this study to select optimal endmembers from the candidate set and improve MESMA computational efficiency and accuracy [19]. We clustered the selected endmembers by the IES algorithm into photosynthetic vegetation (PV), non-photosynthetic vegetation (NPV), and soil spectral libraries. Next, Sentinel-2 pre- and post-fire scenes were unmixed into PV, NPV, soil, and shade fraction images. The performance of all candidate MESMA models was evaluated

using the following requirements for each pixel: (i) the range between minimum and maximum fraction image value was constrained between 0 and 1; (ii) shade fraction values lower than 0.8; and (iii) maximum allowable RMSE equal to 0.025 [23,34]. The constraint in fraction images was selected on the basis of the values physically attainable in the field [32], whereas the shade constraint was chosen to maintain a reasonably high fraction value for the other endmembers [34]. For its part, the RMSE constraint is a standard in the literature [18]. Finally, fraction images were shade-normalized, being the GV shade-normalized fraction of the FVC.

IES, MESMA and shade normalization were implemented in VIPER Tools 2.0, developed by the VIPER Lab at UC Santa Barbara [35].

#### 2.3.2. Hybrid Radiative Transfer Model (RTM) Inversion

The PROSAIL-D RTM, resulting from the coupled PROSPECT-D leaf model [36] and the 4SAIL canopy reflectance model [37], was used to simulate shrubland plant canopy reflectance. PROSPECT-D simulates leaves hemispheric reflectance and transmittance from 400–2500 nm as a function of a number of physiological and biochemical variables at the leaf level. For its part, 4SAIL simulates plant canopy reflectance based on PROSPECT-D simulations, as well as a series of variables related to canopy structure, lighting and viewing conditions [38]. The values of PROSPECT-D and 4SAIL input variables (Table 1) were obtained from satellite imagery metadata, in-depth literature review and field knowledge, considering the variability of the biophysical conditions of the shrubland communities in the study site. Then, PROSAIL-D simulations were executed in direct mode to obtain a dataset of simulated shrubland canopy reflectance and its corresponding FVC, calculated from the viewing conditions, leaf area index and leaf angle of each simulation [11]. The dataset was spectrally resampled to the Sentinel-2 band configuration using its spectral bandwidth and response function. We updated the dataset with 10% of bare soil and nonphotosynthetic vegetation spectra with respect to the total simulations [39]. This amount may be representative in relation to the high variability of the vegetation biophysical conditions in highly heterogeneous plant communities [11].


**Table 1.** Value or range of input variables for PROSPECT-D and 4SAIL models.

The Random Forest (RF; [40]) regression algorithm was used to model the relationship between the simulated Sentinel-2 reflectance at the shrubland canopy level and the corresponding FVC for the dataset generated by PROSAIL-D. In this study, the ntree parameter was set to 500 and the mtry parameter to one third of the number of Sentinel-2 bands, which are the default values. The calibrated RF model was then applied to the reflectance

observed in the Sentinel-2 imagery for obtaining spatially explicit FVC predictions for each of the pixels in the pre- and post-fire time series imagery.

PROSAIL-D parameterization and execution in direct mode, as well as FVC retrieval through RF algorithm were performed in ARTMO [41].

#### *2.4. FVC Retrieval Validation*

In September 2017, we established 60 field plots of 20 m × 20 m in burned areas to evaluate the post-fire FVC retrieval performance through MESMA and PROSAIL-D along the post-fire time series. Likewise, 20 control plots of 20 m × 20 m were established in unburned areas to evaluate pre-fire FVC retrieval. This number of field plots has been shown to be representative of the plant communities' variability in the study site [11]. Plots were stratified by reproductive vegetation strategy (resprouters and facultative seeders). The center of each plot was georeferenced using a sub-meter accuracy GPS receiver in post-processing mode. Both control and burned plots were sampled in September 2017, with burned plots also being surveyed in July 2020. FVC was estimated in each plot as the area of the vertical projection occupied by each community stratum (i.e., herbaceous and shrub strata) in the shrubland communities, using a visual estimation method [42]. The coefficient of determination (R2) was calculated to evaluate the retrieval performance of the FVC through MESMA and PROSAIL-D over the time series.

#### *2.5. Data Analyses*

From the spatially explicit FVC prediction maps for the time series with higher overall accuracy (i.e., PROSAIL-D RTM or MESMA pixel unmixing model), we performed a stratified random sampling of 1000 points, using reproductive vegetation strategy and burn severity categories as strata. A minimum distance of 100 m between points was ensured. For each point, the FVC value was extracted for the considered time series. We computed an impact-normalized resilience index (*Rin*) [43] that represents the recovery of the system property, in this case the FVC, with respect to the impact of the disturbance on that property:

$$R\_{\rm in} = (P\_{\rm trx} - P\_{\rm ti}) / (\mathcal{C}\_0 - P\_{\rm ti})$$

where *Ptx* is the value of the system property at the time point when the resilience is evaluated after the disturbance, *Pti* is the value of the property immediately after the disturbance, and *C0* is its control value. A value of the *Rin* index equal to 1 denotes a full recovery of the property at the considered time point. A linear regression model and subsequent Tukey's HSD test were used to evaluate the effect of regenerative vegetation strategy and burn severity (independent variables), as well as their interaction, on vegetation resilience as measured by the *Rin* index (dependent variable). All statistical analyses were performed in R.4.0.5 [44].

#### **3. Results**

The overall accuracy of the RF algorithm trained with PROSAIL-D model simulations for retrieving FVC from Sentinel-2 imagery (R2 = 0.75–0.79) was substantially higher than that achieved from MESMA models (R<sup>2</sup> = 0.64–0.73), both in the immediate pre- and postfire scenarios, as well as three years after the wildfire (Figure 2). The accuracy of the FVC retrieval for both PROSAIL-D and MESMA models was higher in the pre-fire and long-term post-fire scenarios (R<sup>2</sup> > 0.69) than in the immediate post-fire situation (R2 > 0.64). The FVC estimation for the two shrubland communities considered, dominated by facultative seeder species and by resprouter species, showed no significant under- or overestimation effects over the entire range of field-measured FVC, although the estimates for the resprouters were closer tailored to the 1:1 line (Figure 2).

**Figure 2.** Relationship between the FVC measured in the field and that retrieved from Sentinel-2 imagery for the time series using PROSAIL-D (**A**) and MESMA (**B**) models.

From the PROSAIL-D spatially explicit FVC estimates, we evidenced a significant effect of vegetation reproductive strategy (*p*-value < 0.05) and burn severity (*p*-value < 0.001), as well as their interaction (*p*-value < 0.05), on the resilience of shrubland communities in the study site (Figure 3). High burn severity hindered the short-term resilience of shrubland communities dominated by facultative seeder species, with no FVC recovery to a predisturbance state being observed 3 years after wildfire (Figure 3). In contrast, communities dominated by shrub resprouter species reached pre-fire FVC values 3 years after the fire disturbance, regardless of the burn severity scenario (Figure 3). In this sense, resprouter shrubland communities featured a faster recovery rate in high burn severity scenarios than those communities dominated by facultative seeders.

**Figure 3.** Boxplot showing the relationship between the impact-normalized resilience index (*Rin*) and burn severity in shrubland plant communities dominated by facultative seeders species and resprouter species. Significance of *Rin* predictors (regenerative trait: T; burn severity: S; T × S interaction) is represented as \* (*p*-value < 0.05), \*\* (*p*-value < 0.01), and \*\*\* (*p*-value < 0.001). Lowercase letters denote significant differences in *Rin* at the 0.05 level between burn severity levels within each plant community.

#### **4. Discussion**

The assessment of post-fire recovery trajectories through remote sensing-based techniques is essential for (i) understanding current fire regimes in fire-prone ecosystems of the western Mediterranean Basin [3,45], (ii) providing new insights on the resilience of plant communities at several spatial scales [23,46], and (iii) supporting adaptive management strategies aimed at maintaining ecosystem functions and services endangered by changing fire regimes [47–49].

The present study demonstrated the potentiality of physical-based remote sensing approaches (i.e., radiative transfer and pixel unmixing models) for estimating FVC as an indicator of resilience to fire in several pre- and postfire scenarios of shrubland communities with different vegetation responses. Considerably accurate FVC estimates were achieved with the PROSAIL-D RTM and MESMA approaches considering (i) the complexity of biophysical parameters retrieval in shrubland communities because of the high background signal of non-photosynthetic material and bare soil exposed to the remote sensor, and (ii) the complex mixture of shrub species [11,50].

In any case, we found that hybrid inversion of the PROSAIL-D RTM using the RF algorithm outperformed MESMA models when retrieving FVC in heterogeneous shrubland plant communities. Although both approaches have a solid physical basis and feature an adequate capability to generalize the biophysical parameter retrievals [4], the delineation of spectrally pure endmembers with moderate spatial resolution imagery in pixel unmixing models can be challenging [51], especially in heterogeneous burned landscapes with fine-grained arrangement of vegetation and burned legacies. In addition, spatiotemporal changes in the vegetation's biophysical properties and background features when dealing with time series analysis may not be properly captured when delineating image endmembers in extensive burned landscapes [11]. In contrast, the PROSAIL-D model is parametrized with well-known ranges of model input variables for the plant communities under consideration, providing a strong characterization of the physical relationships between the simulated reflectance and site biophysical variability [11]. However, the use of site-specific field information to parametrize the RTM could provide higher accuracy in the retrieval of the biophysical variable of interest [52].

The higher accuracy of FVC retrieval for both PROSAIL-D and MESMA models in the pre-fire and long-term post-fire scenarios, with stronger vegetation responses than in the immediate post-fire situation, could be related to the increased influence of woody debris and bare soil on the surface reflectance at the first post-fire stages with limited canopy cover [50]. This behavior may also be related to FVC estimates closer to the 1:1 line in the case of shrublands dominated by resprouter species, which feature higher vegetation cover [53] and structural complexity [14] than the communities dominated by seeder species in the study site. In addition, a complex mixture of regenerating grass species, seedling recruitment and resprouting responses in early post-fire stages could be encompassed in a decametric Sentinel-2 pixel, increasing retrieval uncertainty [11]. Also, the non-photosynthetic material and bare soil spectra profile from expected pure pixels may not be accurately collected from decametric satellite imagery [22].

Shrubland communities dominated by both facultative seeder and resprouter species featured a high post-fire recovery ability, especially in low and moderate burn severity scenarios, in line with the results of previous field-based research [42,53–55]. However, under high burn severity scenarios, faster recovery rates were evidenced in shrubland communities dominated by resprouter species compared to those dominated by facultative seeders. In general, surviving resprouting structures confer a rapid recovery of plant biomass and a quick recolonization of the physical space occupied by the vegetation prior to the wildfire [54]. Likewise, the aerial and soil seed bank of shrub seeder species suffers considerable damage at high burn severity [56]. This behavior enabled shrubland communities dominated by resprouter species to achieve resilience in the short term after wildfire, in contrast to communities dominated by facultative seeders [57].

Our results are therefore in agreement with those obtained in previous research based exclusively on field data, which demonstrates the potential of physical-based remote sensing techniques, particularly the hybrid inversion of RTMs, to assess the resilience to fire of shrubland communities in the short term. However, both physical-based approaches (i.e., PROSAIL-D RTM and MESMA) featured several FVC retrieval uncertainties in heterogeneous fire-prone shrubland communities. First, PROSAIL-D is a turbid medium RTM, and, consequently, higher accuracies in the model inversion can be attained by using a geometric RTM to simulate canopy reflectance and transmittance in heterogeneous shrubland communities [21], but at the expense of a more complex model parameterization [58], usually with field data not available short-term after fire. Second, RTM approaches exhibit improved performance from retrieval using passive optical data at high spatial resolution, avoiding the land cover aggregation effect of mixed pixels [11]. This shortcoming would be partially solved by MESMA models, which better capture the ground spectra variability of mixed pixels. Nevertheless, non-linear mixing in sparse canopies, such as shrubland communities in the early post-fire periods, violates MESMA assumptions and has an impact on its performance [59]. However, the main limitation of PROSAIL-D RTM and MESMA lies in the impossibility to determine the post-fire recovery trajectories of specific vegetation types or growth forms within the community [31], both in terms of their composition and structure. In this sense, the fusion of remote sensing data from optical sensors processed through physical-based techniques, with that of active sensors such as synthetic aperture radar (SAR) or light detection and ranging (LiDAR), could provide

valuable insights on the recovery trajectories of biophysical properties at the species level or by height strata [14,60,61].

#### **5. Conclusions**

The assessment of how shrubland communities recover from fire disturbance in fire-prone ecosystems is essential to providing new insights about the resilience of plant communities under changing fire regimes in the Mediterranean Basin. This study novelty compared the potential of radiative transfer and pixel unmixing models for evaluating the short-term resilience to fire of several shrubland communities. We found that the hybrid inversion of the PROSAIL-D RTM outperformed MESMA pixel unmixing models to retrieve FVC in heterogeneous shrubland communities. Adaptations to fire allowed shrub communities dominated by resprouter species to achieve resilience in the short-term period after wildfire, which is consistent with previous studies based exclusively on field data, and thus demonstrates the potential of physical-based remote sensing approaches in fire ecology research.

**Author Contributions:** Conceptualization, J.M.F.-G., S.S.-S., A.F.-M., and L.C.; formal analysis, J.M.F.- G., and C.Q.; funding acquisition, S.S.-S. and L.C.; investigation, L.C.; methodology, J.M.F.-G., C.Q., and A.F.-M.; writing—original draft, J.M.F.-G.; writing—review & editing, S.S.-S., C.Q., A.F.-M., and L.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the Spanish Ministry of Economy and Competitiveness and the European Regional Development Fund (ERDF), in the framework of the FIRESEVES (AGL2017-86075-C2-1-R) project; by the Regional Government of Castilla and León in the framework of the WUIFIRECYL (LE005P20) project; and by the British Ecological Society in the framework of the SR22-100154 project, where José Manuel Fernández-Guisuraga is the Principal Investigator. José Manuel Fernández-Guisuraga was supported by a Ramón Areces Foundation postdoctoral fellowship.

**Data Availability Statement:** The datasets analyzed for this study are available from the corresponding author on reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

