Improving silvicultural practices for Mediterranean forests through fire behaviour modelling using LiDAR-derived canopy fuel characteristics

Additional keywords: airborne laser scanning (ALS), fire management, remote sensing, Spanish PNOA project (Plan Nacional de Ortofotografía Aérea de España).

Fire activity in Mediterranean regions is likely to intensify owing to climate change, aggravating its impacts on the environment and society (Salis et al. 2016). Forest structure affects fire hazard by influencing surface, ladder and crown fuels and their involvement in combustion (Cumming 2001; Castro et al. 2003; Fernandes 2009). Thus, ordinary silvicultural management and forest planning should consider fire risk to minimise potential fire damage. Fire risk and wildfire damage can be reduced by removing or modifying fuels in strategic locations (Finney 2006; Ager et al. 2012). In this context, crown fire behaviour considerations are important when planning or evaluating the effectiveness of fuel treatments (Agee and Skinner 2005; González-Olabarria et al. 2012; Jiménez et al. 2016; Rodríguez y Silva et al. 2017). This has prompted the application of fire behaviour modelling software for both research and operational applications, e.g. FARSITE (Finney 1998) or FlamMap (Finney 2006). The assessment of fire behaviour potential, especially the transition from surface to crown fires, may benefit from detailed three-dimensional information about forest structure. However, few studies in the Mediterranean basin have focused on relevant aspects of crown fire modelling to aid fuel management practices and fire suppression strategies.

The high spatial resolution of remote sensing technologies such as airborne LiDAR (light detection and ranging) and aerial imagery can improve canopy fuel attributes estimation and further describe their spatial heterogeneity (Arroyo et al. 2008; Mutlu et al. 2008; Erdody and Moskal 2010; García et al. 2011). Models for crown fire likelihood and behaviour require canopy fuel characteristics as inputs, namely canopy fuel load (CFL), canopy bulk density (CBD), canopy base height (CBH), stand height (SH) and canopy cover (CC). These metrics are spatially explicit inputs for wildfire behaviour models such as FlamMap and FARSITE and should be estimated with reasonably high accuracy (Ager et al. 2010; Botequim et al. 2017). Quantitative description of canopy fuels is thus needed to define silvicultural prescriptions aimed at reducing crown fire susceptibility (Keyes and O’Hara 2002). Crown fires are serious and dangerous events, as fire suppression efforts are more complex than in the case of surface fires, owing to higher spread rate, fireline intensity, smoke production, spotting and turbulence (Cruz and Alexander 2013). The use of LiDAR-derived forest variables may thus strengthen fire behaviour modelling (Kelly et al. 2017).

Several studies using either discrete return or full waveform data have demonstrated the LiDAR potential for estimating canopy fuel metrics (González-Olabarria et al. 2012; Jakubowksi et al. 2013; González-Ferreiro et al. 2014; Hermosilla et al. 2014; Hevia et al. 2016). Typically, discrete systems derive LIDAR metrics related to canopy height and canopy closure, which subsequently are used as independent variables in regression models (Andersen et al. 2005; Hall et al. 2005; Erdody and Moskal 2010; Skowronski et al. 2011; González-Ferreiro et al. 2014). Variables related to canopy height and closure can be interpreted as a 3D representation of canopy structure, thus of canopy fuel distribution. Nevertheless, the higher vertical resolution of full waveform systems as compared with discrete systems may increase the accuracy of CBD estimation (Hermosilla et al. 2014; Crespo-Peremarch et al. 2016, 2018). Countrywide collection of LiDAR data, e.g. by the PNOA project (Plan Nacional de Ortofotografía Aérea; Ministerio de Fomento 2010), which provides low-density (0.5 first returns m⁻²) airborne LiDAR data for the whole of Spain, enables estimation of pixel-based canopy fuel characteristics (González-Ferreiro et al. 2014, 2017) and is a valuable resource for fire management (Alcasena et al. 2016).

Most studies that have used LIDAR data to describe canopy fuels have focused on coniferous forests, which often are characterised by simple stand structures when compared with hardwood or mixed stands. Few studies have addressed the canopy fuel characteristics of mixed forests (Cao et al. 2014; Latifi et al. 2015; Guerra-Hernández et al. 2016a). Further work is needed to address Mediterranean forest structures because only recently, LIDAR studies found advantages of forest-type strata-specific prediction models over non-stratified ones (Guerra-Hernández et al. 2016b).

Landscape-scale fuel treatments decrease fire intensity and fire growth rate, hence resulting in smaller and less severe fires (Fernandes et al. 2015). Studies linking stand structure and fire behaviour characteristics to establish a quantitative basis for forest fuel treatments are scarce outside North America (Martin et al. 2016; Botequim et al. 2017). Spatially explicit fire simulators, e.g. FlamMap (Finney 2006), offer a semi-empirical quantification of fire behaviour and effects and thus a basis to design fuel treatments (Cruz and Alexander 2010; Martin et al. 2016). Yet fire behaviour simulation across time and space requires dynamic weather information and spatially resolved estimates of fuel characteristics, which are not easily obtained (Botequim et al. 2017). This constrains the use of fire simulators in forest management planning. Indirect approaches are desirable such that local stand- or compartment-level fuel treatments can be prescribed without having to turn to advanced simulation models (Botequim et al. 2017).

The present study develops a methodology to estimate canopy fuel characteristics in complex stand structures by using very-low-density airborne LiDAR data. A public forest case study – Tudia y sus Faldas, extending over 748.2 ha and classified into 22 stands – is used to support the research. Specifically, we show how classification tree analysis can be used to detect significant relationships between stand-level features and fire behaviour characteristics. The resulting identification of stand conditions associated with potential crown fire occurrence and difficult fire suppression circumvents the use of fire simulation in a forest management planning context.

The study was conducted in the Tudia y sus Faldas forest located in the region of Extremadura, SW Spain. This public forest (MUP1, Monte de Utilidad Pública número 1) extends over 748.2 ha and is representative of SW Spain stone pine (Pinus pinea L.) forests, including its mix with maritime pine (Pinus pinaster Ait.) or Pyrenean oak (Quercus pyrenaica Willd.) (Fig. 1); small areas of pure Pyrenean oak and maritime pine stands also occur. The forest is managed for pine nut production, provision of food and shelter for wildlife and wildfire prevention. Additional goals are the preservation of mixed stands and regeneration of Pyrenean oak to maintain scenic quality, and the conservation of native flora and fauna.

Fig. 1.  Tudia y sus Faldas public forest: location, forest types map and study plots (black dots).

The climate is Mediterranean with an oceanic influence, characterised by dry summers with rainfall concentrated in autumn and winter, and mean annual temperature and precipitation of 15.2°C and 909.2 mm respectively. The maximum temperature (40.5°C) is reached in August. Average summer precipitation for July and August is just 4.3 mm (1970–2010, AEMET, Agencia Española de Meteorología). Elevation ranges from 300 to 1100 m above sea level, and the average slope is 25.5%. Between 1980 and 2008, 27.1% of the area burned. As per Extremadura’s fire prevention policies (Province Law 207/2005, 30 August), this public forest is assigned a preferential protection area status given its wildfire risk (Level IV, maximum risk).

Field data were obtained from the forest inventory carried out by the Extremadura Forest Service for forest management purposes. In total, 192 circular sample plots of 11-m radius (380.13 m²) were measured between July and August 2010. All trees with diameter at breast height (DBH) >7.5 cm were measured in each plot. Smaller trees were tallied but did not contribute to basal area and volume calculations. Tree height, height of the live crown base (CBH, defined as the lowest insertion point of live branches in a tree), crown width (in two directions at right angles to each other) and age were measured in a subsample of three trees per plot (the northernmost, southernmost and one dominant tree). The heights and crown length of the remaining trees were estimated from locally adjusted allometric equations. Equations were fitted to estimate the height and crown length of all trees in the sample plots (Table 1).

Table 1.  Non-linear height-diameter and crown length (CL) models
h is total tree height (m), d is diameter at breast height (1.3 m above the ground level, mm), Hm is the mean height, G is the stand basal area, Dg is the quadratic mean diameter and CP is the crown projection. RMSE, root mean square error

The individual tree measurements were used to estimate the following canopy fuel variables in each plot: mean height (Hm, m), SH (m), stand basal area (G, m² ha⁻¹), stand volume over bark (V, m³ ha⁻¹), CBH (m), CFL (kg m⁻²) and CBD (kg m⁻³) (Table 2).

Table 2.  Summary of the mean and range of the derived canopy fuel variables

Stand height was computed as Lorey’s mean height (defined as a basal-area-weighted average height), which is the definition of SH that provides the values closest to the aerodynamic canopy height in pine forest stands (Nakai et al. 2010) and is expected to depict the average height of the dominant and codominant trees in a stand (González-Ferreiro et al. 2014). Crown base height was estimated for each tree as the difference between tree height and crown length; then, crown base height was converted to its stand-level equivalent, CBH. Guerra-Hernández et al. (2016b) details the procedures used to obtain the LiDAR and field data. CBD was calculated by dividing CFL by canopy depth (Reinhardt and Crookston 2003). The latter was estimated as the mean crown length of all trees in the plot (Cruz et al. 2003), where crown length is the difference between tree height and height to live crown. For P. pinea and Q. pyrenaica, no published CFL and CBD stand-level equations were found. CFL was calculated as the biomass of needles and twigs <0.6 cm for P. pinea using an equation from Molina et al. (2011) (Table 1). For P. pinaster, we used allometric biomass equations that estimate needles and <0.6-cm twigs for an individual tree (Gómez-Vázquez et al. 2013). For Q. pyrenaica, we used the foliar biomass model of Salazar Iglesias et al. (2010).

The airborne laser scanning (ALS) data were acquired between July and August 2010 for the PNOA project. The laser equipment was a Leica ALS50 sensor operated with a pulse repetition rate of 83 kHz, maximum scan frequency of 32.1 Hz, maximum scan angle of ±50° and average flying height of 2866 m above sea level, which yielded a theoretical density of 0.5 first returns per m². The equipment operates at a wavelength of 1064 nm and can register up to four returns per pulse. Summary statistics of first return density per m² within plots were as follows: average = 1.76, minimum = 1, maximum = 41 and standard deviation = 1.62.

LiDAR metrics are structural descriptive statistics calculated from the normalised laser-derived point cloud. The metrics for the 192 plots were calculated using the FUSION LiDAR Toolkit (McGaughey 2016). LiDAR metrics were computed for each circular plot after normalising the data by subtracting the digital elevation model (DEM). Details of the procedure used to obtain the LiDAR metrics are in González-Ferreiro et al. (2014). Height distribution and canopy closure variables were obtained, and a set of metrics related to the density of returns enclosed in the vertical space defined by five height intervals was calculated, after establishing the minimum height threshold (MHT) and height break threshold (BHT). A combination of threshold heights (MHT and HBT) for each canopy fuel variable was tested to select the best models (Guerra-Hernández et al. 2016b). LiDAR metrics are summarised and described in Table 3.

Table 3.  Potential explanatory LiDAR metrics related to height distribution

The modelling chain consisted of two steps (Fig. 2):

Fig. 2.  Flowchart of the proposed research methodology. CFL, canopy fuel load; CBD, canopy bulk density; CBH, canopy base height; SH, stand height; LiDAR, light detection and ranging; CC, canopy cover; G, stand basal area; Dg, quadratic mean diameter; N, tree density.

1. Establishment of empirical relationships, by forest type, between metrics derived from LiDAR data and measured stand variables using the 192 sample plots. Forest type-specific models were used to estimate canopy fuel complex characteristics, namely CBH, CFL, CBD, SH and CC.

2. Integration into FlamMap of LiDAR-derived canopy fuel characteristics plus topographic and custom fuel model layers to assess potential fire behaviour. Classification tree analysis (CART) was then used to identify significant relationships between stand and fire behaviour characteristics.

The use of low-density discrete LiDAR data for forest mapping is limited, namely smaller trees and understorey are difficult to map reliably (González-Olabarria et al. 2012; Jakubowksi et al. 2013). These studies found that tree height, CBH, CC and general fuel types were accurately predicted from LiDAR data, but specific fuels were difficult to estimate, especially in dense forest (Jakubowksi et al. 2013). Custom fuel models (UCO40 (University of Córdoba)) for Andalusia (Rodríguez y Silva and Molina-Martínez 2012) were assigned to each stand to obtain a surface fuel map for the study area (Table 4). Fuel models allocation was based on expert local knowledge and visual interpretation of georeferenced digital photographs (Table 4, Fig. 3).

Table 4.  Fuel models for Andalusia (UCO40) present in the study area

Fig. 3.  Stand height (a), Canopy fuel characteristics (b,d,e), fuel models (c) and forest types (f) maps for the study area.

Maps were constructed using the best forest type-specific models for LiDAR data and were supported by a regular grid design covering the entire study site. The cell size (380.13 m²) was chosen to match the spatial unit to which the models were fitted (i.e. circular study plots of diameter 22 m). The FUSION LiDAR Toolkit (McGaughey 2016) was then used to obtain the values (per cell) of the explanatory variables and to export these to raster files. The files were then used in a geographic information system (GIS) to produce SH, CBH, CBD and CC maps (Fig. 2). LiDAR-derived CC was estimated using the ratio of the number of first pulses returned from the upper layer of the canopy (using a 2-m threshold height) to all first returns (throughout the canopy to ground profile) (Smith et al. 2009; Morsdorf et al. 2010; González-Olabarria et al. 2012). The DEM was obtained based on (Kraus and Pfeifer 2001) to 1-m horizontal resolution. Slope, aspect and elevation of the study area were mapped from the country’s digital terrain model (DTM) based on standard GIS procedures (Fig. 4).

Fig. 4.  Topographic layers needed for fire behaviour simulation.

We used a landscape file covering 31 356 pixels (size 22 m) and estimated fire behaviour characteristics under constant wind speed and fuel moisture contents using FlamMap5 (Finney 2006). Fuel and stand management planning to reduce fire hazard should address maximum fire potential. FlamMap estimates maximum fire potential for any given point of a landscape, i.e. the characteristics of the forward section of the fire front when wind and slope directions are aligned.

Mapped LiDAR-derived canopy fuel metrics (Fig. 3) and topographic layers (Fig. 4) plus Mediterranean custom fuel models (Rodríguez y Silva and Molina-Martínez 2012) were incorporated into FlamMap to assess potential fire behaviour, namely rate of spread (ROS, m min⁻¹), fireline intensity (FLI, kW m⁻¹), type of fire (surface, passive crown or active crown fire), crown fraction burned (CFB) and fire suppression difficulty across the study area.

Fuel treatments planning should consider critical fire weather conditions, here defined as the May–October 97th percentile (Ager et al. 2010; Salis et al. 2013) of the weather variables distribution from the nearby (~30 km) Fuente de Cantos weather station (2008–15), specifically an air temperature of 37.0°C and a relative humidity of 8.5%. The corresponding fine fuel moisture content was predicted based on the nomograms of Cruz and Alexander (2010), which correspond to a fine (1-h size class) dead fuel moisture content of 4%. For the remaining simulation inputs, we assumed fuel moisture contents of 5% (10-h dead fuels) and 7% (100-h dead fuels); live fuel moisture contents of 75 and 120% respectively for shrubs and canopy foliage and a 6-m wind speed of 40 km h⁻¹.

The modelled fire behaviour variables along with topographic layers over the whole landscape were combined in ArcGIS with standard stand predictors easily measured or predicted through forest inventories (e.g. tree density, N, trees ha⁻¹), and LiDAR models (SH, CC, CFL, CBD and CBH) resulting in a training database with 14 009 unique combinations of variables as observed at the cell level. The corresponding analysis file was the basis to model fire behaviour from site conditions.

Linear models were used to establish forest-specific empirical relationships between field measurements and LiDAR variables. The general expression is as follows (Eqn 1):

where Y are plot-based estimates of canopy fuel or stand variables; X₁, X₂…, X_n are putative explanatory variables expressing the metrics of heights distributions or canopy closure measurements; β₀, β₁,…, β_n are the parameters to be estimated in the fitting process; and ϵ is the additive error term, assumed to be normally, independent and identically distributed with zero mean. Stepwise selection was used to select the best variables for estimating SH, CFL, CBH and CBD, using the leaps package (Lumley and Miller 2009) of R software (R Core Team 2017) Collinearity between regressors was avoided by checking the condition index (CI) and the variance inflation factor (VIF). Regressions with a CI above 30 or VIF above 10 were disregarded (Belsley et al. 2005). Residuals normality was tested with the Shapiro–Wilk test (Shapiro et al. 1968). In contrast to other studies, e.g. González-Olabarria et al. (2012), log-transformation of the variables was not needed to increase model fit or meet the assumptions of the linear regression model.

Comparison of the estimates for the selected models was based on the adjusted coefficient of determination (adj. R²) and the relative root mean square error (rRMSE).

where y_i is the observed value, is the estimated value, is the mean observed value and n is the number of observations.

The relations between input data of fire simulators or stand-level descriptors and site conditions with fire behaviour were fitted through Classification and Regression Tree (CART) analysis to estimate thresholds for dramatic fire behaviour changes in each unique set of cells. CART analysis divides a dataset into increasingly homogeneous subgroups and is ideally suited to model categorical response variables from multiple variables, particularly when non-parametric or unbalanced data and non-linearity are present (De’ath and Fabricius 2000). CART handles autocorrelated data and can disclose complex interactions among predictor variables and quantify their relative importance.

The type of fire predicted by FlamMap (crown fire activity (CFA)) was coded 1, 2 and 3 for surface fire, passive crown fire and active crown fire respectively. In a first attempt (CART I), several variables defining site (e.g. slope, elevation) and stand conditions (e.g. CBH, CBD) and easily measurable or predictable forest inventory variables (e.g. N, G, quadratic mean diameter (Dg)) were tested with forest types (coded 1, 2, 3 and 4 for Q. pyrenaica, P. pinaster, mixed forest and P. pinea respectively). Then, fuel models (FM) were analysed as potential predictors to derive discrimination rules for specific forest management activities for end users ranging from experts with full access to data (CART II) to forest practitioners with little access to data or lacking expertise (CART III). For the latter, the FM set was classified as litter-dominated (‘Litter’), i.e. P1 or HPM3; dominated by short shrubs (‘S_SHRUB’), i.e. M2 or HPM1; or dominated by moderately tall shrubs (‘MH_SHRUB’), i.e. HPM5 or M7. The selection of independent variables was automatic and model overfitting was prevented by basing the number of splits on a 10-fold cross-validation (Fernandes et al. 2008). All statistical analyses were performed using JMP statistical software (SAS Institute 2012).

Fire suppression difficulty maps were based on the FLI (kW m⁻¹) classes of Alexander and Lanoville (1989), namely low (<500 kW m⁻¹), moderate (500–2000 kW m⁻¹), high (2000–4000 kW m⁻¹) and very high to extreme (>4000 kW m⁻¹).

Multiregression linear models (MRM) performed best for most dependent variables, accounting for a high percentage of the total observed variability for all canopy fuel variables (>60%), except CBD in pure and mixed stand of Pinus pinea (Table 5). The regressions explained 61–85, 70–85, 46–96 and 75–95% of the variability in field-estimated SH, CFL, CBD and CBH respectively, with rRMSE ranges of 6.8–20.6, 24.9–37.9, 13.4–48.7 and 5.9–23.0%. The estimates of CFL in P. pinea and P. pinaster stands were more accurate than those obtained for mixed and Quercus pyrenaica stands with a more complex vegetation structure. The highest rRMSE value was for CBD (rRMSE 48.7%) in mixed stands of P. pinea. The scatterplot of this variable (Fig. 5d) evidences some outliers corresponding to significant underprediction of the highest observed values. Although the best rRMSE results were achieved for CBH, the models clearly overestimated this variable in some plots (Fig. 5c).

Table 5.  Summary of stand variables prediction models and plot-level accuracy assessment results for forest type SH (m, Lorey height), CFL (kg m⁻², canopy fuel load), CBH (m, canopy base height), CBD (kg m⁻³, canopy bulk density)
adj. R², adjusted coefficient of determination; RMSE, root mean square error; rRMSE, relative root mean square error

Fig. 5.  Plot-level observed v. predicted values by forest type: (a) stand height (SH); (b) canopy fuel load (CFL); (c) canopy base height (CBH); (d) canopy bulk density (CBD).

Alternatively, CBD was calculated at pixel level (Riaño et al. 2003; Riaño et al. 2004; González-Ferreiro et al. 2014), owing to the poor performance of its direct estimation from LiDAR data for mixed and pure P. pinea forest. We derived CBDa at each cell by dividing CFL by canopy volume. According to the definition of CBD, we calculated this variable as follows (Eqn 5):

where , and are the estimated values of CFL, CBH and SH, which were previously modelled from LiDAR data (Table 5).

Scatterplots of field-measured v. LIDAR-based plot-level values of canopy fuel variables CFL, CBH, SH and CBD (using the direct estimation methodology for P. pinaster and Q. pyrenaica) and CBDa (using the indirect estimation methodology for pure and mixed stands) are shown in Fig. 5.

The MRM showed only minor improvement in predictive ability when including more than one of the height percentile metrics. However, models based on multiple height bin metrics or multiple density metrics did improve the predictive ability of models by forest type. Models utilising height strata metrics always required more predictor variables than models utilising density metrics. The most notable improvement was achieved with models using height strata metrics, where the use of multiple height bin metrics allowed height bin-based models of CBH to account for 75–85% of field-estimated CBH variability for pure stands. Density metrics by strata exhibit strong spatial correlation, and as such, the inclusion of a large number of density bins leads to multicollinearity.

The simulation of fire behaviour characteristics matched local and regional experience on recent wildfires in the area. All fire behaviour characteristics maps were evaluated to identify stand characteristics and spatial pattern metrics of fire-prone areas (Fig. 6). The FlamMap simulation maps showed that fire spread rate (ROS) and FLI can be very high in Tudia y sus Faldas, reaching maximums of 55 m min⁻¹ and 7256 kW m⁻¹ (Figs 6b and 6c). Approximately half of the area was prone to crowning (Fig. 6e), especially untreated areas of mixed stands of P. pinea and P. pinaster. The other half could be burned by surface fire (fuel models HPM1, HPM5 and HPM3). Mixed stands of P. pinea and P. pinaster displayed more than 63% of CFB (Fig. 6f).

Fig. 6.  Fire behaviour characteristics obtained through FlamMap simulation for an extreme scenario: (a) flame length (FL) (m); (b) rate of spread (ROS) (m min⁻¹); (c) fireline intensity (FLI) (kW m⁻¹); (d) heat per unit area (HUA) (kJ m⁻²); (e) crown fire activity (CFA); and (f) crown fraction burned (CFB) (%).

CART I selected three variables to predict CFA, stand metrics CBD and N combined with forest type. CART I explained 66.7% of the existing variation (after eight splits), with areas under the Receiver Operating Characteristics (ROC) curve of 0.953, 0.973 and 0.952 respectively for surface fire, passive crown fire, and active crown fire. The fire type associated with P. pinaster (2) and P. pinea (4) stands with CBD equal to or higher than 0.08 kg m⁻³ was surface fire (96% of the observations) (Fig. 7). However, active crown fire was slightly predominant (54.6%) for the same stand conditions in mixed forests (3) and Q. pyrenaica. Passive crown fire (96.5% of the cases) characterised mixed forest with CBD <0.08 kg m⁻³ and ≥348 trees ha⁻¹.

Fig. 7.  CART I. Classification tree for crown fire activity at the cell level for Tudia y sus Faldas mixed-forest stands as a function of forest type and stand structure (canopy bulk density (CBD), tree density (N)). Numbers at the ends of the terminal nodes are the percentage of observations per type of fire (SF, surface fire; PCF, passive crown fire; ACF, active crown fire).

The CART II analysis for CFA produced a discrimination based on CBD, CBH, forest type and fuel model (R² 0.86, 10 splits). Areas under the ROC curve were 0.991, 0.993 and 0.989 respectively for surface fire, passive crown fire and active crown fire. CBD determined the first CFA partition (Fig. 8). Partitions on the right side of the tree were based on CBD, CBH, FM and forest type, with increasing crowning activity associated with lower CBH (<1.9 m) in mixed-forest stands and Q. pyrenaica with either HPM5 or M7 as fuel models (active crown fire 83%). Fire activity was rated low when CBH ≥ 1.9 m or CBD > 0.08 kg m⁻³, or CBH < 1.9 m combined with fuel models HMP1, HMP3 or M2, and for both pine species (99% surface fire).

Fig. 8.  CART II. Classification tree for crown fire activity at the cell level for Tudia y sus Faldas mixed-forest stands as a function of forest type and fuel model combined with stand structure (canopy bulk density (CBD), canopy base height (CBH)). Numbers at the ends of terminal nodes are the percentage of observations per type of fire (SF, surface fire; PCF, passive crown fire; ACF, active crown fire).

The CART III analysis achieved an R² of 0.85 with 10 splits. Areas under the ROC curve were 0.991, 0.993 and 0.989 respectively for surface fire, passive crown fire and active crown fire. Active crowning occurred in 99% of the observations, all corresponding to Litter and S_SHRUB-dominated stands of maritime and stone pine, coded 2 and 4 respectively. Hence, FM and CBD were the major determinants of the likelihood of crowning (Fig. 9). Fire spread actively in the canopy for 83.5% of MH_SHRUB-dominated understories, mixed forest and Q. pyrenaica stands <1.9 m with CBD ≥ 0.08 kg m⁻³ (Fig. 9).

Fig. 9.  CART III. Classification tree for crown fire activity (CFA) at the cell level for Tudia y sus Faldas mixed-forest stands as a function of forest type and simplified manager fuel models (Litter, S_SHRUB or MH_SHRUB) combined with stand structure (canopy bulk density (CBD), canopy base height (CBH)). Numbers at the ends of the terminal nodes are the percentage of observations per type of fire (SF, surface fire; PCF, passive crown fire; ACF, fire).

A stand-level fire risk map was provided for forest management and risk mitigation purposes (Fig. 10). The area associated with a maximum FLI of 7256 kW m⁻¹ was classified into as extreme in terms of fire suppression difficulty (Alexander et al. 1989). Fire simulation outputs for fire suppression difficulty showed complex patterns that were generally related to the predominant fuel model and topography in the study area (Fig. 10). The shrubland fuel model (M7) was associated with the highest fire suppression difficulty in the south-western and north-eastern areas. The areas that had burned previously (north-east of the study area, occupied by young mixed stands of P. pinea and Q. pyrenaica) and were replanted after the fire carried heavy fuel loads that created favourable conditions for high-intensity fire.

Fig. 10.  Fire suppression difficulty map for Tudia y sus Faldas forest.

This study highlights the value of LiDAR models for characterising the structural properties of canopy fuel layers and helping to address wildfire management concerns. LIDAR provides 3D vegetation models to accurately map fuel properties, which contributes to more reliable fire behaviour estimates. Our approach is potentially valuable for decision support, policymaking and risk mitigation in southern European forests likely to experience large wildfires.

The CFL model fit results in Pinus pinea pure stands are similar to those of González-Olabarria et al. (2012) (adj. R² 0.69, rRMSE 26.7%), but are slightly worse than those reported by González-Ferreiro et al. (2014) (adj. R² 0.82, rRMSE 12.4%). We did not find studies for mixed and Quercus pyrenaica forest, although the values are slightly better than those provided by Cao et al. (2014) for estimating foliage biomass using a discrete LiDAR system in a mixture of coniferous and broad-leaved species (adj. R² 0.44, rRMSE 44. 6%) and for broad-leaved forest (adj R² 0.18, rRMSE 36.4%). The results obtained for Pinus pinaster forest must be treated with caution, given the very small sample size, although they provide some information about the variables that potentially explain the biomass components. The R² values obtained were similar to those obtained by Hevia et al. (2016).

Despite the low-density data of the PNOA project, CBH models were the most accurately modelled of the analysed canopy fuel variables. The adj. R² values (0.75–0.98) confirm that CBH estimation was similar or better than achieved by other studies, where R² ranged from 0.53 to 0.93 (Andersen et al. 2005; Hall et al. 2005; Erdody and Moskal 2010; Zhao et al. 2011; González-Olabarria et al. 2012; González-Ferreiro et al. 2014; Hermosilla et al. 2014). Results confirm that the discrete-pulse LIDAR system systematically overestimates CBH at plot level (Riaño et al. 2004; Hall et al. 2005). The reason is related to vegetation characteristics, such as canopy permeability or CC that could limit LiDAR ability to penetrate the lower parts of the canopy, thus overestimating CBH (Maltamo et al. 2014).

Our CBD modelling results are slightly better than those of González-Ferreiro et al. (2014) probably because first return average density per m² within plots (1.76) was greater than in their study area (0.476). Both studies suggest that the quality of CBD models may be affected by LiDAR data density rather than by plot size. CBD yielded poor results for pure P. pinea stands. This may be also a consequence of varied structural differences between the plots within a large dataset. In the case of mixed forest, the structural complexity and spatial variability of vegetation within the canopy may be another reason for obtaining less accurate models.

The forest type-specific LiDAR models performed better than the model developed from all plots (Cao et al. 2014), which was more obvious on the wall-to-wall mapped area-based predictions. This is consistent with results of Hermosilla et al. (2014) and Latifi et al. (2015).

CBD was highest in P. pinaster stands, particularly in the older stands (Fig. 3). Additionally, CBD and CFL increased in unthinned mixed stands of P. pinea and P. pinaster. Patterns in canopy fuel distribution for pure P. pinea stands were mostly driven by stand structure variables including SH and stand density where predicted CBD pixel values ranged between 0.1 and 0.2 kg m⁻³ across stands.

The lowest CBH and lowest canopy fuel variables are located within mixed stands of young P. pinea and Q. pyrenaica, including along streams and in some riparian buffers. The absence of defined structure within these mixed stands hinders the correct estimation of canopy fuel metrics.

A set of numerical patterns derived from biometric variables was identified based on the classification tree approach to facilitate the identification of thresholds for radical change in fire behaviour and further support hazard-reduction treatments. The CART approach is straightforward and identified thresholds for sudden changes in fire behaviour, providing quantitative evidence for implementing the hierarchy of fuel-reduction treatment principles (Agee and Skinner 2005). Although the general findings of the study are expected to hold in forests with mixed stands elsewhere, model coefficients and variable thresholds are specific for the ranges and combinations of fuel, stand and slope conditions found in the study area. This information may be used to integrate wildfire considerations into long-term forest planning, albeit with the inherent limitations of the fire behaviour models implemented in FlamMap (Cruz and Alexander 2010, 2013).

Fuel-related metrics such as surface fuel characteristics, stand density, vertical continuity and tree size are controllable and can potentially be used as predictors in forest planning systems (Garcia-Gonzalo et al. 2014; Ferreira et al. 2012, 2014).

Forest type in CART I is a predictor variable that combining with lower CBD (<0.08 kg m⁻³) and forest density (<348 trees ha⁻¹) increases crowning likelihood. Previous studies in the Mediterranean Basin considered stand density and height as potential predictors of fire type (Fernandes 2009; Alvarez et al. 2012; Fernández-Alonso et al. 2013). However, only 66.7% of the existing variation was accounted for (Fig. 7). Forest type is consistent in all three CART models. Mixed-forest stands and Quercus pyrenaica remained crown-fire prone, maybe because of their status in the south-west of the study area, especially in the untreated areas of mixed stands of P. pinea and P. pinaster, which is consistent with studies that indicate that older pine stands are less prone to crowning (Fernandes et al. 2015).

The CART analyses differ in their predictive fuel model variables to address the needs of different end users. Two compatible models (CART II and III) were developed to predict the type of fire. We explored the significance of fuel models and found that this predictor increases the accuracy of the models up to a maximum of 86% of variation explained (Fig. 8 and Fig. 9). CART analysis also identified a 0.08-kg·m⁻³ CBD threshold, which is within the 0.05–0.10 kg m⁻³ range for active crown fire development (Agee 1996). Fuel model HPM5 (0.5–1.5-m shrubs) and M7 (1.5–3-m shrubs) combined with low CBH (e.g. 1.5 m) facilitate active crowning.

CBH had the greatest impact on crown fire occurrence as expected from FlamMap adopting the Van Wagner (1977) model for crown fire initiation. A CBH of 1.9 m emerged as a relevant threshold for dramatic changes in the type of fire over the Tudia y sus Faldas forested landscape. Stands with these characteristics should be prioritised for treatment, i.e. reducing surface fuels, increasing CBH and decreasing canopy density through thinning. Our results are consistent with the role played by fuel and stand structure in determining the type of fire (González-Olabarria et al. 2012; Fernández-Alonso et al. 2013; Gómez-Vázquez et al. 2014; Fernandes et al. 2015; Botequim et al. 2017).

The present study achieved a better understanding of how context factors such as forest type and forest management influence fire behaviour and fire suppression strategies. Thus, in light of the need to preserve mixed-forest stands and Q. pyrenaica regeneration, stands with the highest fuel hazard should be considered as priority intervention areas in Tudia y sus Faldas, specially the untreated mixed stands of P. pinea and P. pinaster (south-west of the study area) and previously burned areas characterised by young mixed stands of P. pinea and Q. pyrenaica (north-east of the study area). This is owing to the overcrowded, unthinned young stands creating high fuel continuity and heavy fuel accumulation favouring high-intensity fires.

The aim of this study was to provide multiple forest actors with practical and global discrimination rules to address wildfire management planning. The methods presented have potentially wider application to explore fuel and fire management options in Tudia y sus Faldas. We found a strong relationship between LiDAR-derived metrics and field-based fuel estimates. These results suggest that the freely available low-density LiDAR data provided by the PNOA project are a valuable source of information for estimating canopy fuel variables in Mediterranean forests. The developed models can be used to periodically describe the canopy fuel stratum, whenever the PNOA releases new data, thus reducing financial costs as well as computation, storage and handling efforts. LiDAR data from the PNOA project and the proposed models will be available to all stakeholders and the research community, thus facilitating mapping efforts and enabling a more realistic and accurate description of fire behaviour potential. Future work should enable more detailed fuel models mapping from LiDAR at the scale of the Extremadura region to obtain high-resolution maps for different fuel classification systems based on fuel complex structure from LiDAR-derived metrics in combination with vegetation composition retrieved from LiDAR data and moderate–high-resolution satellite imagery (Sentinel 2A or Worlview-4) and using field data from the Fourth Spanish National Forest Inventory (IFN4).

The second component of this study is the first attempt to develop crown fire behaviour classification rules and to map fire suppression difficulty for a mosaic of Mediterranean forest stands. We attempted to minimise errors with the application of FlamMap by using regional custom fuel models to represent surface fuel beds for the study area. The current research allows crown fire potential to be assessed throughout stand development based solely on forest type or fairly simple inventory data. This allows targeting high-hazard stands for silvicultural interventions and designing appropriate density management alternatives. The CART rules presented allow the development of site-specific prescriptions to treat surface and crown fuels, either emphasising tree survival or fire control, or achieving a compromise between these two objectives. Additionally, maps that combine potential fire behaviour provide the information needed to plan fuel treatments and well-structured firefighting operations, by taking into consideration the existing distinct levels of fire hazard.

Our set of classification rules is useful as (i) it precludes the need to use spatially explicit fire simulation tools, i.e. fire behaviour is assessed solely from fuel and forest structure data; (ii) it helps define effective hazard-reduction silvicultural practices; and (iii) it contributes to developing management guidelines for fuel and stand structure modification in these fire-prone forest stands.

The authors declare no conflicts of interest.