Operational fuel model map for Atlantic landscapes using ALS and Sentinel-2 images

In the new era of large, high-intensity wildfire events, new fire prevention and extinction strategies are emerging. Software that simulates fire behavior can play a leading role. In order for these simulators to provide reliable results, updated fuel model maps are required. Previous studies have shown that remote sensing is a useful tool for obtaining information about vegetation structures and types. However, remote sensing technologies have not been evaluated for operational purposes in Atlantic environments. In this study, we describe a methodology based on remote sensing data (Sentinel-2 images and aerial point clouds) to obtain updated fuel model maps of an Atlantic area. These maps could be used directly in wildfire simulation software. An automated methodology has been developed that allows for the efficient identification and mapping of fuel models in an Atlantic environment. It mainly consists of processing remote sensing data using supervised classifications to obtain a map with the geographical distribution of the species in the study area and maps with the geographical distribution of the structural characteristics of the forest covers. The relationships between the vegetation species and structures in the study area and the Rothermel fuel models were identified. These relationships enabled the generation of the final fuel model map by combining the different previously obtained maps. The resulting map provides essential information about the geographical distribution of fuels; 32.92% of the study area corresponds to models 4 and 7, which are the two models that tend to develop more dangerous behaviors. The accuracy of the final map is evaluated through validation of the maps that are used to obtain it. The user and producer accuracy ranged between 70 and 100%. This paper describes an automated methodology for obtaining updated fuel model maps in Atlantic landscapes using remote sensing data. These maps are crucial in wildfire simulation, which supports the modern wildfire suppression and prevention strategies. Sentinel-2 is a global open access source, and LiDAR is an extensively used technology, meaning that the approach proposed in this study represents a step forward in the efficient transformation of remote sensing data into operational tools for wildfire prevention.


Background
Large wildfires are an increasingly frequent phenomenon around the globe (Tedim et al. 2018;Linley et al. 2022).The 2017 and 2022 wildfire seasons are examples of this trend (Nature, 2017;Rodrigues et al. 2023).In 2017, a total of 622,000 ha burned in Chile, Russia, the USA, Canada, Greenland, and the Mediterranean (Nature, 2017).Similar examples, such as the recent wildfires recorded in California, have arisen in subsequent years (CalFire, 2023).The 2019/20 Australian wildfire season is another example: the wildfires affected, either directly or indirectly, nearly 80% of Australian citizens (Biddle et al. 2020).These large wildfires have severe direct and indirect impacts on human lives, as pointed out by Nature (2017) and Tedim et al. (2018), but they also have significant ecological and socio-economic implications (Li et al. 2022;Taboada et al. 2021;Wang et al. 2021a).According to different scientific records, these events represent a trend that is expected to continue in the future (Dong et al. 2022;Senande-Rivera et al. 2022;Flannigan et al. 2013).
Multiple high-intensity wildfire events have demonstrated that traditional fire extinction and prevention strategies are insufficient for controlling large wildfires.(Moreira et al. 2020;Quílez et al. 2020;Dupuy et al. 2020).For example, the fire prevention infrastructures in Spain are currently still designed for low-intensity fires, e.g., underestimated firebreaks, (Quílez et al. 2020).The need to renovate the traditional management approaches is driving a search for new prevention and extinction strategies that would be effective in controlling these high-intensity wildfires (McWethy et al. 2019;Romero-Vivó et al. 2019;Quílez et al. 2020).In the Mediterranean basin, where large wildfires pose a great threat to some forest ecosystems and human lives and assets, the concept of strategic management areas has been proposed to improve wildfire prevention (Costa Alcubierre et al. 2011;Royé et al. 2019;Oliveira et al. 2018).For instance, in Spain, the recent law, Royal Decree-Law 15/2022, of August 1 st , mandates the identification and design of annually updated strategic management areas throughout the country (RDL 15/2022).
Among the new approaches to wildfire prevention, simulation software is becoming paramount since it allows for the prediction of wildfire behavior in a certain area (Romero-Vivó et al. 2019;Quílez et al. 2020;Iglesias et al. 2022).Wildfire simulation software uses several types of input data to determine the rate of spread, intensity, flame length, and affected area, among other variables.The input data mainly include weather conditions, topographic information, and fuel models (Ellsworth et al. 2022;Rothermel 1983).Fuel models are sets of data that numerically describe the structure and composition of the surface organic matter with the potential to burn, leading to different wildland fire behaviors (Anderson 1982;Benali et al. 2017).These fuel models are the primary input when it comes to the reliability of virtual simulations since fire behavior is highly dependent upon vegetation type and availability as well as on the horizontal and vertical distribution of biomass (Jain et al. 2020;Gale et al. 2021).Accurate and up-to-date fuel model maps of a study area are, thus, essential for effective wildfire prevention and the identification of strategic areas.
According to Simons et al. (2013), most of the existing simulation software relies on Rothermel's (1972) fire spread model and is updated by Albini (1976) using fuel models described by Anderson (1982) and Scott and Burgan (2005).Traditionally, identifying the corresponding fuel model for the specific vegetation types present in an area has been done through field surveys (Anderson 1982;Arroyo et al. 2008).However, using field surveys to obtain updated fuel model maps for large areas can be arduous and potentially not cost-effective (Gale et al. 2021).
In recent decades, new remote sensing technologies have been developed that are providing useful data for fire ecology studies in general (Szpakowski and Jensen 2019) and for the characterization and mapping of fuels in particular (Wallace et al. 2022).Multispectral and optical data can be useful for obtaining geographic information about different types of vegetation in different places and ecosystems (Erinjery et al. 2018;Zeng et al. 2020).For example, Hościło and Lewandowska (2019) used Sentinel-2 images to differentiate between forest types (coniferous and broadleaf ) and between eight different tree species (oak, alder, birch, spruce, pine, beech, larch, and fir).Very high-resolution satellite images (from WorldView-3) were also used successfully to differentiate between species by Fang et al. (2020).This capacity to distinguish between forest types or species is interesting from the forest-firemanagement point of view, since different species and types of forest exhibit different fire behavior when wildfires occur (Wang et al. 2021a, b).However, a limitation of multispectral data is that they do not allow for the detailed characterization of the structure of tree stands.For this purpose, SAR (Synthetic Aperture Radar) data and LiDAR (Light Detection And Ranging) data are being considered.Whereas the potential of SAR data has yet to be fully assessed, the capabilities of LiDAR data to study the structure of forests are well established and understood (Moran et al. 2018;Wiggins et al. 2019;Tello et al. 2018).For example, Jarron et al. (2020) assessed the capability of aerial LiDAR data to estimate sub-canopy attributes and developed a methodology to map these attributes.
In the field of forest fire management, there are also several methodologies based on remote sensing techniques that have been presented in specialized scientific publications.Gale et al. (2021) recently published a review of this topic.They highlighted the existence of two main approaches: using remote sensing data for the characterization of different fuels and using remote sensing data in the automated identification of predefined fuel models.Of the studies analyzed, only 30% use the latter approach (Gale et al. 2021).Some examples are Domingo et al. (2020) and Huesca et al. (2019).They used remote sensing data to identify, at the pixel level, models included in the Prometheus fire model, a model that was developed specifically for Mediterranean ecosystems (Arroyo et al. 2008).To be able to use the obtained fuel maps operationally in simulation software, they designed a procedure to assimilate the Prometheus models to the Rothermel fuel models.Another example is the study of Marino et al. (2016) which developed a fuel model map for the Canary Islands.They used ALS data and Landsat images to create the map, using the Canary Islands Fuel Model (CIFM), a local adaptation of the Rothermel fuel models (Marino et al. 2016).Ferrer Palomino and Rodríguez y Silva (2021), used LiDAR data to characterize, identify, and map fuel models according to the UCO40 fuel models, thus adapting Scott and Burgan's fuel models to a specific Mediterranean environment (Andalucía).An example of a different type of study is the work of Heisig et al. (2022) that, in central European forests, used remote sensing, field work, and modeling techniques to develop a customized fuel map that can be used with the FlamMap software.
An important consideration inferred from the Gale et al. (2021) review is that most of the studies that use remote sensing for fuel mapping focus on Mediterranean forests and woodlands.These kinds of studies, as is pointed out by Cardil et al. (2021), should be performed in a variety of different environments so that wildfire simulation software can be used for operational decision-making in a wider range of areas.For instance, Atlantic climates involve singularities that require specific analysis: vegetation grows faster than in a Mediterranean climate leading to a faster accumulation of large fuel loads which can lead to extreme fire behavior (Arellano et al. 2017).
Galicia is an area representative of an Atlantic climate.It is also one of the Spanish regions that is most prone to wildfires (de Diego et al. 2021).Various efforts have been made to identify and characterize the most common fuel models in Galicia; the most prominent is the fuels photoguide published by Lourizán CIF (photo-guide from now on) (Arellano et al. 2017).Other researchers have characterized some fuel variables through remote sensing; some examples are Alonso-Rego et al. (2021), Alonso-Rego et al. (2020), Fidalgo-González et al. (2019) and Arellano-Pérez et al. (2018).However, these methodologies do not result in a map with the distribution of fuel models; these maps are needed to be used with simulation software to obtain reliable results about the behavior of wildfires in a certain area.In this sense, methodologies that produce updated fuel model maps for Atlantic landscapes are still needed.
In this study, we describe a remote-sensing-data-based methodology for obtaining an operational fuel model map for an Atlantic-vegetation-covered area in Galicia (Northwestern Spain).We used Sentinel-2 images and ALS (Aerial Laser Scanner) data.

Study area
This study was conducted in the region of Galicia, Spain.This region represents only 6% of the total surface area of Spain; however, it encompassed 22.51% of the total burned area in Spain between 2006 and 2015 (Ministerio de Agricultura, Pesca y Alimentación, 2019).This trend has continued in subsequent years (López-Rodriguez et al. 2021;Xunta de Galicia, 2022) with several large wildfires and extensive burned areas being reported (Rodrigues et al. 2023;Xunta de Galicia, 2022).For this reason, improved wildfire prevention plans are needed.For this research, a pilot area was selected.It is shown in Fig. 1.
Galicia has an Atlantic climate (Royé et al. 2016).According to the most recent Spanish Forest Map, approximately 69% of the area studied is covered by forest or shrublands (MITECO, 2011).The main tree species present in the study area are Pinus pinaster Aiton, Pinus radiata D. Don, Eucalyptus globulus Labill., and Quercus robur L., along with riparian forest tree species (Alnus glutinosa L., Fraxinus excelsior L., and Salix spp.L. among others) (MITECO, 2011).The main shrubs in the area are Ulex europaeus L. and Erica spp.L. Obtaining fuel model maps through remote sensing is especially challenging in this region due to the rapid rates of vegetation growth (Arellano et al. 2017), which, in conjunction with an active forestry sector, leads to highly dynamic land cover paradigms.

Satellite images
A time-series of Sentinel-2 images dating from the year 2019 was used in this study.Sentinel-2 is a constellation of two satellites (Sentinel-2A and 2B), each equipped with a Multispectral Instrument (MSI) (European Space Agency, 2015).The MSIs can obtain multispectral images with 13 multispectral bands of different spatial resolutions, ranging from 10 to 60 m.The specifications of the bands are presented in Table 1.Together, the satellites provide a revisit time of five days at the equator and two to three days at mid-latitudes.
Level-2A Sentinel-2 products were used.These products incorporate radiometric, geometric, and atmospheric corrections (European Space Agency, 2023).A total of 12 images were used in this study, one image per month.The images were downloaded from the Copernicus Open Access Hub (https:// scihub.coper nicus.eu/).For each month, the image with the least cloud cover was selected, with the stipulation that no image should have a cloud cover greater than 50%.

Airborne LiDAR data
The LiDAR data were acquired in November of 2019.The point clouds have a mean spatial resolution of 7.2 points/ m 2 .The structure of the data consists of 7786 1 km × 1 km tiles that cover a total area of 7426.35 km 2 of coastal municipalities.The study area includes 378 tiles, corresponding to 6 different municipalities.These 378 tiles cover 366.25 km 2 , which represents 1.2% of the total surface area of Galicia.

Reference data
Open-access high-resolution aerial orthorectified images were used as reference images for the land cover mapping verification process.They were obtained from the Spanish National Plan of Aerial Orthophotography, PNOA by its acronym in Spanish (MTMAU, 2023).In particular, two sets of images were used: (a) images obtained from photogrammetric flights performed between the 30th of May and the 1st of September in 2017 (2017 PNOA); and (b) images obtained from photogrammetric flights performed between the 30th of May and the 1st of September in 2020 (2020 PNOA).The 2017 PNOA images have a spatial resolution of 0.25 m and a georeferencing mean square error of ≤ 0.5 m (MTMAU, 2023).The 2020 PNOA images have a spatial resolution of 0.15 m and a georeferencing mean square error of ≤ 0.20 m (MTMAU, 2023).Both sets of images include 4 bands (Red, Blue, Green, and NIR).

Fuel models and fuel situations
The fuel models that were used were the Rothermel fuel models; these are the models that were created by Rothermel (1972), updated by Albini (1976), and then summarized by Anderson (1982).To become acquainted with the fuel situations that are present in the study area, a fuels photo-guide (Arellano et al. 2017) was used.Table 11 in Appendix 1 presents a detailed description of the Rothermel fuel models, one of the most commonly used fuel model sets in the world (Ascoli et al. 2015;Vacchiano and Ascoli 2015;Simons 2013).
The photo-guide presents a set of real fuel situations in the study area, obtained from destructive field sampling.It offers sample graphic representations of the most representative fuel situations that exist in Galicia.For each fuel situation, it includes quantitative data about the structure of the vegetation.It includes information about dominant species, mean vegetation heights, and fuel loads, as well as representative photographs.Additionally, it provides an estimation of fire behavior in each fuel situation, including information such as the estimated speed of propagation and flame lengths for different wind speeds and terrain slopes under a dead fuel humidity scenario that represents low to moderate risk in Galicia.The data from this fuel photo-guide was obtained through destructive field sampling at the Lourizán CIF (a forestry research station).Table 12 in Appendix 2 summarizes the photo-guide fuel situations present in the study area.

Methodology
The methodology aims to create an automated process that uses remote sensing data to generate a complete, objective fuel map that can be used with wildfire behavior simulation software.The Rothermel fuel models were selected for this purpose since they are the most commonly used in simulation software (Simons et al. 2013).Furthermroe, they are mostly used in the study area and in the whole country since they are the models currently used in the Forest Map of Spain elaboration (MITECO, 2011).However, the Rothermel models are qualitative, and they were developed for a climate that is quite different from that of the study area.Therefore, the relationships between the vegetation types and structures in the study area and the Rothermel models cannot be derived systematically.The fuels photo-guide was used for this purpose.It includes the most common vegetation types in Galicia, classified into fuel situations that are characterized according to qualitative and quantitative parameters.These parameters could potentially be estimated for the study area using remote sensing data.Additionally, the guide provides wildfire behavior data for each fuel situation, facilitating the identification of the relationship between the fuels photo-guide and the Rothermel fuel models.This process allows for the objective and systematic mapping of Rothermel fuel models using remote sensing data from the study area.
The methodology consists of four main steps: (1) identification of the vegetation types in the study area using Sentinel-2 images; (2) characterization of the vegetation height (Canopy Height Model, CHM) and structure (LiDAR metrics) for the study area using ALS data; (3) identification of the correspondences between the fuel situations in the photo-guide, the Rothermel fuel models, and the remote sensing data; (4) mapping of Rothermel fuel models.An overview of the method is shown in Fig. 2.

Land cover mapping
According to the fuels photo-guide, the principal variable that determines the assignment of a fuel situation to a specific plot or forest stand is its dominant species.Sentinel-2 images were used to obtain a land cover map of the study area specifying the principal species and groups of species.The target classes were defined by considering the vegetation in the study area, the aim of the study, and the authors' previous experience in multispectral classification.These classes are detailed in Table 2.The map was created following the methodology described in Alonso et al. ( 2021), which optimizes forest-oriented land cover mapping in Galicia.First, the images were downloaded and pre-processed to mask out the cloudy pixels from each image.This step was performed by applying the cloud mask provided by the Sentinel-2 Level-2A product.
Next, training data were collected through the photointerpretation of 2017 and 2020 PNOA images to perform a supervised classification.Polygons for each target class were manually delineated over the entire study region.Then, a supervised classification was performed for each image in the time series.This time series consisted of 12 images, one per month, spanning the whole year of 2019.This selection of images was performed in such a way as to consider the entire phenological period for the vegetation in the multispectral analysis.According to Alonso et al. (2021), this optimizes the multi-temporal approach in Galicia.The supervised classifications mentioned above were performed using the random forest algorithm (Breiman 2001).This step was performed using the "randomForest" R package (Liaw andWiener 2002, R Core team, 2022) with the default configuration parameters, as indicated in the methodology section of Alonso et al. (2021).The Sentinel-2 bands used in this step were the B02, B03, B04, B05, B06, B07, B08, B8A, B11, and B12 bands.All bands were resampled at 5 m using nearest neighborhood interpolation.This resampling was  performed to adapt the spatial resolution of the multispectral bands of Sentinel-2 to the spatial resolution that can be obtained from LiDAR data.Thus, it is possible to combine the two types of data while at the same time maintaining the details provided by LiDAR data, whose spatial resolution can be higher than the spatial resolution of Sentinel-2.Once single-date classifications were obtained, they were aggregated by plurality voting (Lewiński et al. 2017).The study by Alonso et al. (2021) indicates that this is the aggregation method that best fits the study area, as it reports high accuracy metrics and outperforms other methods in terms of computation time.It consists of selecting as the final class for a pixel the most common class from among all of the classes detected in the single-date classifications.

Structural characterization of vegetation
According to the fuels photo-guide, the assignment of a fuel to a certain area is determined not only by the species present in that area but also by the height and structure of the vegetation.LiDAR data were used to characterize the vegetation height in the study area and to generate an estimation of the forest structure.The process consisted of two main steps: the first one, which includes noise removal, height noise removal, classification, and normalization, is oriented at preparing the data.The second step is the data processing; it includes canopy height modeling, computation of statistics, and rasterization.

LiDAR data pre-processing
The first step of the LiDAR data pre-processing is noise removal, that is, the elimination of isolated points.This was performed using the "lasnoise" tool in the LAStools software (Isenburg 2021).Afterwards, the LiDAR data were classified into two classes: ground points and non-ground points.This was done with the "lasground" tool from the LAStools software (Isenburg 2021).In order to optimize the performance of the tool, the configuration of the parameters of the lasground tool was determined through systematic tests in a representative subset of points from the point cloud of the study area.Several extreme and intermediate values were tested for each of the parameters.In each test, only one parameter was modified to clearly evaluate any gain or loss in efficiency.The optimal classification was considered to be the one in which there were no points classified as ground points above points classified as non-ground points or vice versa.This step was necessary for the subsequent normalization of the point cloud.The normalization transforms the z coordinates.These coordinates, which originally represented the altitude of each point above the Earth's reference surface (according to the official vertical Datum), after normalization, represent the height of each point above the ground.The normalization algorithm applied consisted of a nearest neighborhood (KNN) interpolation with an inverse distance weighting (IDW).It was implemented using the "lidR" package in the R software (Roussel and Auty 2021;Roussel et al. 2020, R Core Team, 2022).The final pre-processing step was the vegetation height segmentation.Ground points and aerial noise were removed so that the data processing could focus solely on the vegetation points.For this step, all points below 0.5 m and above 60 m were removed from the point cloud.

LiDAR data processing
To estimate vegetation height and obtain information about the vertical structure of stands, several statistics were obtained from the normalized point cloud.Each statistic was represented as a raster layer to facilitate its individual analysis, its interpretation using the reference images and its comparison with the other statistics.
First, a canopy height model (CHM) was obtained.The CHM is a raster layer where the digital value of each pixel corresponds with the height of the vegetation at that pixel.The vegetation height is estimated as the maximum value of the z-coordinate, selected from among all the points within the vertical prism defined by the pixel's contour.These values were obtained using the "LidR" package available in the R software (Roussel and Auty 2021;Roussel et al. 2020, R Core Team, 2022).The CHM was produced using a 5-m pixel size.
Additionally, a set of 36 statistics were calculated from the point clouds.Each of the statistics was computed using the normalized values of the z-coordinates of the points; meaning that they ought to provide an approximation of the vertical structure of the point cloud.They were obtained and rasterized using the "LidR" package available in the R software (Roussel and Auty 2021;Roussel et al. 2020, R Core Team, 2022).A raster layer with a spatial resolution of 5 m was obtained for each metric.The complete list of metrics is detailed in Table 3.A detailed explanation of the metric "Cumulative percentage of returns in the X th layer, " can be found in Woods et al. (2008).To obtain these layers, the maximum height of the LiDAR data was divided into 10 intervals.The cumulative percentage of LiDAR returns was calculated for each interval.The final interval was excluded as the cumulative percentage of points in this interval is always equal to 1.Further information on the rest of the metrics can be found in Roussel and Auty (2021).
Finally, horizontal sections were obtained to evaluate the number of points at different levels and the vertical continuity of the vegetation.Specifically, 7 vertical sections were considered.The height thresholds of the sections are summarized in Table 4.It should be noted that the first threshold, 0.5 m, was fixed for the purpose of discarding points at ground level.The last threshold, 12 m, was established since points above that value correspond to high stand canopy.These statistics were transformed into raster layers with a spatial resolution of 5 m.In these layers, the digital value of each pixel corresponds to the number of points within the prism defined by the cell contour and the minimum and maximum z values for the section in question.The rasters were obtained by using the "counter" switch that is available in the "lasgrid" tool in LAStools (Isenburg 2021).

Identification of relationships between photo-guide fuel situations and fuel models
As previously mentioned, the photo-guide presents a set of real fuel situations in the study area.It provides data about fire behavior for each of the situations.Therefore, it was possible to identify relationships between these situations and the fuel models of Rothermel (1972) which were updated by Anderson (1982).In this way, the fuel models that are present in the study area could be established.These relationships were identified as follows: For fuel situations where the main species was a grassland species, models 1, 2, and 3 (Rothermel 1972;Anderson 1982) were considered as possible correspondences.For each grassland situation, the photo-guide provides the following parameters: the main species, the fuel load, and the weighted height.These parameters were compared with the parameters provided for each fuel model.In this way, a relationship was identified between each grassland fuel situation and the most closely corresponding fuel model.
For fuel situations where the major species was a shrub species, models 4, 5, 6, and 7 (Rothermel 1972;Anderson 1982) were considered to be the possible correspondences.Firstly, the parameters provided in the photo-guide for each situation were analyzed: shrub species, fuel loads, and weighted heights.These parameters were compared with the corresponding parameters of the different potential fuel models.Thus, a first link was established between each fuel situation and the most similar fuel model.In order to evaluate the correctness of these links and confirm or discard the correspondence, the fire behavior of each fuel situation was evaluated to see if it was indeed analogous to the fire behavior of the corresponding fuel model.For this purpose, fire behavior charts were created for each potential fuel model with the VisualFuego software (LABIF-UCO, 2019).The charts were built using the same meteorological and topographical variables that were used to develop the photo-guide.They were then compared with the fire behavior charts provided in the photo-guide for the different corresponding fuel situations.In this way, the similarity between each fuel situation and its linked fuel model was confirmed.Where there were no conclusive similarities, the principal of caution was followed: the fuel model representing the worst-case scenario, in terms of fire behavior, was adopted.Where shrub fuel models did not correspond with any shrub fuel situations, in terms of fuel loads and heights or fire behavior, the fuel model was discarded.
Figure 3 shows an example of one of the charts provided by the photo-guide and one of the charts created in Visual-Fuego (LABIF-UCO, 2019) from the fuel model data.The rate of spread for the fuel situation in the photo-guide with a 40% slope is quite similar to the rate of spread for the fuel model with a 40% slope (The 40% slope is a parameter configured by the user).
For fuel situations where the main species is a tree species, the identification of correspondences was performed similarly.In this case, the fuel models considered were 4, 5, 6, 7, 8, 9, 10, and 11.First, the parameters of interest provided in the photo-guide relating to each fuel situation were analyzed: weighted tree height, weighted understory heights, and dead fuel load-dead understory load, slash, and litter loads.These parameters were compared with the  corresponding parameters of the fuel models, and the correspondences were established between each situation and the most similar model.To ensure correspondence between model and situation, fire behavior charts were generated for the fuel models.These charts were compared with those provided in the photo-guide for each fuel situation.Thus, each fuel situation was definitively assigned to the fuel model with the most similar behavior.Where there were no conclusive similarities, the principle of caution was followed as in the case of shrubs.Where a tree fuel model did not correspond with any tree fuel situation in terms of fuel loads and heights or fire behavior, the fuel model was discarded.An extra fuel model, fuel model 0, was added to assign to areas that did not correspond with any particular fuel model, typically built-up areas or bodies of water.
These charts were compared with those provided in the photo-guide for each fuel situation (Table 5).

Identification of correspondences between fuel models and remote sensing classes
Once the fuel models in the study area were identified, the correspondences between them and the land cover classes from remote sensing data could be identified.
The rocky area class from the land cover map was assigned to Rothermel fuel model 1.According to the Rothermel descriptions, fuel model 1 is driven by fine, dry, low grasslands with a light shrub load.Rocky areas correspond to areas mainly covered by rocks, but they usually include patches of scattered grasslands and shrubs.Therefore, as a principle of caution, rocky areas were included in fuel model 1.This particular landscape was not identified in the fuels photo-guide; therefore, this correspondence could not be cross-referenced.The crops and pasture class covers grassland and crop areas in the study area.This class was assigned to Rothermel fuel model 2 since, according to the fuels photoguide and the previously mentioned analysis, it presents a fire behavior typical of this model.
The shrubs class corresponds to shrubs fuel models.Models 4 and 5 can be differentiated by establishing a height threshold according to the Rothermel fuel models height and load parameters (Anderson 1982); the threshold was fixed at 1 m.
Tree-covered vegetation classes included in the land cover map corresponded to the tree fuel models contemplated.Rothermel fuel model 7 is associated to very flammable shrub species with heights of less than 2 m or pine stands such as those found in the state of Florida.However, according to the analysis of the fuels photoguide, diverse tree stands in the study area can present fire behaviors similar to the behavior described for this Rothermel fuel model.Since the main characteristic of these stands is that they present vertical continuity, fuel model 7 was assigned to areas covered by trees that presented a continuous vertical structure.On the other hand, areas covered by trees with a discontinuous vertical structure were associated with Rothermel fuel model 10.A stand was considered discontinuous when the crown base height was at least twice the height of the understory.
Fuel models 8 and 9 were not classified in this study since none of the behaviors of the fuel situations present in the photo-guide corresponded to these models.
Finally, areas potentially covered by significant amounts of slash were assigned to Rothermel fuel model 12. Rothermel fuel models 11 and 13 were not used since they were not found in the study area.
Once the fuel models in the study area were identified, the correspondences between them and the land cover classes from remote sensing data could be identified (Table 6).

Mapping of Rothermel fuel models
Once the correspondences between the remote sensing data and the fuel models were defined, the next step was the generation of a map for each Rothermel fuel model.Different models required different procedures, which are described below.The final step was the integration of all of the models into a single Rothermel fuel model map.a) Models 0, 1, and 2.
These three fuel models were obtained by establishing a direct correspondence between the model and some of the classes from the land cover map.Specifically, all pixels classified as Anthropogenic areas and bodies of water in the land cover map were assigned to fuel model 0. Pixels corresponding to Rocky areas were assigned to fuel model 1 and Crops and pastures were assigned to fuel model 2.
These models were identified through the land cover map and the LiDAR-derived CHM.Models 4 and 5 corresponded with shrubs: model 4 with tall shrubs and model 5 with short shrubs.Therefore, all the pixels depicting shrubs in the land cover map were classified as model 4 or 5. Furthermore, shrub pixels with a CHM value of greater than or equal to 1 m were reclassified as fuel model 4. Shrub pixels with a CHM of less than 1 m were reclassified as fuel model 5.This procedure was performed using the "raster" R package (Hijmans 2021, R Core Team, 2022).c) Models 7, 10 and 12.
Models 7 and 10 corresponded with areas covered by trees, with and without continuous vertical structures, respectively.First, all the pixels depicting trees (Eucalyptus sp., conifers, and broadleaves) on the land cover map were selected.
Once the tree-covered portions of the study area were identified in the land cover map, the vertical continuity of the vegetation in these zones was evaluated.LiDAR metrics were needed for this purpose.Since the vertical structure of vegetation is a complex parameter, it cannot be identified by simply establishing threshold values, as was possible for vegetation height, for example.The automated identification of Rothermel fuel models 7 and 10 was performed through supervised classification.Training and verification areas were defined for this purpose.A set of 90, 25 m × 25 m regular polygons, was designed and distributed across the previously identified tree-covered areas.Half of the polygons represented model 7, and the other half represented model 10.The identification of fuel models in these polygons was carried through visual interpretation of the point clouds.In order to systematize this visual analysis; the 7 horizontal raster layers containing the total number of points per pixel, obtained as described in the "LiDAR data processing" section, were used in this process.Pixels where all 7 raster layers had a high number of points were considered to correspond to stands with vertical continuity and therefore assigned to fuel model 7.
Conversely, pixels where the upper and lower layers presented a much greater number of points than the intermediate layers were considered to correspond to stands with vertical discontinuity, and therefore assigned to fuel model 10.Also, to make this step more objective, vertical sections of the point clouds were generated for the training area polygons.These sections were used to verify that the parameter defined to distinguish discontinuity was met.It was thus confirmed that in the fuel model 10 training polygons, the crown base height was twice or more the height of the understory, whereas in the fuel model 7 training polygons this difference in height was shorter or nonexistent.
The combination of the vertical and horizontal sections of the point clouds allowed for efficient photointerpretation of the samples of point clouds.The 90 regular polygons were randomly divided into training and verification data.A total of 60 polygons were used as training areas, encompassing a total of 1500 pixels.The other 30 polygons were kept aside for posterior validation.
Once the training data were obtained, the supervised classification was performed.This was performed using the machine learning algorithm Random Forest (Breiman 2001); it was applied using the "randomForest" package in R (Liaw andWiener 2002, R Core Team, 2022) with the pre-defined configuration parameters.The predictor variables used in the implementation of the Random Forest predictive model were the 36 metric rasters obtained from the LiDAR point cloud (see the "LiDAR data processing" section).As a result, the Random Forest algorithm creates a model that can be applied to each individual pixel in the study area to predict its corresponding class, or in this case, its fuel model.Afterwards, the model was applied to the selected pixels (pixels representing tree-related classes), resulting in a fuel model map showing fuel models 7 and 10.The importance of the different variables, based on their mean decrease in Gini, was also calculated (Breiman and Cutler 2022) to identify which variables were the most vital in predicting the fuel model.
Fuel model 12 is associated to recently harvested areas.All pixels with tree land cover (either Eucalyptus sp., Conifers or Broadleaves) on the land cover map but that have no data in the LiDAR metrics layers were assigned to this model.As explained previously, the land cover map is obtained by assigning to each pixel the most frequently detected land cover throughout the year in question (in this case, 2019); harvested areas are therefore classified as a tree-related cover.However, since the LiDAR pointcloud was obtained in November of 2019, areas where harvest events had taken place prior to the flight did not present enough returns to calculate the statistics due to the absence of trees or any other vertical structure.

Fuel model map validation
The final map was verified through the validation of the maps that are used to generate it and field observations.The land cover map, which provides the geographical distribution of the types of vegetation in the study area, was cross-verified through a sample of 400 points.They were distributed throughout the study area following a stratified random approach.The ground truth of these points was obtained through the photointerpretation of 2017 and 2020 PNOA images.Sentinel-2 images were used to provide supporting information, to address which PNOA image should ultimately be used to assign the land cover class.For example, if the 2017 PNOA image shows an area covered by conifer trees and the 2020 PNOA image shows that that same area is now covered by shrubs, Sentinel-2 images were used to confirm when the harvesting event took place, allowing us to determine if the land cover class for the year 2019 should be "Shrubs" or "Conifers".A confusion matrix was then constructed and accuracy metrics were calculated.The accuracy metrics calculated were: the Overall Accuracy (OA), the Producer's accuracy (PA), User's Accuracy (UA), and Kappa Index.
The map representing the structure of tree covers (the map with Models 7 and 10 obtained through the classification of LiDAR-derived statistics) was validated through a cross-verification procedure.In this case, the verification dataset comprised of 30 regular polygons.A confusion matrix was then obtained and accuracy metrics were calculated.The accuracy metrics that were calculated are the same as the ones that were used to validate the land cover map.
Finally, all the models identified in the study area were evaluated visually through in situ observations along a track.An experienced operator in wildfire behavior evaluated the compatibility of current situations with the mapped fuel models.

Land cover mapping
The land cover map obtained is presented in Fig. 4. Table 7, presents a summary of the areas occupied by each land cover class.The cross-verification of the map resulted in high accuracy metrics: an Overall Accuracy of 90%, a Kappa Index of 0.88, and a UA and a PA of between 80 and 100% (except in the case of the PA for shrubs, which was lower: 70%).These results are detailed in Table 8, which shows the obtained confusion matrix.It should be noted that cloudy pixels include areas permanently covered by the Sentinel-cloud cover.This includes areas permanently covered by clouds as well as some anthropic areas with very reflective materials such as metallic roofs.

Structure of vegetation
A detailed example of the CHM is presented in Fig. 5.This figure helps visualize the coherence of the CHM.
Lower values are obtained for shrub areas and higher ones for Eucalyptus spp.and conifers.
Examples of the horizontal slices obtained and the number of points in each slice are presented in Figs. 6 and 7.Each example includes a vertical section of the LiDAR point cloud for the selected stand to provide additional context and aid in the comprehension of the examples.Figure 6 corresponds to a tree stand with vertical continuity while Fig. 7 corresponds to a stand with a discontinuous vertical structure.Both figures aid in understanding the relationship between the number of points in the successive slices and the vertical continuity of the stand.Of particular interest are the lower slices (the 2-4-m and 6-8-m slices): Fig. 6 presents more points in these slices than Fig. 7.
Examples of the set of statistics obtained for fuel models 7 and 10 are presented in Figs. 8 and 9, respectively.Together, the two figures illustrate how entropy and the percentage of points in the intermediate layers (the 4th to 6th layers) reflect the vertical structure of the stand.The fuel model 7 stand reports higher values for all of these metrics than the fuel model 10 stand.

Correspondences with Rothermel fuel models
Field work enabled confirmation that rocky areas corresponded to Rothermel fuel model 1.These soils include   rocky areas, but also rocky areas partially covered by grass or very small shrubs (see Fig. 10a).
In relation to Rothermel fuel model 2, field work confirmed that Crops and pastures in the study area do not usually exceed 1 m in height (see Fig. 10b).Furthermore, there is an absence of large grain crops that would lead to fire behavior characteristic of Rothermel fuel model 3.
The shrubs in the study area were assigned to Rothermel fuel models 4 and 5, tall ones and short ones respectively.The most common shrubs in the study area are Ulex spp., Cytisus spp., and Erica spp.Forest areas that are covered predominantly by Cytisus spp.have mainly been assigned to fuel model 4 (see Fig. 10c), while Erica spp.have been assigned to fuel model 5 (see Fig. 10d).Ulex spp.can be found in both models 4 and 5.
In this study, Rothermel fuel model 7 was assigned to areas covered by trees with a continuous vertical structure.These areas are mostly dense young eucalyptus or conifer plantations and stands with mixed species.Broadleaves are included in this category when they are found alongside shrubs or ferns.Figure 10e shows an example of this model.
Conversely, fuel model 10 was assigned to areas covered by trees with a discontinuous vertical structure.These areas mainly correspond to single-species plantations where silviculture actions are commonly or periodically performed.They are usually eucalyptus or conifer plantations or broadleaf stands without a shrub layer.Figure 10f shows an example of this model.
Fuel model 12 was assigned to areas that were identified as harvested areas with forest harvesting residues.An example of this model can be seen in Fig. 10g.

Fuel model mapping
Figure 11 shows the fuel model map obtained for the study area.Table 9 shows a summary of the areas occupied by each fuel.Since a high percentage of the area is covered by agricultural land, fuel model 2 is one of the most predominant models in the study area.Regarding forest-land-related fuel models, the map highlights a greater presence of models 4 and 7 than of 5 and 10.Additionally, fuel models 7 and 10 do not present a mosaic structure but rather, they tend to be concentrated.For example, there is a greater presence of fuel model 10 in the southeast than in the northeast where fuel model 7 prevails among the forested areas.Additionally, when this map is compared with the land cover map it is evident that fuel model 10 is more often associated with areas covered by conifers than areas covered by broadleaves or eucalyptus.
The accuracies of the fuel model map for models 0, 1, and 2 are directly related to the accuracy of the land cover classes used for assigning the models.The accuracy of the fuel model map for fuel models 3, 4, and 12 is also related to the accuracy of the land cover classes used for assigning the models as well as to the precision of the CHM.The accuracy for fuel models 7 and 10 is related to the land cover classes used for assigning the models as well as to the performance of the random forest model built for mapping the models.The performance of the random forest model is presented in Table 10, which shows the confusion matrix performed and the accuracy metrics obtained.The accuracy metrics calculated were all quite close to 100%, and the Kappa Index was 0.9771.
Figure 12 shows a comparison of variables in terms of their importance, and their mean decrease in Gini, for constructing the model.The four most valuable bands in order of importance were entropy followed by the cumulative percentage of returns in the 5 th , 6 th , and 4 th layers.All these variables are tied closely to the vertical distribution of the points in the point cloud, therefore they were able to capture and reflect the vertical continuity or discontinuity of the forest stand.

Discussion
The described methodology uses remote sensing data to obtain an updated operational fuel model map for an area characterized by an Atlantic climate.The map is compatible with most fire behavior simulation software that uses Rothermel's rate of spread equation (Simons et al. 2013).The methodology includes the definition of objective criteria for systematically associating the Rothermel fuel models to the types of vegetation present in the study area.This is not a simple task since the Rothermel Models were described for a different climate region with quite different vegetation types and structures.This adaptation could lead to a standardized way of mapping fuel models in this area.Automating this process makes it more objective and less dependent on the criteria of each individual technician.These types of adaptations are not needed when the fuel model classification is available for the study area (Huesca et al. 2019;Domingo et al. 2020), as was the case in Marino et al. (2016) for the Canary Islands and in Ferrer Palomino and Rodríguez y Silva, (2021) for Andalucía.Various efforts are being made to establish the fire spread equations for Galicia (Arellano et al. 2017).Until they are completed and made ready for use in wildfire simulation software, the 40 fuel models of Scott and Burgan (2005) or the fuel models developed by Portuguese fire behavior researchers (Fernandes et al. 2009;Fernandes and Loureiro 2022) could be an interesting alternative.Both of them include several wooded classes that could help to better describe the varying conditions of the stands in the study area and their wildfire behavior.The usefulness of the Sentinel-2 images together with LiDAR data in the detection and mapping of these fuel models could be explored by taking this study as a starting point.
In the present study, the land cover classification of the study area separates tree stands into three different classes, whereas the fuel model classification divides the three different classes of tree stands into only two fuel  (Arellano et al. 2017).Consequently, this distinction has been made to ensure the maximum potential benefit be derived from remote sensing capabilities in the future in the event that different fire spread equations should be developed for each tree class.
The proposed methodology is similar to the approach taken by Marino et al. (2016).The study of Marino et al. also uses a land cover map derived from satellite imagery (Landsat-8) and ALS data to subsequently build a decision tree to obtain fuel models.However, the methodology they used to create the land cover map was quite different from the method used in this study.Like in this study, they performed a random forest-supervised classification, but rather than using the digital values of spectral bands, like in the present study, they used spectral indices.Neither do they rely on multi-temporal data.The overall accuracy of the land cover of the present study is Fig. 11 Fuel model map of the study area higher (9 percentage points higher) than the OA obtained by Marino et al. (2016).They remarked that the reliability of these kinds of methodologies is highly dependent upon the accuracy of the land cover maps (Marino et al. 2016).However, since the target land cover map classes in this study greatly differ from those used by Marino et al. (2016), a comparison of the accuracies obtained in the two studies might not be fully appropriate.
The most remarkable difference in relation to previous studies is the procedure used for discriminating areas with vertical continuity from those with vertical discontinuity.Marino et al. (2016) and studies relying on Prometheus models (Arroyo et al. 2008;Huesca et al. 2019;Domingo et al. 2020) fix thresholds for height and percentages of points in the point clouds.However, the specific threshold values to be designated are dependent on the LiDAR density, occlusions, and the type of species present (García-Cimarras et al. 2022).As a result, the replicability of the method might be compromised when it is applied in different contexts.In this case, stands are differentiated through a random forest classification, whose results do not depend on local parameters.Furthermore, studies that rely on the Prometheus system depend greatly on the study of the canopy cover (Arroyo et al. 2008).This type of analysis is not feasible when using a 5 × 5-m pixel.
Field work is still commonly used in operational and scientific work to identify fuel models visually (Arroyo et al. 2008).However certain conditions (e.g., weather, stand location, and conditions) prevent operators from performing in situ observations.When forest fuels of vast areas are to be mapped, the execution of field work may compromise affordability and time efficiency.For example, the methods presented by Ferrer Palomino and Rodríguez y Silva (2021) or Heisig et al. (2022) might be difficult to apply on a broader scale since they are designed for pilot areas or for local conditions and therefore the models developed cannot easily be extrapolated to other conditions (Huesca et al. 2019;Ferrer Palomino and Rodríguez y Silva, 2021;Heisig et al. 2022).
In this study, the identification of fuel models in the study area is carried out through the photo-guide that includes the fuel situations in the area.The proposed methodology produces a fuel model map for large areas through remote sensing, minimizing field work.The multispectral data allow for the identification of species and the LiDAR data provide points that penetrate the canopy, bouncing off the vegetation or tree trunks that are under the canopy.Thus, these data allow for the analysis of the vertical structure of the understory and of forest masses.Since field work is not needed, the method ought to be easily applied in the rest of Galicia or in areas with similar conditions.
Another significant advantage of our proposed methodology is that it produces a map with a spatial resolution of 5 m.This is thanks to the LiDAR data, which allows for   Although this methodology is presented as a methodology for obtaining updated operational fuel model maps, its dependence upon LiDAR data could mean that its updatability is compromised.Therefore, between one LiDAR data update and the next, it may be necessary to update the fuel model map using open-access satellite imagery as was done, for example, by DeCastro et al. (2022).A methodology that may prove valuable for this task is the methodology developed by Alonso et al. (2021) for detecting forest disturbances.However, this methodology would be insufficient for detecting changes between models that differ solely in the vertical structure of the vegetation (i.e., it can't differentiate between model 7 and model 10), as this differentiation has thus far only been accomplished with LiDAR data.Having updated fuel model maps is of the utmost importance, and it would be advisable to increase the frequency at which open-access ALS data is updated in Spain (currently it is updated every 7 to 11 years (MTMAU, 2023).Additionally, it is worth exploring in detail for how long model 12 can indeed be assigned to describe an area since harvesting slash is never permanent.
Bearing in mind the high-risk fire regime in this region, the map produced allows for an in-depth study of the fuel model distributions and patterns in the study area.This map may also be useful for stakeholders in the area to enhance fire prevention using new strategies based on the use of fire behavior simulation software (Simons 2013;Romero-Vivó et al. 2019;Moreira et al. 2020;Quílez et al. 2020;Iglesias et al. 2022).Additionally, and in terms of fire prevention, this would help these stakeholders fulfill the new legal requirements set forth by the Spanish government (RDL 15/2022).

Conclusion
This study allowed for the creation of an updated fuel model map for an Atlantic landscape using remote sensing.The fuel map obtained is adapted to the Rothermel fuel models and therefore it can be used operationally with diverse fire simulation software.Additionally, it was based specifically on the fuel types present in the study area and could serve as the basis for future fuel models designed explicitly for this environment.These would aid in the future replicability of the methodology and mark a step towards mapping customized fuel models.Future studies might explore the development of these fuel model maps with other types of fuel classifications such as the Scot and Burgan fuel models or the fuel models developed by Portuguese fire behavior researchers.This might enable the differentiation between stands with the same vegetation class.
The availability of these maps will allow stakeholders in this region to make a shift towards better suppression and fuel management actions and to address the new phenomena of ever-worsening fire scenarios.This is thanks to better, more up-to-date fuel maps for use with fire behavior simulation software.Additionally, these maps would also allow stakeholders to fulfill the legal requirements aimed at fire prevention that are set forth by the Spanish government.

Fig. 2
Fig. 2 Overview of the method

Fig. 3 Table 5
Fig. 3 Example of the identification of the relationship between a fuel situation and a fuel model.The red boxes denote the rate of spread data (in m/min).a Chart of the fire behavior of fuel situation Ea-05 from the photo-guide.b Chart of the fire behavior of Rothermel fuel model 4 obtained with VisualFuego (LABIF-UCO, 2019).c Table showing the results of the calculations performed in VisualFuego (LABIF-UCO, 2019) for Rothermel fuel model 4

Fig. 4
Fig. 4 Land cover map of the whole study area

Fig. 5
Fig. 5 Detailed view of the CHM in a forest stand (left) and the corresponding PNOA image 2020 (MTMAU, 2023) (right)

Fig. 6 Fig. 8
Fig. 6 Raster layers representing the number of points in three different horizontal sections of the LiDAR point cloud in a forest stand that corresponds to fuel model 7; the corresponding vertical section of the LiDAR point cloud for the stand is shown on the right

Fig. 9
Fig. 9 Example of the main predictive LiDAR variables in a forest stand corresponding to fuel model 10.The variables are shown in order of importance according to their mean decrease in Gini: a entropy, b cumulative percentage of returns in the 5th layer, c cumulative percentage of returns in the 6th layer, d cumulative percentage of returns in the 4th layer.As additional context, e shows the corresponding discontinuous vertical section of the LiDAR point cloud

Fig. 10
Fig. 10 Graphic examples of the defined fuel models.a Fuel model 1: Rocky areas with low grass.b Fuel model 2: areas with crops.c Fuel model 4: high shrub areas.d Fuel model 5: low shrub areas.e Fuel model 7: tree-covered areas with continuous vertical structure.f Fuel Model 10: tree-covered areas with discontinuous vertical structure.g Fuel model 12: harvested areas the characterization of the vertical arrangement of vegetation at this resolution.In contrast,Huesca et al. (2019) produced a map at 30 m,Marino et al. (2016) at 25 m andDomingo et al. (2020) at 10 m.According toTaneja et al. (2021), different spatial resolutions of fuel models provide different results in fire simulations.They found that when using a fuel model map with higher resolutions in simulation, better prevention strategies might be established.Nevertheless, spatial resolution is limited by the resolution of the ALS data.The commented studies used point clouds with a density of 0.5 pts/m 2 and 1 pt/ m 2(Huesca et al. 2019;Marino et al. 2016;Domingo et al. 2020), while the ALS resolution for this study was much greater: 7.2 pts/m 2 .

Fig. 12
Fig. 12 Plot of the mean decrease in Gini for each variable

Table 2
Legend of the land cover map created using Sentinel-2 images ShrubsAreas of shrubs (mainly Ulex and Cytisus genera and ericaceae family)Rocky areasLand covered by rocks and very small shrubs or non-anthropogenic, non-vegetated areas Anthropic Built-up areas Water Bodies of water, both marine and continental

Table 3
Metrics obtained from the normalized LiDAR point clouds

Table 4
Height thresholds of the horizontal sections of the point cloud in meters

Table 6
Fuel model correspondences between the fuels photoguide, the Rothermel fuel models, and the remote sensing land cover classes

Table 7
Area occupied by each land cover class in the land cover map obtained

Table 8
Confusion matrix of land cover classification Euc Eucalyptus sp, Con Conifers, Bro Broadleaves, Cro Crops and pastures, Shru Shrubs, Rock Rocky areas, Ant Anthropogenic areas, Wat Water, T Total, UA User's Accuracy, PA Producer´s Accuracy, OA Overall Accuracy, KI Kappa Index

Table 9
Area occupied by each fuel model

Table 10
Confusion matrix for fuel models 7 and 10.The numbers are the number of pixels analyzed

Table 12
(Arellano et al. 2017) to be present in the study area Description of Lourizán CIF(Arellano et al. 2017)fuel situations detected in the study area