Jos′eMnuelFern′nez-Guisurg,SusnSu′rez-Seone, Pulo M.Fernnes,Vítor Fern′nez-Grí, AlfonsoFern′nez-Mnso, Crmen Quintno,LeonorClvo

a Area of Ecology,Department of Biodiversity and Environmental Management,Faculty of Biological and Environmental Sciences,University of Le′on,24071,Le′on,Spain

b Centro de Investiga?～ao e de Tecnologias Agroambientais e Biol′ogicas, Universidade de Tr′as-os-Montes e Alto Douro, 5000-801, Vila Real, Portugal

c Department of Organisms and Systems Biology (Ecology Unit) and Research Unit of Biodiversity (IMIB; UO-CSIC-PA), University of Oviedo, Oviedo, Mieres, Spain

d Agrarian Science and Engineering Department, School of Agricultural and Forestry Engineering, University of Le′on, 24400, Ponferrada, Spain

e Electronic Technology Department, School of Industrial Engineering, University of Valladolid, 47011, Valladolid, Spain

f Sustainable Forest Management Research Institute,University of Valladolid-Spanish National Institute for Agriculture and Food Research and Technology(INIA),34004,Palencia, Spain

Keywords:Aboveground biomass Burn severity Landsat LiDAR Pinus pinaster

ABSTRACT Background: The characterization of surface and canopy fuel loadings in fire-prone pine ecosystems is critical for understanding fire behavior and anticipating the most harmful ecological effects of fire. Nevertheless, the joint consideration of both overstory and understory strata in burn severity assessments is often dismissed.The aim of this work was to assess the role of total,overstory and understory pre-fire aboveground biomass(AGB),estimated by means of airborne Light Detection and Ranging (LiDAR) and Landsat data, as drivers of burn severity in a megafire occurred in a pine ecosystem dominated by Pinus pinaster Ait. in the western Mediterranean Basin.Results: Total and overstory AGB were more accurately estimated (R2 equal to 0.72 and 0.68, respectively) from LiDAR and spectral data than understory AGB (R2 ＝0.26). Density and height percentile LiDAR metrics for several strata were found to be important predictors of AGB.Burn severity responded markedly and non-linearly to total (R2 ＝0.60) and overstory (R2 ＝0.53) AGB, whereas the relationship with understory AGB was weaker(R2 ＝0.21). Nevertheless, the overstory plus understory AGB contribution led to the highest ability to predict burn severity(RMSE ＝122.46 in dNBR scale),instead of the joint consideration as total AGB (RMSE ＝158.41).Conclusions: This study novelty evaluated the potential of pre-fire AGB, as a vegetation biophysical property derived from LiDAR, spectral and field plot inventory data, for predicting burn severity, separating the contribution of the fuel loads in the understory and overstory strata in Pinus pinaster stands.The evidenced relationships between burn severity and pre-fire AGB distribution in Pinus pinaster stands would allow the implementation of threshold criteria to support decision making in fuel treatments designed to minimize crown fire hazard.

1. Background

Burn severity, quantified by the loss of organic matter, both aboveground and belowground (Keeley, 2009), is a critical factor shaping ecosystem responses in the early stages of succession(Bastos et al.,2011;Fern′andez-Guisuraga et al., 2021a) that influences long-term ecological effects of wildfires in fire-prone ecosystems(Harris and Taylor,2017).In recent decades, land use changes (Pausas, 2004; Sagra et al., 2019),anthropogenic climate warming (Gonz′alez-De Vega et al., 2016) and a lack of appropriate forest policies aimed at enhancing global change adaptation in the long term (Vil‵a-Cabrera et al., 2018), have promoted the development of fire-prone forest stands with increasing fuel load and continuity(Fern′andez-Guisuraga et al.,2021a),resulting in abrupt shifts in the fire regime of Mediterranean Basin ecosystems (Fernandes et al.,2014; Pausas and Keeley, 2014; Fern′andez-García et al., 2021). In this context, high-severity wildfires are expected to increase (Pausas and Fern′andez-Mu～noz,2012;Chergui et al.,2018),potentially hindering the post-fire recovery of vegetation and decreasing ecosystem resilience(Seidl et al., 2014; Turetsky et al., 2017; Taboada et al., 2018; Fern′andez-Guisuraga et al.,2021a).

Burn severity prediction is emerging as a topic of crucial importance(S′anchez-Pinillos et al., 2021). Among bottom-up fire controls (e.g.vegetation fuel, topography and disturbance history), the nature and spatial variation of fuel structure and loading dictate fire intensity and burn severity at fine scales(Weise and Wright,2014;Keane,2016).Fuel characteristics are the only control that can be handled through adaptative management to reduce high burn severity likelihood in fire-prone ecosystems (Fern′andez-Guisuraga et al., 2021b; S′anchez-Pinillos et al.,2021). Indeed, several studies have evidenced that vegetation structure strongly influence fire spread and burn severity(e.g.Safford et al.,2009;Fernandes et al.,2010;Coppoletta et al.,2016;Harris and Taylor,2017;García-Llamas et al.,2019; Fern′andez-Guisuraga et al.,2021b).

Ecosystems dominated by maritime pine(Pinus pinaster Ait.) are one of the most widely distributed lowland pine woodlands in the Mediterranean Basin(Tapias et al.,2004).Maritime pine stands are intrinsically flammable owing to litter structure and accumulation(Fern′andez-García,2019), resulting in potentially extreme fire behavior, particularly in dense stands with closed canopies focused on maximizing biomass productivity (Fernandes and Rigolot, 2007; G′omez-V′azquez et al., 2013).Additionally,the shrub layer is often tall and dense(Taboada et al.,2018)and is dominated by fine-fuel rich species (Fern′andez-García, 2019),significantly contributing to the surface fuel complex (Castedo-Dorado et al.,2012)and promoting crown fire spread through increased surface fire intensity and vertical continuity(Fernandes and Rigolot,2007).The characterization of both surface and canopy fuels in these ecosystems is then crucial to understand fire behavior and anticipate detrimental ecological effects(Nunes et al.,2019).However,the joint consideration of the overstory and understory layers on burn severity prediction in fire-prone ecosystems is often disregarded (Keane, 2016;S′anchez-Pinillos et al.,2021).

In this context, remote sensing techniques offer nowadays the most feasible alternative for deriving pre-fire fuel load over extensive areas(Garcia et al., 2017) compared to traditional field-based approaches(Saarela et al., 2020; Fern′andez-Guisuraga et al., 2021b), and for modeling the relationships between fuel load and ecosystem wildfire impacts (Viedma et al., 2015). Particularly, active remote sensing data,such as those provided by airborne Light Detection and Ranging(LiDAR)sensors,have been successfully used to estimate three-dimensional forest structure and vegetation biophysical properties, including forest aboveground biomass (AGB), at both plot and stand levels (e.g. Hudak et al.,2012; N?sset et al., 2013; Sheridan et al., 2015; Garcia et al., 2017;Montealegre-Gracia et al., 2017; Csillik et al., 2019; Guerra-Hern′andez and Pascual, 2021). Remarkably, low-density LiDAR is a useful tool to estimate total biomass in forest ecosystems using area-based approaches(Kleinn et al., 2020), including the contribution of overstory and understory layers(N?sset and Gobakken,2008;Domingo et al.,2018).Plot and stand level metrics computed from LiDAR data have also been used to clarify the contribution of pre-fire forest structure to burn severity(e.g.Kane et al.,2015a;Fernandez-Manso et al.,2019).For these applications,the fusion of single-wavelength LiDAR data with passive remote sensing data, such as those provided by optical satellite missions (e.g. Landsat and Sentinel-2), have also yielded satisfactory results to estimate the contribution of pre-fire fuel load to burn severity (Viedma et al., 2020;Fern′andez-Guisuraga et al.,2021b).In this context,satellite optical data can deliver additional information about top of the canopy traits under canopy closure conditions, such as shadowing or moisture content(Avitabile et al., 2012; Vogeler and Cohen, 2016), which can be subsequently correlated with ecosystem three-dimensional structure (Healey et al.,2020).

The traditional approach to model burn severity is based on computing a battery of LiDAR metrics and/or multispectral/hyperspectral products as proxies for pre-fire vegetation structure in the vertical stand profile (Kane et al., 2015a; Viedma et al., 2015, 2020;García-Llamas et al., 2019, 2020; Fern′andez-Guisuraga et al., 2021b),and establishing relationships with burn severity through statistical or machine learning models. These studies indicated that several pre-fire LiDAR or spectral products are robustly correlated with wildfire effects,but they are not intrinsic vegetation biophysical variables than can be generalizable (Fern′andez-Guisuraga et al., 2021a). In turn, pre-fire vegetation biomass distribution in the stand provides a much clear ecological interpretation. However, to the best of our knowledge, no study has directly evaluated the role of pre-fire AGB as a specific fuel biophysical property, derived from LiDAR, spectral and field plot inventory data,on burn severity,and namely separating the contribution of the understory and overstory strata.Here,we evaluated the contribution of pre-fire fuel load (understory, overstory and total biomass components,as intrinsic vegetation biophysical characteristics)to burn severity in a maritime pine ecosystem of the western Mediterranean Basin, by means of pre-fire LiDAR and Landsat data.Specifically,we investigate(i)whether the fusion of low pulse density LiDAR and optical Landsat data allows accurate estimates of pre-fire total AGB and its overstory and understory compartments; and (ii) whether pre-fire AGB (overstory,understory and total)significantly influence burn severity.

2. Methods

The methods comprised five steps: (i) burn severity mapping and validation immediately after fire;(ii)field data acquisition for estimating pre-fire overstory,understory and total AGB at plot level;(iii)LiDAR data acquisition and processing; (iv) pre-fire AGB modeling at plot level through Random Forest regression and AGB prediction at stand level;and(v)data analysis(Fig.1).

2.1. Study site and burn severity mapping

The study site is located in the Sierra del Teleno mountain range(northwest Spain;Fig.2).The site has a heterogeneous relief made up of wide valleys, prominent crests and sedimentary plains, with an altitude ranging between 836 and 1,499 m above the sea level.The climate of the region corresponds to an Atlantic-Mediterranean transition area, with mean annual temperature of 10°C and mean annual precipitation of 640 mm,featuring less than two months of summer drought(Ninyerola et al.,2005).Wildfires are recurrent in this site(free fire interval between 1 and 34 years; García-Llamas et al., 2019) and are mainly associated to spring-summer storms (Santamaría, 2015). In August 2012, a wildfire burned 11,602 ha predominantly occupied by Pinus pinaster stands, as well as isolated stands dominated by Quercus pyrenaica Willd. and Quercus ilex L. The understory (up to 1.5 m in height) and the open shrubland formations were both dominated by Erica australis L.,Halimium lasianthum subsp. alyssoides (Lam.)Greuter and Pterospartum tridentatum

(L.)Willk.Extreme fire weather conditions,including the Haines Index at its maximum value, were recorded during the initiation and growth of the wildfire,resulting in plume-dominated fire behavior(Quintano et al.,2015;García-Llamas et al.,2019).

Landsat 7 Enhanced Thematic Mapper Plus (ETM+) Collection 1 Level-1 pre-fire (September 20th, 2011) and post-fire (September 6th,2012) images (Path 203/Row 31) were acquired from the USGS Earth Explorer data portal(http://earthexplorer.usgs.gov/).The optical bands of both Landsat 7 ETM + scenes were corrected for atmospheric and topographic effects using the ATCOR algorithm (Richter and Schl¨apfer,2018),which was parametrized by means of MODIS water vapor product(MOD05) and meteorological data from both the National Oceanic and Atmospheric Administration(NOAA)and the State Meteorology Agency of Spain (AEMET) (Fern′andez-Guisuraga et al., 2021a), to obtain a surface reflectance product. Burn severity, as a descriptor of total biomass consumption (Keeley, 2009; Morgan et al., 2014), was calculated using the differenced Normalized Burn Ratio(dNBR;Eq.(1)and Eq.(2))index(Key,2006)from surface reflectance data of bands 4(near infrared-NIR-)and 7 (short-wave infrared -SWIR-) of the pre and post-fire Landsat 7 ETM+scenes.

Fig. 1. Methodology flowchart of the present study.

Fig.2. Location of the study site in the western Mediterranean Basin(Iberian Peninsula),burn severity spatial patterns of the 2012 wildfire according to the difference of the Normalized Burn Ratio (dNBR) index and spatial distribution of the field plots of the Third Spanish National Forest Inventory (SNFI-3).

The offset term in Eq. (2) corresponds to the mean dNBR value from pixels in homogeneous and unchanged areas outside the wildfire perimeter(Parks et al.,2014).This index was selected because it is used operationally as the primary spectral index within the Rapid Damage Assessment (RDA) module of the European Forest Fire Information System (EFFIS). In addition, computation of the dNBR index is the most widely used approach and a methodological reference for burn severity initial assessment (Soverel et al., 2010). The remote sensing-based estimation of burn severity was validated with an adapted version(Fern′andez-García et al., 2018) of the Composite Burn Index (CBI; Key and Benson, 2006). This approach is a field standard measurement of burn severity for validating remote sensing products (Holden et al.,2009), and requires the visual assessment of several metrics in four vegetation strata and soil,providing an overall idea of the damage caused by fire(see Fern′andez-García et al.,2018 for details on the modified CBI sampling procedure).The CBI was assessed in the field three months after the wildfire in 54 plots of 30 m×30 m established following a random sampling design. The coefficient of determination of the linear relationship between dNBR and CBI was 0.86.

2.2. Field plot data

The Third Spanish National Forest Inventory(SNFI-3)was used as the field data source to estimate AGB. Within the study site, 40 variableradius plots were established between November 2002 and February 2003. Each circular plot included four nested subplots of 5–25 m radii where trees were inventoried according to their diameter at the breast height (DBH) (S′anchez-Pinillos et al., 2021) (Fig. 3). Height and DBH were measured for each living Pinus pinaster individual with a DBH>7.5 cm. The temporal mismatch between SNFI-3 tree measurements and pre-fire LiDAR data acquisition within the study site (eight years) was addressed using height and DBH growth rates per diametric class provided by the SNFI-3(Domingo et al.,2018).Then,AGB(kg)of each tree in the plot was calculated using Pinus pinaster allometric equations(Ruiz-Peinado et al., 2011) and the sum of the plot values was extrapolated to overstory AGB (t?ha-1) using expansion factors for each tree diametric class. These factors represent the number of trees per hectare that each measured tree represents in relation to the subplot radii(Jim′enez et al.,2017).In addition,the percent cover and mean height of understory shrub species were measured in the 10-m-radius subplot(Fig. 3). Allometric equations for shrub species and shrublands of the Iberian Peninsula (Montero et al., 2020) were used to calculate understory shrub AGB. An annual biomass growing equation (Montero et al.,2020)was applied to match the date of pre-fire LiDAR data acquisition.Finally,total AGB was calculated as the sum of overstory plus understory AGB.

2.3. LiDAR data acquisition and processing

Pre-fire LiDAR data were acquired by the Spanish National Plan for Aerial Orthophotography (PNOA) on October 21st, 2010, using a Leica ALS50 sensor aboard a fixed-wing aircraft.The sensor captured up to four returns per pulse. The mean point cloud density was 0.76 m?2(pulse spacing of 1.15 m). The overall vertical accuracy (RMSEZ) reported by PNOA was lower than 0.2 m. There were no disturbance events within the two years between pre-fire LIDAR data acquisition and the fire date.However, the spatial distribution of fuels was confirmed to be similar during the time lag through photointerpretation of PNOA pre-fire orthophotographs (Fern′andez-Guisuraga et al., 2021b), assuming that LiDAR data were representative of pre-fire vegetation conditions according to Fernandez-Manso et al.(2019).

The LiDAR point cloud was filtered using the multiscale curvature classification algorithm (Evans and Hudak, 2007) implemented in MCC-LIDAR 2.1 software. This algorithm was selected because of its balanced performance when classifying ground points and non-ground points in forest environments (Montealegre et al., 2015, 2016). A digital terrain model (DTM) with a spatial resolution of 2 m was computed from a triangulated irregular network interpolated from the filtered ground returns using LAStools software (rapidlasso GmbH, Germany).Then, LiDAR returns were height-normalized by DTM subtraction, and the normalized point cloud was clipped to the spatial extent of the SNFI-3 field plots.A canopy height model(CHM)with a spatial resolution of 2 m was generated from the interpolation of the highest LiDAR returns with a triangulated irregular network. US Forest Service's FUSION software package version 3.80 (McGaughey, 2018) was used to compute several LiDAR metrics at plot level, ecologically related with overstory and understory AGB.Additionally,metrics were computed at stand level across the study site with a grid size equivalent to that of the burn severity product(30 m).A minimum height threshold of 0.2 m was implemented to remove the influence of non-ground misclassified points and other non-interest ground features such as rocks (Domingo et al., 2018). The metrics comprised: (i) mean height and standard deviation of LiDAR returns;(ii)height percentile values(1st,5th,10th,25th,50th,75th and 95th);(iii)canopy density for several strata(0.2–2,2–4,4–8,8–16 and>16 m) (Kane et al., 2013); (iv) canopy cover (Ma et al., 2017); and (v)vertical complexity index(VCI,Eq.(3);van Ewijk et al.,2011).

where ST is the total number of strata, and piis the proportional abundance of LiDAR returns in stratum i.

Fig. 3. Sampling design within the variable-radius plots of the Third Spanish National Forest Inventory (SNFI-3). Trees are inventoried in the subplots according to their diameter at the breast height (DBH).

2.4. Biomass estimation

Random Forest (RF) regression (Breiman, 2001) was used to model AGB (total, overstory and understory) using as predictors (i) LiDAR metrics at plot level, and (ii) pre-fire surface reflectance of Landsat 7 ETM + optical bands (B1–B5 and B7). RF is an ensemble learning algorithm based on classification and regression trees(CART;Oliveira et al.,2012).The algorithm fits multiple CARTs through bootstrap aggregating techniques, which help to reduce overfitting issues and improve the stability and accuracy of the algorithm(Cutler et al.,2007).RF can also handle spatial autocorrelation in the predictors (García-Llamas et al.,2020). The internal out of bag error rate was used to compute the variance explained (pseudo-R2) by the RF, minimizing the reliance on independent validation datasets(Cutler et al.,2007).RF model parameter mtry was tuned to find the optimal value that minimize the error in the estimation(Liaw and Wiener,2002;Oliveira et al.,2012).Additionally,ntree RF parameter was set to the recommended value of 1000 for obtaining stable model predictions (Probst and Boulesteix, 2018). The relative importance of the predictors in the model was evaluated by means of the percentage increase in mean square error (%IncMSE),which represents the decrease in explained variance if a variable is dropped from the model (Fern′andez-Guisuraga et al., 2020). A forward model selection technique proposed by Kane et al. (2015b) for RF regression was used to select the key LiDAR and Landsat ETM+variables that maximize AGB prediction performance for each stand layer (total,overstory and understory),providing a much clearer interpretation and a more robust model. Partial dependence plots were generated for each selected predictor to disentangle their role on AGB estimation and its ecological significance.Finally,RF regression model objects were used to generate total,overstory and understory AGB prediction maps across the study site by means of the LiDAR metrics at stand level and pre-fire Landsat 7 ETM+reflectance bands.

All analyses were implemented in R (R Core Team, 2020) using the“RandomForest”(Liaw and Wiener,2002),“caret”(Kuhn,2020),“raster”(Hijmans,2021) and“rgdal” (Bivand et al.,2021) packages.

2.5. Pre-fire biomass contribution to burn severity

Pre-fire areas dominated by open shrublands were masked using the CHM(Alonzo et al.,2020),discarding areas with vegetation height lower than 3 m based on field knowledge. Then, pre-fire PNOA orthophotographs were used to sample two hundred polygons stratified within homogeneous patches of Quercus sp.and Pinus pinaster stands.Normalized difference vegetation index (NDVI) values, computed from surface reflectance data of bands 4 (NIR) and 3 (red) of the pre-fire Landsat 7 ETM + scene, were extracted from the polygons of each stand class. A fixed-threshold approach based on standard deviation to average subtraction (e.g. Pu et al., 2008; Martín-Sotoca et al., 2019) were used to separate both classes and mask Quercus sp. stands (NDVI ≥0.74). The accuracy of the non-pine mask was measured through one thousand pixels randomly sampled within the study site and a confusion matrix,using as reference dataset the pre-fire PNOA orthophotographs. The overall accuracy of the mask was 95.20%. Commission and omission errors were equal to 2.63%and 6.77%,respectively.

One thousand pixels were randomly sampled within the non-masked areas of the study site, ensuring a minimum distance of 30 m, for extracting the pre-fire total, overstory and understory AGB, as well as dNBR values. The relationships between burn severity (dependent variable)and AGB at the three levels were evaluated through scatterplots and the coefficient of determination (R2) of the corresponding univariate linear regression models.Models were fitted using a linear or a quadratic term to account for potential non-linear relationships. Finally,the point sample was randomly split into training (70%) and validation (30%)subsets. RF regression was used to generate dNBR predictions based on(i)total AGB,(ii)overstory AGB,(iii)understory AGB,and(iv)overstory and understory AGB as two independent predictors in the same model to account for total AGB differentiated by strata. The root-mean-squared error(RMSE)was computed to measure the prediction performance.

Although the Landsat pre-fire image was used for the estimation of both (i) burn severity through the dNBR thresholding approach(together with the post-fire image) and (ii) pre-fire AGB, a high correlation between the contribution of pre-fire Landsat-derived AGB and dNBR would not necessarily be anticipated because of the following considerations. First, pre-fire AGB and dNBR estimation was based on different temporal data sources.Second,the dNBR index estimates the loss of aboveground biomass as a consequence of fire (Key, 2006), but pre-fire AGB is not directly related to the degree of biomass consumption, i.e. high biomass is not always conductive to high burn severity(Viedma et al.,2020;Fern′andez-Guisuraga et al.,2021b).For example,mature Ponderosa pine and mixed-conifer forests of western North America typically exhibit high pre-fire biomass values but are characterized by low-moderate burn severity regimes where loss of biomass is low because of the fuel distribution in the stand(e.g.Odion et al.,2014;Lesmeister et al., 2019). Third, burn severity was independently measured in the field through the CBI for validating dNBR index burn severity estimates, featuring a high coefficient of determination of the linear relationship.

3. Results

Total and overstory layer AGB were accurately predicted from LiDAR metrics and spectral information using RF regression models, with pseudo-R2of 0.72 and 0.68, respectively. By contrast, the prediction of the understory layer AGB was more limited(R2＝0.26)(Table 1).Based on prediction performance maximization and model parsimony, four to five predictors were selected in total and overstory AGB models, while only one predictor was included in the understory model.Density LiDAR metrics for several strata were the most important AGB predictors in the RF models(IncMSE>25%)(Table 1),featuring a direct relationship with total AGB (Fig. 4) and AGB of the overstory (Fig. 5) and understory compartments(Fig.6).LiDAR height percentile metrics were selected as predictors of total and overstory AGB,the 95th percentile being the most important in both RF models(IncMSE>20%)(Table 1).AGB was directly correlated with upper height distribution percentiles (75th and 95th)(Figs. 4 and 5), whereas the relationship with the 25th percentile of vegetation returns in the total AGB model was inverse (Fig. 4). In the latter model, an even distribution of LiDAR returns over the entire vertical profile, defined by increasing VCI values, was closely associated(IncMSE ＝16%) (Table 1) with stands characterized by a high AGB(Fig. 4). Among Landsat 7 ETM + spectral variables, only surface reflectance of band 5 (SWIR) was selected as a predictor in total and overstory RF models (IncMSE >10%) (Table 1), featuring an inverse relationship with AGB(Figs.4 and 5).

Spatially explicit and superimposed prediction maps of total, overstory and understory AGB of Pinus pinaster stands (Fig. 7) show a wide strip along a northwest-southeast axis over the study site, as well as several southerly areas,dominated by mature Pinus pinaster stands with total AGB values over 80 t?ha-1in pre-fire situation.The contribution of the understory to the total AGB is substantial in these areas. The blank strip in the map, parallel to the former, corresponds to a masked area dominated by shrub vegetation and affected by a wildfire that occurred fifteen years earlier than the 2012 wildfire.

In general, pre-fire AGB was a significant controlling factor of burn severity (p < 0.01) in Pinus pinaster stands. Burn severity, estimated through dNBR index, increased strongly with pre-fire total (R2＝0.60)and overstory(R2＝0.53) AGB(Fig. 8a and b), whereas the correlation with understory AGB was weaker(R2＝0.21)(Fig.8c).The relationships were considerably non-linear.

Prediction of burn severity based on total AGB produced higher error(RMSE ＝158.4 in dNBR scale)(Fig.9a),than predictions considering the AGB of overstory and understory compartments as independent variables in the same model (RMSE ＝122.5) (Fig. 9b). Burn severity predictive capacity of each AGB compartment in separate models was weaker,especially in the case of the understory layer(Fig.9c and d).

Table 1 Pseudo-R2 and variable importance of the LiDAR and spectral predictors included in RF models of total,overstory and understory AGB of pine stands,as well as tuned RF model parameter mtry.

Fig.4. Partial dependence plots of the selected LiDAR and spectral predictors in the total AGB model.Red line represents a LOESS smooth curve.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5. Partial dependence plots of the selected LiDAR and spectral predictors in the overstory AGB model. Red line represents a LOESS smooth curve. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

4. Discussion

LiDAR plot metrics describing three-dimensional vegetation structure, coupled with spectral products related with top-of-canopy biophysical traits,have been evidenced to produce reliable AGB models in a wide variety of woodlands around the globe(e.g.Swatantran et al.,2011;Estornell et al.,2012;Latifi et al.,2012;Greaves et al.,2016;Bell et al.,2018; Heiskanen et al., 2019; L′opez-Serrano et al., 2020). However, to our knowledge, this is the first study that demonstrates the potential of low pulse density LiDAR data, complemented with multispectral information, for estimating accurately AGB through the vertical vegetation profile as a driver of burn severity in Pinus pinaster stands. This knowledge is crucial for supporting adequate environmental policies and fuel treatment strategies aimed at reducing the most adverse environmental effects of wildfires (Corona et al., 2015; Fern′andez-Guisuraga et al.,2021b).

4.1. LiDAR data as a proxy for overstory, understory and total AGB

Total and overstory AGB model accuracy were comparable to those from similar research in conifer ecosystems using LiDAR technology(e.g.Skowronski et al., 2007; García et al., 2010; Sheridan et al., 2015),including low pulse density LiDAR data (N?sset and Gobakken, 2008;Domingo et al., 2018; Tojal et al., 2019), in which the correlation coefficient typically ranges between 0.6 and 0.8 (Skowronski et al., 2007).The variables selected in our models were ecologically representative according to the considered stand compartment,which increases model generalization and predictive ability (Bouvier et al., 2015). Despite the ecological relevance of LiDAR density metrics, since they reflect the distribution of fuel loading and canopy openness per strata(Kane et al.,2010; García-Llamas et al., 2019), preference is usually given to percentiles of the height distribution in AGB estimation (Sheridan et al.,2015). However, density metrics in this study were important in RF models, supplementing upper height distribution percentiles and enabling better characterization of the vertical structure in intermediate parts of the canopy (N?sset and Gobakken, 2005; Chen, 2013) and the stand growth stage (van Ewijk et al., 2011). Indeed, previous research found that the inclusion of LiDAR density metrics in AGB models increased model parsimony due to their high relative importance(Sheridan et al.,2015), and,therefore,model robustness and generality(Bouvier et al., 2015). Upper LiDAR percentiles were systematically included in the models due to their stability and direct relationship with mean canopy height in the plot as a proxy for stand AGB (Kane et al.,2010).In contrast,the 25th percentile LiDAR metric can serve as a proxy for shrub and sapling abundance in the understory compartment(Domingo et al., 2018). In fact, lower 25th percentile return heights describe stands in which a substantial proportion of LiDAR returns are concentrated in the lowest bins (Chen, 2013), which is consistent with the inverse relationship between total AGB and the 25th percentile of the height distribution in the RF model.Despite LiDAR metrics were ranked as the most important variables for explaining total and overstory AGB,Landsat 7 ETM+band 5(SWIR)was also selected in RF models owing to the sensitivity of this spectral region to canopy biophysical traits,such as moisture content and shadowing(Fensholt and Sandholt,2003;Karlson et al., 2015), and, therefore, to AGB, as evidenced in previous studies(Avitabile et al., 2012; Dube and Mutanga, 2015; Deo et al., 2017). In fact,these traits tend to increase when moving from open single-layered to closed multi-layered forest canopies with higher AGB(Avitabile et al.,2012), the latter presenting on average lower SWIR reflectance values.

Fig. 6. Partial dependence plot of the selected LiDAR predictor in the understory AGB model.Red line represents a LOESS smooth curve.(For interpretation of the references to colour in this figure legend,the reader is referred to the Web version of this article.)

Fig.7. Spatially explicit maps of predicted AGB by the RF model in the overstory and understory compartments of Pinus pinaster stands,as well as total predicted AGB in pre-fire situation. Blank regions not mapped correspond to masked areas dominated by open shrublands.

Fig.8. Relationship between burn severity(dNBR)index and total(a),overstory(b)and understory(c)pre-fire AGB.The dotted red line represents the fitted line of the linear model. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig.9. Relationship between observed and predicted burn severity(dNBR)from total(a),overstory plus understory(b),overstory(c)and understory(d)pre-fire AGB using RF regression algorithm. The dotted black line represents the 1:1 line.

The strength of the relationships of LiDAR and spectral data with AGB in the understory compartment was weak compared to total and overstory AGB.Several factors can explain this behavior.First,closed canopy conditions,determined by the dominant trees in the stand,decrease the proportion of laser pulses that penetrate the overstory and reach the lower strata (e.g. Skowronski et al., 2007; Su and Bork, 2007; Hill and Broughton,2009;Martinuzzi et al.,2009;van Ewijk et al.,2011),being the understory compartment undersampled in the LiDAR point cloud to be properly characterized (Ferraz et al., 2016). The current mean point cloud density of 0.76 m?2is clearly insufficient to reduce estimation uncertainty of underneath vegetation strata in Pinus pinaster stands.Second,underestimation of canopy structure variables and AGB has been widely reported when dealing with low stature vegetation, especially with understory shrubs (e.g. Estornell et al., 2011a; Glenn et al., 2011;Mitchell et al.,2011;Li et al.,2017),because LiDAR pulses do not usually intercept the top foliage of understory individuals(Mitchell et al.,2011)and ground returns are frequently misclassified where low growth-form vegetation dominates (Estornell et al., 2011b), among other reasons.Third, the species composition in the shrub understory community is heterogeneous and,therefore,shrub AGB allometric equations are much less reliable than in the case of the overstory where each individual tree is surveyed (Ferraz et al., 2016). It is worth mentioning that the model selection routine in the understory RF model specifically considered the LiDAR density metric in the 0.2–2 m stratum, entirely dominated by shrub vegetation,instead of the 25th percentile of the height distribution as in the total AGB model,so that the main spatial patterns in understory biomass distribution may have been captured(Latifi et al.,2016).Despite these limitations, when overstory and understory AGB were evaluated together, i.e., total AGB, the RF model was able to properly capture the pre-fire AGB variability within the Pinus pinaster stands,attaining similar accuracy to the overstory AGB model. As in the study of Ferraz et al.(2016),this behavior is explainable by overstory-driven AGB variability(52.19 ± 32.49 t?ha-1) and, therefore, total AGB estimates may not be significantly biased by the understory stratum(14.16±5.08 t?ha-1).

It should also be emphasized that the temporal mismatch between SNFI-3 field measurements and pre-fire LiDAR data acquisition constitutes a source of uncertainty in AGB estimations despite being addressed by biomass growing equations(Fern′andez-Guisuraga et al.,2022),since vegetation growth rates in fire-prone ecosystems are influenced by the fire regime(Botequim et al.,2015).However,this issue is not expected to cause a significant bias in the present study due to the absence of fire disturbances in the site for more than ten years (Fern′andez-Guisuraga et al., 2019). In addition, biomass growing equations were developed using populations subject to the same fire regime(Montero et al.,2020).

4.2. Role of pre-fire AGB on burn severity in Pinus pinaster forest

Pre-fire AGB was shown to be a meaningful driver of burn severity in a large and high-intensity crown fire in Pinus pinaster forest, consistent with generic expectations (e.g. Amato et al., 2013; Lecina-Diaz et al.,2014; Fernandes et al., 2015; Viedma et al., 2015). The observed non-linear relationships between burn severity and pre-fire AGB distribution in the stand compartments suggests the definition of fuel treatment thresholds in pre-fire decision-making processes. For instance,overstory AGB density values above 50 t?ha-1were associated to severe wildfire damage(Fig.8b),regardless of higher stand biomass values.The observed dNBR saturation at higher AGB is probably an outcome of weak variation in canopy consumption beyond a certain threshold of canopy density (Stocks et al., 2004). Similar non-linear relationships were also found by García-Llamas et al. (2019) between burn severity and vegetation spectral indices computed from optical satellite imagery.Although vegetation indices(and passive optical data in general)are correlated to a certain extent with AGB and are widely used as surrogates in burn severity assessments(Parks et al.,2018;Zald and Dunn,2018),they are not intrinsic physical quantities (Carlson and Ripley, 1997) and are highly context-dependent(Arkle et al.,2012;Fern′andez-Guisuraga et al.,2021c), hindering results comparison with those from this study(Fern′andez-Guisuraga et al.,2021a).

Several field-based studies conducted in conifer forests have shown strong positive correlations between burn severity and overstory AGB as a measure of pre-fire fuel loading(Cocke et al.,2005;Harris and Taylor,2015), which is consistent with the relationships found in this study.Remarkably, the individual contributions in the same model of pre-fire AGB for the overstory plus understory layers led to the highest predictive capacity of burn severity since they better reflect the distribution of biomass in the stand, instead of their joint consideration in the form of total AGB. In fact, understory fuel characteristics are of remarkable interest for modeling fire behavior,including their role in the transition to crown fire and thus high burn severity (Van Wagner, 1977; Cruz et al.,2004). As noted above in this section, the understory AGB model is unlikely to be accurate in an absolute sense,but the variability captured in the spatial patterns of understory fuel load may be enough to identify the triggering of the crowning process.Despite the main trends observed in this study, further efforts are needed to address burn severity in conifer ecosystems by means of remote sensing data with higher quality to enable a more accurate characterization of the surface fuel complex(Xiao et al., 2019). In this context, the use of LiDAR data with higher pulse density becomes mandatory, in spite of the increased acquisition costs(Bouvier et al.,2015).

5. Conclusions

Characterizing surface and canopy fuel loadings in fire-prone pine ecosystems is a crucial endeavor for anticipating the most harmful ecological effects of fire. This study novelty evaluated the potential of pre-fire AGB, as a vegetation biophysical property derived from LiDAR,spectral and field plot inventory data, for predicting burn severity,separating the contribution of the fuel loads in the understory and overstory strata in Pinus pinaster stands. Although remote sensing data characterized total and overstory AGB better than understory pre-fire AGB, primarily due to the low density of the LiDAR point cloud and the undersampling of the understory compartment in the LiDAR point cloud under closed canopy conditions, combining the individual contributions of pre-fire overstory and understory AGB resulted in the highest capacity to predict burn severity. High pre-fire AGB values were associated to severe wildfire damage, but a weak variation in canopy consumption beyond a certain threshold of canopy density was observed.The evidenced non-linear relationships between burn severity and prefire AGB distribution in Pinus pinaster stands would allow the implementation of threshold criteria to support decision making in fuel treatments designed to minimize crown fire hazard. Suggested improvements of the proposed approach for future research would involve the use of high pulse density LiDAR to minimize uncertainties in the characterization of surface fuel loading and structure.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

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 GESFIRE(AGL2013-48189-C2-1-R) and FIRESEVES (AGL2017-86075-C2-1-R) projects; and by the Regional Government of Castilla and Le′on in the framework of the FIRECYL (LE033U14), SEFIRECYL (LE001P17) and WUIFIRECYL(LE005P20) projects. P.M. Fernandes contributed to this article within the framework of the UIDB/04033/2020 project, funded by the Portuguese Foundation for Science and Technology (FCT). J.M. Fern′andez-Guisuraga is supported by a predoctoral fellowship(FPU16/03070)and a research stay grant (EST19/00310) from the Spanish Ministry of Education.

Authors' contributions

J.M.F.-G., S.S.-S. and L.C. conceived and designed the experiment;J.M.F.-G., P.M.F., V.F.-G., A.F.-M. and C.Q. analyzed the data; J.M.F.-G.wrote the first draft of the paper and S.S.-S., P.M.F., V.F.-G., A.F.-M.,C.Q. and L.C. contributed to the writing; L.C. and S.S.-S. acquired the funding and coordinated the study. The authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors would like to thank the Spanish Ministry for the Ecological Transition and the Demographic challenge for providing the Third Spanish National Forest Inventory(SNFI-3) data.

List of Abbreviations

%IncMSE Percentage increase in mean square error

AGB Aboveground biomass

CARTs Classification and regression trees

CBI Composite Burn Index

CHM Canopy height model

It was an immediate30 hit. Within six months sheet-music sales topped half a million copies. The haulting ballad31 expressed all the fears of a soldier far from home, yearning32 to be in his love s arms. Adaptations of Leip s poem have appeared in more than 40 languages

DBH Diameter at the breast height

dNBR difference of the Normalized Burn Ratio

DTM Digital terrain model

ETM+ Enhanced Thematic Mapper Plus

LiDAR Light Detection and Ranging

NDVI Normalized difference vegetation index

NIR Near infrared

PNOA Spanish National Plan for Aerial Orthophotography

RF Random Forest

RMSE Root-mean-squared error

SNFI-3 Third Spanish National Forest Inventory

SWIR Short-wave infrared

VCI Vertical complexity index