Developing customized fuel models for shrub and bracken communities in Galicia (NW Spain)

Journal of Environmental Management 351 (2024) 119831 Available online 21 December 2023 0301-4797/© 2023 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Research article Developing customized fuel models for shrub and bracken communities in Galicia (NW Spain) Jos´ e A. Vega a , 1 , Juan Gabriel ´ Alvarez-Gonz´ alez b , 1 , St´ efano Arellano-P´ erez b , c , Cristina Fern´ andez a , Pedro Cui˜ nas a , Enrique Jim´ enez a , Jos´ e M. Fern´ andez-Alonso a , Teresa Fontúrbel a , Cecilia Alonso-Rego b , Ana Daría Ruiz-Gonz´ alez b , * , 1 a Centro de Investigaci´ on Forestal de Louriz´ an, PO Box 127. 36080, Pontevedra, Spain b Unidad de Gesti´ on Ambiental y Forestal Sostenible (UXAFORES). Departamento de Ingeniería Agroforestal, Campus Terra, Universidad de Santiago de Compostela, Campus Universitario S/n, 27002, Lugo, Spain c AGRESTA Sociedad Cooperativa, C/ Duque de Fern´ an Nu˜ nez 2, 28012, Madrid, Spain ARTICLE INFO Keywords: Wildland fuels Fuel models Medoids Cluster analysis Fire behaviour Fire management Fuels classification ABSTRACT Geospatial fire behaviour and fire hazard simulators, fire effects models and smoke emission software commonly use standard fuel models in order to simplify data collection and the inclusion of complex fuel scenarios. These fuel models are often mapped using remotely sensed data. However, given the great complexity of fuelbeds, with properties that vary widely in both time and space, the use of these standard fuel models can greatly limit accurate fuel mapping. This affects fuel hazard assessment, fuel reduction treatment plans, fire management decision-making and evaluation of the environmental impact of wildfire. In this study, we developed unique customized fire behaviour fuel models for shrub and bracken communities, by using k-medoids clustering analysis based on both fuel structural characteristics and potential fire behaviour. We used an original database of 722 destructive sample plots in nine different shrub and bracken communities covering the entire distribution area in Galicia (NW Spain), one of the regions in Europe most affected by forest fires. Measurements of cover, height and fuel fractions loads differentiated by size and vegetative state (live or dead) were used to estimate the potential rate of fire spread with five different models including fireline intensity, heat per unit area and the flame length for each sampling site and considering extreme environmental conditions. The optimal number of clusters was established by combining practical knowledge about the shrubland communities under study and their associated fire behaviour, with maximization of the mean value of the silhouette variable and minimization of the within-cluster sum of squares. The structural characteristics of the medoids derived from the analysis were associated with each of the proposed customized fuel models. Finally, a simple dichotomous classification based only on shrub height was developed to enable construction of spatially explicit fuel model maps based on remotely sensed data. Thus, the methodology applied allows generation of a more realistic representation of fuel distribution in the landscape, based on fuel structure measurements of natural regional ecosystems rather than on the use of standard models. We believe that the proposed methodology is generally applicable to communities composed of other shrub and fern species in different biogeographical regions. 1. Introduction Detailed and accurate wildland fuel parameterization across landscapes is an essential support in fire hazard assessment and fire management decision-making. In particular, fire behaviour prediction, a critical factor in wildfire management, requires adequate spatial characterization of vegetation in terms of fuel properties (Finney, 1998, 2006). Land-use planning, prioritisation of fuel reduction treatments * Corresponding author. E-mail addresses: [email protected] (J.A. Vega), [email protected] (J.G. ´ Alvarez-Gonz´ alez), [email protected] (S. Arellano-P´ erez), [email protected] (C. Fern´ andez), [email protected] (P. Cui˜ nas), [email protected] (E. Jim´ enez), josemfernandez@uvigo. gal (J.M. Fern´ andez-Alonso), [email protected] (T. Fontúrbel), [email protected] (C. Alonso-Rego), [email protected] (A.D. RuizGonz´ alez). 1 These authors contributed equally to this work. Contents lists available at ScienceDirect Journal of Environmental Management journal homepage: www.elsevier.com/locate/jenvman https://doi.org/10.1016/j.jenvman.2023.119831 Received 13 July 2023; Received in revised form 29 November 2023; Accepted 3 December 2023 Journal of Environmental Management 351 (2024) 119831 2 and assessment of the environmental impact of wildfire can also benefit from such information. Given the great complexity of fuelbeds, with properties that vary widely in both time and space, numerous forest fuel classification systems have been devised (e.g. Sandberg et al., 2007; Ottmar et al., 2007; Arroyo et al., 2008; Fernandes et al., 2009b; Weise and Wright, 2014; Duff et al., 2017; Gale et al., 2021; Rego et al., 2021; S´ a et al., 2023). In particular, Keane (2013, 2015) analyzed and discussed the advantages and limitations of the main approaches used to create fuel classifications for fire behaviour models: association, opportunistic, classification and abstraction. In the first approach, fuel information is linked to the different categories of vegetation-based classifications. This method is used in the Canadian Forest Fire Behaviour Prediction System (Forestry Canada Fire Danger Group, 1992), in which fuel types are associated with the main vegetation cover classes in the region (Phelps and Beverly, 2022), although shrublands are not considered. In Australia, the Bushfire Fuel Classification (BFC) uses an associative method, with a top-down hierarchical approach (Gould and Cruz, 2012; Hollis et al., 2015). The first level corresponds to major vegetation types, including shrublands, and the second level is defined by biophysical fuel properties, which are direct inputs of specific fire behaviour models (Cruz et al., 2018). In the opportunistic method, unique fuelbeds are subjectively identified in the field and selected as a new category for inclusion in the classification on the basis of their representativeness. The Natural Fuels Photo Series (Ottmar et al., 1998, 2004, 2007) and the Fuel Characteristics Classification System (FCCS) (Riccardi et al., 2007a, 2007b; Sandberg et al., 2007; Prichard et al., 2013) are examples of this approach that are not restricted to fire behaviour prediction but also assess fire effects, fuel consumption, smoke production and carbon storage. The classification approach is based on the systematic, comprehensive clustering of fuelbeds into unique groups based on selected attributes, mainly loads quantified by fuel components (e.g. Dimitrakopoulos, 2002; Miller et al., 2003; Parresol et al., 2012; Ch´ avez et al., 2014; Bright et al., 2016). Finally, the abstraction approach relates fuelbed attributes to fire behaviour. This leads to the concept of a fuel model as a simplified quantitative representation of a fuelbed, constituted of a set of physical parameters used as inputs for fire behaviour models (Rothermel, 1972; Burgan and Rothermel, 1984; Andrews and Queen, 2001). The four shrub fuel models developed by the Northern Forest Fire Laboratory (NFFL) in the United States (Albini, 1976a; Anderson, 1982) and later expanded to nine models in the Fire Behaviour Fuel Models (FBFM) (Scott and Burgan, 2005) are the most commonly used in this approach. These were developed as categorical representations of fuel attributes for use as inputs in the Rothermel (1972) surface fire spread model. Customized fuel models try to represent local conditions more accurately (Andrews, 1986; Mallinis et al., 2008), and they are created from fuel data and inventories or by adjusting existing models (Burgan and Rothermel, 1984; Burgan, 1987; Ascoli et al., 2015). Customized shrub models are particularly required for Mediterranean ecosystems, due to the high structural heterogeneity and complexity of these plant communities (Arca et al., 2007), and cluster analysis has frequently been used for this purpose (e.g. Dimitrakopoulos, 2002; Fernandes et al., 2009b; Duce et al., 2012; Elia et al., 2015). These models provide the inputs for the most widely used fire behaviour modelling systems in use today at operational level, such as BehavePlus (Andrews et al., 2008), Farsite (Finney, 1998) and FlamMap (Finney, 2006), all of which are based on the equations developed by Rothermel (1972). Other hierarchical approaches used in fuel classification, frequently supported by multispectral imagery and LiDAR data, and in which vegetation is categorized and fuelbeds are identified, have led to the fuel type concept. This term is assumed to represent vegetation classes with similar fire behaviour (Pyne et al., 1996; CIFFC Canadian Interagency Forest Fire Centre, 2003; Scott and Burgan, 2005). However, the basis of this assumption has not been fully developed and, in fact, different vegetation types can display similar fire behaviour characteristics. According to this interpretation, fuel types (considered “general fire-carrying fuel categories”) include different fuel models (Scott and Burgan, 2005). The term fuel type is also frequently used as an equivalent of fuel model, although it provides more limited information on the physical properties than provided by fire behaviour fuel models. One such example is the Prometheus system (Prometheus, 2000), which considers three shrubland fuel types, according to their height ranges. Although no quantitative justification is given for selection of these ranges in terms of fire behaviour (Ria˜ no et al., 2002), and no other fuel physical characteristic is provided, this fuel typology is widely used in Mediterranean fuels mapping (e.g. Lasaponara and Lanorte, 2007; Huesca et al., 2019; Viedma et al., 2020). The proposed FirEurope Fire Risk Map (Aragoneses et al., 2023) also identifies three shrub fuel types based on height ranges. These authors suggest subjectively linking these fuel types to six FBFM shrub fuel models, depending on biogeographical features (arid or humid) and directly assuming the physical parameters of FBFM. However, these authors also question the appropriateness of this approach, in line with others (e.g. Santoni et al., 2011; Salis et al., 2016). Finally, both association and abstraction approaches require assignment of the particular vegetation scenario under consideration to a fuel model. This process generally requires decision rules, particularly for constructing customised models. Overall, the assignment is subject to various sources of error and uncertainty, such as subjective perception of potential fire behaviour in the fuel complex and spatial variation of fuel models, among others (Jakubowksi et al., 2013; Perrakis et al., 2018; S´ a et al., 2023). It is therefore strongly recommended that the process should rely on robust classification methods, supported by fuel structure variables obtained in sufficiently representative fuel inventories (Ascoli et al., 2020) and linked to fire behaviour prediction models. Finally, these models should be tested with real fire behaviour data or expert judgement should at least be applied. Obviously, the whole task is challenging, and the current fuel classification systems do not generally meet all of these requirements. Galicia is a particularly interesting area for developing customized shrub models because 29% of all forest fires in Spain occur in this region, and the burnt area represents 22% of the total (MAPA, 2019), even though forest land in the region represents only 10% of the Spanish forest area. Furthermore, although the mean area burnt annually in Galicia over the last 10 years has decreased to around 16,000 ha (Consellería de Medio Rural, 2022), the potential for catastrophic wildfires in the region is highlighted by some examples of severe fire seasons, with about 62,000 and 52,000 ha burnt in 2017 and 2022 respectively, in line with the unfavourable climate change projections for the area (Vega et al., 2009a, 2009b). The problem particularly affects shrublands, as 63% of the burned area in the region occurs in these communities (Consellería de Medio Rural, 2022), even though shrubland represent only about 30% of the total forest land area (MARM, 2011b). The fast shrub growth, associated with mild annual temperatures and abundant rainfall in most of the region (Martinez Cortizas and P´ erez Alberti, 1999), promotes large fuel accumulations in shrublands (Vega et al., 2022a), favouring high intensity wildfires (Arellano-P´ erez et al., 2017). The location of Galicia between the Atlantic and Mediterranean biogeographical domains (European Environment Agency, 2016) generates a wide variety of shrubland floristic composition with its associated broad range in fuel attributes. Currently, the Forest Fire Service of Galicia is using NFFL models for the characterization of shrublands. Nonetheless, in a previous study based on destructive sampling in Galicia (Vega et al., 1998), fire spread rates according to the Rothermel (1972) model mismatched those observed in experimental fires in the field when using customized shrub fuel models obtained following the methodology proposed by Burgan and Rothermel (1984). This suggests the need for further research to develop more suitable fuel models for fire management operational uses. For the above reasons, a shrub and bracken fuel classification method that allows unique fuel categories to be differentiated on the basis of selected attributes, and in association with their potential fire J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 3 behaviour is essential for making key decisions about the best alternatives for fuel and fire management. Accordingly, the aim of this study was to develop customized fuel models for shrub and bracken communities of Galicia on the basis of both structural characteristics of the fuelbeds and potential associated fire behaviour variables. 2. Material and methods 2.1. Study area, shrub communities considered, inventory plots and fuel structure variables The study was conducted in Galicia (NW Spain), where shrublands cover an area of just over 600,000 ha (approximately 20% of the total area and 30% of the forest area of the region). Communities dominated by gorse (Ulex sp.) occupy 55% of the shrubland area, heather-type shrubs (Ericacea sp.) cover 22% of the area, and brooms (Cytisus sp. or Pterospartum tridentatum) occupy 11%. The remaining treeless land is distributed between areas with little or no vegetation and high fire recurrence, areas dominated by brackens, rockroses (Cistus sp. and Halimium sp.) and grasslands (MARM, 2011a, 2011b). The climate is transitional, with a large area in the west being humid and temperate with oceanic influence, shifting to Mediterranean climate eastwards and inland (Retuerto and Carballeira, 1992), with slight continental features regulated by the complex relief (Rodríguez-Guiti´ an and Ramil-Rego, 2007). The mean annual precipitation is 1200 mm, varying between about 500 and 1800 mm; the mean annual temperature is 13.3 ◦C, ranging seasonally between 8.5 and 19 ◦C (Martinez-Cortizas and P´ erez-Alberti, 1999). The most frequent soils on which the shrublands are established are Regosols and Leptosols and to a lesser extent Cambisols and Umbrisols developed mainly on granitic and schist substrates (Carballas et al., 2016). In the present study, eight shrub communities were initially considered: gum rockrose, dominated by Cistus ladanifer (Cl); low broom, dominated by Cytisus multiflorus (Cm); high broom, with Cytisus striatus, C. scoparius as the main species; high heath, dominated by Erica australis, Erica arborea or E. scoparia (Ea); low heath, with Erica umbellata or E. mackaiana (Eu); prickled broom, dominated by Pterospartum tridentatum (Pt); high gorse, formed by pure stands of Ulex europaeus (Eu) or mixed stands, with Ericacea species; and low gorse, dominated by Ulex gallii or U. minor (Ug). These communities cover about 90% of the Galician shrublands (Izco et al., 1999) and are mainly composed of perennial, multi-stemmed evergreen sclerophyll woody species which usually form closed stands of medium to moderately high height (0.5–3 m). Bracken communities consisting of Pteridium aquilinum (Pa) and generally including a mix of herbaceous and to a lesser extent woody species, were also considered in this study. This type of community shows stronger seasonal variations in the distribution of live and dead biomass than shrub communities, forming temporally extensive patches of vegetation in the area and frequently being involved in fires in the region. A more complete description of each community is provided in Vega et al. (2022a). These communities were destructively sampled to estimate the main structural characteristics. A circular plot of variable radius (20–30 m in length, depending on the height of the shrubs) was established in each sampling site. Two perpendicular diameters were randomly selected and four-square sampling subplots were established in the centre of the four radii of the circular plot. The area of each subplot varied between 4 and 36 m 2 , depending on shrub height (2 ×2 m for shrubs smaller than 1.0 m in height and 3 ×3 or 6 ×6 m for shrubs taller than 1 m). Destructive sampling was carried out in each subplot. The material in the standing shrub stratum was physically separated by size class into fine fuels (diameter<0.6 cm, hereafter G1), medium fuels (0.6 cm ≤diameter <2.5 cm, hereafter G2) and coarse fuels (2.5 cm ≤diameter <7.5 cm, hereafter G3), and it was also further subdivided by vegetative condition (live or dead). Finally, the values of the following structural characteristics were obtained in each of the 722 inventories: shrub cover (Cov Shr ); mean shrub height (hShr); dead fine shrub load (W Shr_G1_dead ); live fine shrub load (W Shr_G1_live ); fine shrub load (W Shr_G1 ); coarse shrub load (W Shr_G23 ) and total shrub load (W Shr =W Shr_G1 +W Shr_G23 ). The destructive sampling technique is described in detail by Vega et al. (2022a)(Vega et al., 2022b) The basic descriptive statistics of shrub strata for the main structural characteristics of each shrub community are shown in Table 1. 2.2. Variables related to surface fire behaviour As the aim of the study was to develop customized fuel models based on both structural characteristics of the shrubland and associated fire behaviour variables, the potential rate of fire spread (R), fireline intensity (I), heat per unit area (H) and flame length (fl), were estimated for each shrubland sampling site. In determining the fire behaviour variables, a distinction was made between inventories of shrubland communities (674 sample plots) and inventories of bracken-dominated communities (49 sample plots) because of the strong seasonal variations in biomass condition (live or dead) of the bracken-dominated communities and the different limitations of the models for estimating surface fire rate of spread in both cases. For shrubland, the potential rate of fire spread (R) was estimated using five different equations to develop the customized fuel models: three empirical equations, a semi-physical approach and a physical approach. 1) The empirical model proposed by Anderson et al. (2015), developed using shrubland fire behaviour data from 79 experimental fires. These burns covered a wide range of heathlands and shrubland species associations and vegetation structures in Australia, New Zealand, South Africa, Portugal and Spain. The model was evaluated using data from 32 fire spread rates observations in 24 wildfires and 32 prescribed fires in mixed heathlands in the study area: R=5.6715 ⋅ (wf⋅U10)0.9102 ⋅h0.2227 wgh ⋅exp(− 0.0762 ⋅ Md)(1) where R is the rate of spread (m/min); w f is the wind correction factor used to estimate the wind speed at midflame from the wind speed at 10 m (U 10 , km/h), which for these communities, according to the authors, takes a value of 0.67; h wgh is the height of the shrubland weighted by shrubland cover (m), and M d is the moisture content of the dead fine fuels (%). 2) The empirical model proposed by Vega et al. (1998), developed using data from 41 experimental fires carried out on the same shrubland communities considered in this study: R=0.249 ⋅ (wf⋅U10 /3.6)1.193 ⋅(hwgh⋅100)0.658⋅exp(0.01088 ⋅ S)(2) where S is the slope expressed in percent and the other parameters are as previously defined. 3) The empirical model proposed by Fernandes (2001), based on data from 29 experimental fires in shrubland and six prescribed burns for fuel hazard reduction in Portugal on similar shrublands communities considered in this study: R=7.255 ⋅ exp(0.092⋅wf⋅U10)⋅h0.932 wgh ⋅exp(− 0.067 ⋅ Md)(3) 4) The Balbi physical fire spread model (Balbi et al., 2020; Chatelon et al., 2022) with the values proposed by the authors for the specific parameters of the model (see Table 1 in Chatelon et al., 2022). Because the equation used to estimate the surface fire rate of spread (R) depends on itself, among other factors, the R values were estimated by iterative procedures. The model was evaluated with 109 J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 4 data from experimental and prescribed fires in shrublands (Anderson et al., 2015) and 14 experimental fires in fynbos vegetation (van Wilgen et al., 1985). 5) The semi-physical approach proposed by Rothermel (1972), with the modifications proposed by Frandsen (1973), Albini (1976b) and Andrews et al. (2013). Structural data from the shrub inventories related to the fuel load sorted by size and physiological status (live or dead) classes and fuel depth were used to construct customized fuel models as inputs for the fire spread model following Burgan and Rothermel (1984), and the “ros” function of the “Rothermel” package (Vacchiano and Ascoli, 2015) for R (R Core Team, 2020) was used. The Rothermel model requires several additional parameters such as the moisture of extinction (M x , %), the low heat content for the fuel classes (h, kJ/kg), and the surface-to-volume ratio for the fuel classes ( σ , m 2 /m 3 ). Furthermore, default constant values for the total mineral content fraction (0.055), effective mineral content (0.010) and oven-dry particle density (513 kg/m 3 ) were used. M x was set at 40% according to the values proposed by Scott and Burgan (2005) for subhumid to humid climate; the low heat content (h) was set at 18,000 kJ/kg for all of the shrubland communities and fuel classes according to Keane (2015), and the values of the surface-to-volume ratio ( σ , m 2 /m 3 ) were set at those proposed by Hernando et al. (2004) for the shrubland communities in the study area. Moreover, the midflame wind speed (U), required as input in this (semi-physical) model, was estimated from U 10 using the wind adjustment factor (Waf) proposed by Andrews (2012) for formations with no tree cover or fraction of covered less than 5%, as in this study. For bracken-dominated communities, only Rothermel’s model was used, as there were no empirical models for these formations, and the Balbi physical model (Balbi et al., 2020; Chatelon et al., 2022) only applies to shrublands. Most of the values of the input parameters used for these communities were the same as those used for the shrubland communities. The exceptions were for the surface-to-volume ratio, which was assigned the value proposed by Hernando et al. (2004) for the same community, and the moisture of extinction (M x ), which was assigned a value of 25%, as this is the mean value proposed by Scott and Burgan (2005) for conditions that are potentially similar to those in bracken-dominated communities. Since the morphological and physiological characteristics of these bracken-dominated communities change throughout the seasons, four different situations were considered: 1) full vegetative development (similar to those at the time of the inventory); 2) autumn conditions, during which the bracken is dead but standing at the same height as when sampled (the total load is considered dead fuel); 3) winter conditions, during which bracken is dead and moderately compacted (total load is considered dead fuel and mean height is adjusted to 40 cm according to field observations); and 4) conditions of full vegetative development but assuming that only the dead fuel burns, due to the high moisture contents of the live fine fuels and considering a mean height for this fuel of 40 cm according to field observations (Ruiz-Gonz´ alez, 2022). Extreme burning conditions were assumed for simulation of fire behaviour with all models, although within the range of U 10 and M d values for which both the fuel effect and the weather effect are relevant in fire behaviour (Cruz et al., 2022). A wind speed of 28 km/h at 10 m (U 10 ) was considered. This value corresponds to the 97% percentile of the local wind speeds observed at 14 h local time during the fire season (122 days from June to September) at the meteorological stations in Santiago (42◦53 ′ 17 ″ N; 8◦24 ′ 38 ″ W), Lugo (43◦6 ′ 41 ″ N; 7◦27 ′ 27 ″ W), Ourense (42◦19 ′ 31 ″ N; 7◦51 ′ 35 ″ W) and Vigo (42◦14 ′ 19 ″ N; 8◦37 ′ 26 ″ W), as the most representative of the study area and for which the largest time series of data are available (see Arellano-P´ erez et al., 2020). The moisture content of dead fine fuels and live fine fuels were assumed to be respectively 7% and 90%, which are the minimum values observed in different shrubland communities in the study area (Vega et al., 2009c; Ruiz-Gonz´ alez, 2022). For bracken communities values of 12 and 150% were assigned to the moisture content of the dead and live fine fuels, respectively, according to the field values for these communities (Ruiz-Gonz´ alez, 2022). Finally, the fireline intensity (I, kW/m), the heat per unit area (H, kJ/ m 2 ) and the flame length (fl, m) were estimated using the equations of Byram (1959) and considering only fine fuels as those consumed in the active flaming front. 2.3. Statistical analysis The non-hierarchical k-medoids clustering method (Kaufman and Rousseeuw, 1987) was used to develop the new customized fuel models. K-medoids clustering is a variant of k-means clustering. K-means clustering iteratively finds the k centroids and assigns each observation to the nearest centroid, where the coordinate of each centroid is the mean of the coordinates of the observations in the cluster. By contrast, k-medoids clustering finds the k most centered observations, the so-called medoids. Because it is based on the most centrally located observation in a cluster, it is less sensitive to outliers than k-means clustering (Park and Jun 2009; Jin and Han, 2011). The partitioning around medoids (PAM) algorithm proposed by Kaufman and Rousseeuw (1990), was used to define the k clusters, with the absolute distance between observations as a dissimilarity measure. Table 1 Mean values of the standing shrub fuel strata characteristics. Std. dev. =standard deviation, n =number of plots, hShr =shrub height, Cov Shr =shrub cover, W Shr =total shrub fuel load, W Shr_G23 =coarse shrub fuel load, W Shr_G1 =fine shrub fuel load, W Shr_G1_dead =dead fine shrub fuel load, W Shr_G1_live =live fine shrub fuel load. Variable Statistic Cl Cm Cs Ea Eu Pa Pt Ue Ug n 23 47 44 125 68 49 69 191 106 hShr mean 120.04 114.49 241.70 110.40 51.75 105.24 90.54 115.09 74.99 (cm) std. dev. 38.33 52.91 157.43 76.30 19.09 31.63 52.13 60.08 27.77 Cov Shr mean 69.26 85.09 84.59 89.77 92.18 83.49 83.76 84.62 93.38 (%) std. dev. 14.90 16.29 17.26 17.38 14.46 13.43 17.47 22.74 14.40 W Shr mean 1.11 2.37 5.37 2.47 1.97 1.05 2.49 3.36 2.84 (kg m −2 ) std. dev. 0.34 1.06 3.60 1.84 0.88 0.52 1.30 1.45 0.89 W Shr_G23 mean 0.29 0.84 3.54 0.99 0.26 0.17 0.63 1.26 0.56 (kg m −2 ) std. dev. 0.18 0.78 3.18 1.23 0.26 0.14 0.69 1.06 0.54 W Shr_G1 mean 0.82 1.52 1.83 1.48 1.71 0.88 1.86 2.10 2.28 (kg m −2 ) std. dev. 0.20 0.34 0.61 0.77 0.71 0.40 0.74 0.68 0.65 W Shr_G1_dead mean 0.07 0.49 0.48 0.40 0.56 0.43 0.70 0.82 0.77 (kg m −2 ) std. dev. 0.04 0.36 0.25 0.28 0.33 0.32 0.29 0.37 0.39 W Shr_G1_live mean 0.74 1.04 1.36 1.07 1.15 0.45 1.16 1.28 1.51 (kg m −2 ) std. dev. 0.20 0.28 0.46 0.53 0.44 0.24 0.56 0.43 0.41 Cl =Cistus ladanifer, Cm =Cytisus multiflorus, Cs =Cytisus striatus, Ea =Erica australis, Eu =Erica umbellata, Pa =Pteridium aquilinum, Pt =Pterospartum tridentatum, Ue =Ulex europaeus and Ug =Ulex gallii. J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 5 All of the fuel complex structural variables and the associated fire behaviour variables were initially considered as possible classification features; however, analysis of the correlation matrix indicated (as expected) strong correlations between fl and R, H and I (r >0.80) and between I and R (r =0.87); fl and I were therefore excluded from the analysis to avoid multicollinearity problems that could affect the results. On the other hand, W Shr_G1live and W Shr_G23 were also excluded from the analysis owing to the additivity relationships between the fuel loads of the different fractions (W Shr_G1 =W Shr_G1_dead +W Shr_G1_live and W Shr = W Shr_G1 +W Shr_G23 ). Variables remaining were standardized to account for differences in means and variances (Elia et al., 2015; Alhaj-Khalaf et al., 2021). The number of clusters was optimized by maximizing the silhouette values (Rousseeuw, 1987) and minimizing the sum of squares of within-cluster distances, i.e. the so called “Elbow method” (Ketchen and Shook, 1996). Finally, flexible discriminant analysis was performed to assess the accuracy of the cluster analysis. Analyses were performed using R software (R Core Team, 2020). The “pam” function of the “cluster” package (Maechler et al., 2019) was used to classify the observations, the “fviz_nbclust” function of the “factoextra” package (Kassambara and Mundt, 2020) was used to optimize the number of clusters, and the “fda” function of the “mda” package (Leisch et al., 2022) was used to run discriminant analysis. 3. Results and discussion 3.1. Fuel models for shrubland communities Analysis of the correlation between the values of the fire rate of spread (R) obtained with each of the five equations proposed for the extreme burning conditions evaluated indicated significant ( α =0.01) and positive correlations for all comparisons, especially for the three empirical models (Anderson et al., 2015; Vega et al., 1998; Fernandes, 2001), with Pearson’s linear correlation ( ρ ) values above 0.95 (see Fig. S-1 in supplementary material). The equation proposed by Anderson et al. (2015) yielded the highest mean value of R (47.13 m/min) and, at the same time, the lowest variability (s.d. =6.64 m/min). The model proposed by Rothermel (1972) yielded the lowest mean value (23.57 m/min) but the highest variability (s.d. =24.29 m/min). Comparison of only the two “local” models (Vega et al., 1998; Fernandes, 2001) (which are based on studies of fire behaviour in shrub communities similar to those considered in the present study and which also share biogeographical characteristics) showed that the estimated mean R values are significantly different ( α = 0.05; see Fig. S-2 in supplementary material), although the ρ value is 0.99 and the standard deviations are very similar: 17.71 m/min for Fernandes (2001) and 16.34 m/min for Vega et al. (1998). However, the estimated R values correspond to values obtained in a single simulation for extreme burning conditions, so the comparisons are merely descriptive. The previous analysis for determining the optimum number of clusters produced the same result regardless of the modelling approach used to estimate R. According to the criterion of maximizing the mean value of the silhouette variable, the optimum number of clusters is 2, with a local maximum at 4, which stabilizes as the number of clusters increases. When considering the number of clusters that minimizes the within cluster sum of squares, the optimal number ranges from 7 to 12 clusters, although the pattern is clear stabilization from four clusters onwards. These results, combined with practical knowledge about the shrubland communities under consideration and the associated fire behaviour, led to the development of four customized fuel models for these formations. Fig. 1 shows the pattern of these two selection criteria for estimation of the fire rate of spread (R) using the model proposed by Vega et al. (1998), considered the most appropriate since it is based on fire behaviour studies carried out in Galicia and in the same shrubland communities as those analyzed in this study. Once the optimal number of clusters was determined, the four clusters obtained for each of the five rate of fire spread (R) estimation approaches were compared (Table 2). A total of 67.16% of the plots were classified similarly, regardless of the R estimation method used (Table 3). When considering only the empirical models (Anderson et al., 2015; Vega et al., 1998; Fernandes, 2001) the percentage agreement was 87.37%. Although these three empirical models have a similar structure and common independent variables, the classification based on the model proposed by Anderson et al. (2015) yields discrepancies of about 12% in the classification relative to the other two empirical models. This contrasts with the 1.2% discrepancy between the two “local” models (Vega et al., 1998; Fernandes, 2001). On the other hand, the empirical model of Anderson et al. (2015) also yields the largest discrepancies relative to the semi-physical (>25%) and physical approaches (>22%), whereas the discrepancies between the other two empirical models and the semi-physical or physical approaches do not exceed 18% or 15%, respectively. This can be explained by the lower variability in the estimated R values obtained with the empirical model of Anderson et al. (2015) than in the values obtained with the remaining approaches (see Fig. S-2), especially the other two empirical models and the semi-physical model of Rothermel (1972). Table 3 shows the percentage concordance for the five R estimation approaches, for the three empirical models and for the two local models (Vega et al., 1998; Fernandes, 2001) applied to different shrubland communities. The greatest discordance occurs in the classification of sample plots of low broom-dominated communities (Cm) or of low heath-dominated communities (Eu) and, to a lesser extent, in sample plots of low gorse-dominated communities (Ug). The percentages of sample plots classified with negative silhouette values, i.e. sample plots in which the mean distance from sample plots in some other cluster is less than the mean distance from sample plots in Fig. 1. Graphs for determining the optimal number of clusters for shrubland communities by maximizing the mean silhouette value and minimizing the within cluster sum of square of distances when the rate of fire spread (R) is estimated using the model proposed by Vega et al. (1998). Table 2 Percentage concordance between the classifications based on the different approaches for estimating the fire rate of spread in shrubland communities. R estimation approach Vega et al. (1998) Fernandes (2001) Balbi et al. (2020) Rothermel (1972) Anderson et al. (2015) 87.73% 88.56% 77.56% 74.15% Vega et al. (1998) 98.81% 86.03% 82.32% Fernandes (2001) 85.74% 82.02% Balbi et al. (2020) 82.06% J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 6 the same cluster, were low in all cases, ranging from 5.1% for estimation of R with the model proposed by Vega et al. (1998) to 8.5% for estimation of R with the model proposed by Rothermel (1972). Considering the five models based on the different approaches for estimating the fire rate of spread overall, the medoids obtained in each cluster for each classification were similar. For both the first and fourth clusters there were only two possible medoids in each case, while for the second and third clusters there were three possible medoids. The results shown below were obtained in the cluster analysis of the sample plots in which the R was estimated using the empirical model proposed by Vega et al. (1998). Table 4 shows the results of the discriminant analysis of the four clusters obtained. The overall classification error rate was 9.09%, and the discriminant analysis indicates that the clusters are effectively unique groups that can be feasibly distinguished on the basis of structural variables and associated fire behaviour variables. The first two canonical variables (CV1 and CV2) explained 99.69% (86.88 and 12.81%, respectively) of the discriminatory capacity of the variables used in the cluster analysis. The structural variables with the greatest weighting in the formulation of the canonical variables were hShr, (13.5% in CV1 and 27.0% in CV2), W Shr_G1 (15.2% in CV1 and 23.5% in CV2) and W Shr_G1_dead (12.4% in CV1 ang 5.4% in CV2), whereas R was the fire behaviour variable with the greatest weighting in the formulations (26.3% in CV1 and 22.5% in CV2). Fig. 2 shows a scatter plot of the clusters using the first two canonical variables. The spatial arrangement of each cluster was clearly defined and there were only some slight overlaps corresponding to the error rate of 9.01%. The positions of the medoids obtained in the final cluster analysis, as well as the alternative medoids obtained in the other four cluster analyses using the other two empirical models (Anderson et al., 2015; Fernandes, 2001), the semi-physical approach (Rothermel, 1972) and the physical approach (Balbi et al., 2020) for R estimation, are shown in this scatter plot. As previously mentioned, the distance in canonical space between the medoids finally selected and the alternative options in each cluster is narrow, especially for clusters 1, 2 and 4. The first customized fuel model (Shrub-1; Fig. S-3), with the lowest height and lowest fuel loads; includes 197 sample plots (29.3% of the total sample plots) and corresponds to young communities with dominance of any of the species considered or non-senescent communities dominated predominantly by the shortest species in the study area, such as Erica umbellata or E. mackaiana (Ea) or Cistus ladanifer (Cl). In any case, this customized fuel model should be assigned to fuelbeds in which coarse fuel (G23) and dead fine fuel (G1dead) are poorly represented, with a predominance of the fine live fraction. The second customized fuel model (Shrub-2; Fig. S-4) includes 198 sample plots (29.4% of the total sample plots) and corresponds to shrub communities with relatively small mean shrub heights, although higher than the first model, but with much larger loads, especially of fine fuels, which results in larger fine fuel bulk densities. Also noteworthy is the increase in the dead fine fraction relative to Shrub-1, as Shrub-2 generally corresponds to adult communities whose leaves are dying but still hanging on the plants. The third customized fuel model (Shrub-3; Fig. S-5) includes 191 sample plots (28.4% of the total). This model has higher heights and fuel loads than the two previous ones and is the model with the highest load of both live and dead fine fuels, the latter representing about 40% of the total fine fuel load. About 67% of the sample plots (128) assigned to Shrub-3 correspond to communities dominated by gorse (Ulex europaeus, U. gallii or U. minor), although this model can be associated with any of the shrubland communities, except those of lower height in the adult stage, such as those dominated by Erica umbellata or E. mackaiana (Eu). Finally, the fourth customized fuel model (Shrub-4; Fig. S-6) includes 87 sample plots (12.9%) and corresponds to adult communities mainly dominated by species of the genus Cytisus, (Cs, Cm), Erica australis or E. arborea (Ea) and Ulex europaeus (Ue), which have the highest heights and largest total (W Shr ) and coarse fuel loads (W Shr_G23 ). In these formations, dead leaves and thin branches are frequently detached, especially from the lower and middle parts of the vertical profile, contributing to increasing the depth and load of the litter and duff layers, while generating gaps favouring wind circulation and thus fire propagation. According to the criteria reported by (Rodríguez y Silva et al. (2014; 2020), the proposed models represent an increasing gradient of fire suppression effort (Shrub-1 to Shrub-4) due to structural characteristics that affect the associated fire behaviour, the penetrability for firefighters and the fireline production rate. Fig. 3 illustrates this behaviour and shows the box plots with the distribution of R values estimated using the Table 3 Percentage concordance for the classifications based on the different approaches for estimating the fire rate of spread for the different shrubland communities analyzed. Shrubland community Percentage concordance All approaches 3 Empirical models 2 Local models Cl 91.30% 91.30% 100% Cm 51.06% 72.34% 97.87% Cs 79.55% 93.18% 100% Ea 77.60% 91.20% 99.20% Eu 48.53% 80.88% 98.53% Pt 66.67% 92.75% 98.55% Ue 68.59% 89.53% 98.95% Ug 61.32% 83.02% 98.11% Total 67.16% 87.37% 98.81% Table 4 Discriminant analysis results with the number of sample plots (top) and the percentage classified (bottom) in each cluster and the total percentage of classification error when using the empirical model proposed by Vega et al. (1998) to estimate the R values. From cluster To cluster 1 2 3 4 Total 1 180 17 0 0 197 (91.37%) (8.63%) (0%) (0%) 2 6 184 8 0 198 (3.03%) (92.93%) (4.04) (0%) 3 0 11 177 3 191 (0%) (5.76%) (92.67%) (1.57%) 4 0 0 16 71 87 (0%) (0%) (18.39%) (81.61%) Error 9.09% Fig. 2. Scatter plot of the four clusters obtained using the empirical model of Vega et al. (1998) to estimate the values of R. The positions of the medoids obtained in this cluster analysis (“Final Medoids”) as well as the medoids obtained in the four cluster analyses using the other approaches for the estimation of R (“Alternative Medoids”) are also shown. J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 7 equation proposed by Vega et al. (1998) included in each of the four proposed customized fuel models. In addition, values corresponding to the standard NFFL fuel models -FM4, FM5 and FM7- (Rothermel, 1972; Anderson, 1982) and the climate humid shrub models -SH3, SH4, SH6, SH8 and SH9- (Scott and Burgan, 2005) have been added. The FM6 model and the dry climate shrub models of Scott and Burgan (2005) were not included in the comparison, since the fuel situations described by these models are scarce in the study area. Fig. 3 shows how the NFFL and FBFM models analyzed do not adequately represent the extreme R situations simulated in the study area, which, however, are characterized by the Shrub-1 and Shrub-4 models proposed in this study. This could have serious consequences particularly due to the underestimation that occurs when using the attributes of the NFFL and FBFM models in the situations associated with the Shrub-4 model. Table 5 shows the mean value and standard deviation of the main structural variables of the shrubland in the sample plots included in each cluster as well as the values that correspond to the medoid of each cluster and which represent the proposed customized fuel model. The number of fuel models proposed seems consistent with the relatively wide range of variability of the fuel attributes of the shrub communities analyzed. This variability was mainly associated with the fact that the study area included two biogeographic regions (European Environment Agency, 2016), with contrasting bioclimatic domains: Atlantic (temperate-sub-humid/humid) and Mediterranean (arid/semiarid), according to the criteria of Zepner et al. (2020). Furthermore, the heterogeneous physiography, long history of disturbance and land use (Manuel and Gil, 2001; Kaal et al., 2011; L´ opez-Merino et al., 2012) and the wide floristic diversity of the shrub communities considered (Ramil et al., 2013) probably contribute to increasing this variability. The number of fuel models in studies classifying shrublands for fire behaviour and management purposes varies widely from one (e.g. Wu et al., 2011; Cai et al., 2014; Elia et al., 2015) to nine (Scott and Burgan, 2005) or to a figure open to the inclusion of new cases, as in the FCCS (Ottmar et al., 2007; Ottmar, 2014). The number we proposed was in the range of the frequently used fuel models/fuel types for shrublands (e.g. Anderson, 1982; Ria˜ no et al., 2002; Dimitrakopoulos, 2002; Fernandes et al., 2009b; Aragoneses et al., 2023). From a practical point of view, assignment of a customized fuel model to a surface fuelbed in the field should be based on structural variables that are easy to evaluate without carrying out complex inventories or measurements. Therefore, considering the most influential structural variables in the discriminant analysis and the ranges of these variables in the different customized fuel models, the mean shrub height (hShr) and the dead fine shrub load (W Shr_G1_dead ) were proposed as classifying variables. According to the results of the discriminant analysis, the overall error rate increased from 9.09% on the basis of structural variables and associated fire behaviour variables, to 19.62% on the basis of only hShr and W Shr_G1_dead . Fig. 4 shows the box plots with the distribution of the values of these two classificatory variables in each of the four customized fuel models proposed. The combination of both variables allows accurate classification with a limited degree of overlap between models. In addition, values corresponding to the standard NFFL fuel models -FM4, FM5 and FM7- (Rothermel, 1972; Anderson, 1982) and the climate humid shrub models -SH3, SH4, SH6, SH8 and SH9- (Scott and Burgan, 2005) have been added for comparison. The mean shrub heights associated with the NFFL and FBFM models analyzed do not adequately represent the extreme values observed in the sample plots (Fig. 4 upper), similar to the R values in Fig. 3. These results were expected since mean height is the fuel structural variable that drives the R estimates in the model of Vega et al. (1998). Regarding dead fine fuel loads (Fig. 4 lower), the opposite situation is observed, i.e., the NFFL models basically characterize extreme situations similar to the Shrub-3 and Shrub-1 models proposed in this study, without Fig. 3. Box plots of fire rate of spread (R) estimated with the model proposed by Vega et al. (1998) for the extreme burning conditions analyzed for each customized fuel model proposed for shrubland communities. Different letters indicate significant differences between mean values. The blue lines correspond to the R values estimated for the standard NFFL fuel models (FM4, FM5 and FM7) of Rothermel (1972) and Anderson (1982) and the red lines to the R values estimated for the humid climate shrub models (SH3, SH4, SH6 SH8 and SH9) of Scott and Burgan (2005). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Table 5 Medoid values and mean and standard deviation of the cluster for the main standing shrub fuel strata variables. Sd. =standard deviation, n =number of plots, hShr =shrub height, Cov Shr =shrub cover, W Shr =total shrub fuel load, W Shr_G23 =coarse shrub fuel load, W Shr_G1 =fine shrub fuel load, W Shr_G1_dead = dead fine shrub fuel load, W Shr_G1_live =live fine shrub fuel load. See definitions in the text. Different letters indicate significant differences between mean values ( α =0.95). Variable Cluster (Model) 1 (Shrub1) 2 (Shrub2) 3 (Shrub3) 4 (Shrub4) n 197 198 192 87 hShr (cm) Medoid 50.9 79.0 116.0 233.2 Mean (cluster) 53.3 d 82.2 c 117.9 b 264.1 a Sd. (cluster) 31.41 28.81 32.45 89.71 Cov Shr (%) Medoid 96.79 100 100 100 Mean (cluster) 73.63 b 90.66 a 94.37 a 93.81 a Sd. (cluster) 23.69 14.32 12.07 10.57 W Shr (kg m −2 ) Medoid 1.113 2.447 3.714 5.705 Mean (cluster) 1.109 d 2.394 c 3.687 b 6.127 a Sd. (cluster) 0.391 0.474 0.611 1.986 W Shr_G23 (kg m −2 ) Medoid 0.166 0.629 1.156 3.346 Mean (cluster) 0.159 d 0.529 c 1.171 b 3.765 a Sd. (cluster) 0.132 0.421 0.630 2.018 W Shr_G1 (kg m −2 ) Medoid 0.947 1.818 2.558 2.359 Mean (cluster) 0.951 d 1.865 c 2.516 a 2.362 b Sd. (cluster) 0.337 0.353 0.452 0.561 W Shr_G1_dead (kg m −2 ) Medoid 0.194 0.620 0.994 0.687 Mean (cluster) 0.208 d 0.611 c 0.977 a 0.824 b Sd. (cluster) 0.162 0.199 0.299 0.308 W Shr_G1_live (kg m −2 ) Medoid 0.753 1.198 1.564 1.672 Mean (cluster) 0.742 c 1.254 b 1.539 a 1.538 a Sd. (cluster) 0.279 0.340 0.396 0.426 J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 8 intermediate values that are represented by the proposed Shrub-2 and Shrub-4 models. In the case of fuel characterization using remote sensing, the dead fine shrub load is difficult to estimate accurately (e.g. Arellano-P´ erez et al., 2018; D’Este et al., 2021; Fern´ andez-Alonso et al., 2022), and dichotomous classification criterion based exclusively on mean shrub height was therefore developed (see Fig. 5). Considering this classification criterion, 67.2% of the 673 sample plots were correctly assigned to a customized fuel model; 19.4% of the sample plots were assigned to a customized fuel model of higher suppression complexity, and the remaining 13.4% were incorrectly classified in a customized fuel model with lower suppression complexity, the latter being the most unfavourable situation by underestimating the risk in fuel management decision-making. Selection of shrub height as a key classification variable for our shrublands is consistent with the results of previous studies on global sensitivity analysis of fire behaviour in different shrub communities. Such studies have identified height as one of the most influential fuel structural parameters affecting the potential rate of fire spread (R), in different prediction models. Regarding Rothermel’s model, Salvador et al. (2001) found that R was more sensitive to height than any other fuel or environmental variable and even similar to the wind effect, in scenarios with high wind velocity, in Mediterranean shrublands assimilated to NFFL fuel models 5 and 6. For these fuel models, Song et al. (2016), Liu et al. (2016), Ervilha et al. (2017), Ballester-Ripoll et al. (2018) and Plischke et al. (2021) also recognized height as the most influential fuel structural parameter. On the other hand, Ujjwal et al. (2021) found that R was more sensitive to the ratio of dead fine fuel load to total fuel load than fuel depth for the above-mentioned fuel models, and Cai et al. (2014) found that the fine fuel load was more important than fuel depth in a customized shrub model. The latter results are consistent with the discriminatory capacity shown by W Shr_G1_dead in the present study. In addition, the different sensitivity assessment methods used in the aforementioned studies may partly explain the apparent discrepancies regarding the most influential fuel parameter (Ujjwal et al., 2021). In taller shrub (chaparral, NFFL model 4), ¨ Okten and Liu (2021) also recognized shrub height as the most important variable for R, while Liu et al. (2015) and Ervilha et al. (2017) identified height as the most important structural variable influencing R. This contrasts with Clark et al. (2008), who found that fuel load was the most important variable, while Wadhwani et al. (2021), who used a neural network approach to predict R, confirmed the pre-eminent relative position of fuel depth in shrub communities of mixed species with markedly different structural characteristics. On the other hand, the uncertainty of W Shr_G1_dead and fuel depth can exert the greatest influence on the uncertainty of R, even in areas with mixed fuels, including shrublands (Cai et al., 2019). Regarding the sensitivity of other fire spread rate models to fuel structural variables, in the model of Anderson et al. (2015), a 10% error in fuel height and bulk density input variables resulted in a smaller (i.e. <10%) or proportional (i.e., ~10%) variation in the predicted rate of spread, respectively. Chatelon et al. (2022) indicated that the Balbi fire spread model showed a degree of sensitivity to vegetation height that depended on wind velocity. For instance, for low wind speeds (<2 m/s at shrub height), a doubling in shrub height can result in a 50% increase in R, while the same doubling in fuel height at high wind speeds could result in a 25–35% increase in R. Overall, the above observations highlight the importance of the accuracy in the measurements of fuel structural variables to improve the Fig. 4. Box plots of mean shrub height and dead fine shrub load (W Shr_G1_dead ) for each customized fuel model proposed for shrubland communities. Different letters indicate significant differences between mean values. The blue lines correspond to the values of the standard NFFL fuel models (FM4, FM5 and FM7) of Rothermel (1972) and Anderson (1982) and the red lines to the humid climate shrub models (SH3, SH4, SH6 SH8 and SH9) of Scott and Burgan (2005). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Fig. 5. Dichotomous classification criteria for customized fuel models based exclusively on mean shrub height. J.A. Vega et al. Journal of Environmental Management 351 (2024) 119831 9 performance of fire behaviour prediction models. In addition, the uncertainty in R can greatly affect the uncertainty of other fire behaviour variables to which it is closely related and thus affect fire management decisions. Most shrubland fuel classification systems use stand height as the first structural classifying variable, either as a mean value or as a range (e.g. Anderson, 1982; Scott and Burgan, 2005; Prometheus, 2000; Fernandes et al., 2009b; Cruz et al., 2018; S´ a et al., 2023; Solares-Canal et al., 2023). The limits of the height ranges in our models are greater than in more typically Mediterranean communities (e.g. Ria˜ no et al., 2002; Dimitrakopoulos, 2002; Rodríguez y Silva and Molina-Martínez, 2012), which seems consistent with a comparatively higher shrub development in our region. Shrub cover, a surrogate for fuel continuity, is also frequently used jointly together with height in shrub fuel classification, either as a rather qualitative descriptor (e.g. Scott and Burgan, 2005) or in quantitative ranges (e.g. FCCS, Prichard et al., 2013; Bushfire Fuel Classification of Australia, Kenny et al., 2019). Different thresholds of percentage shrub cover have been proposed in two of the most commonly used fuel classification systems: >50% in FBFM (Scott and Burgan, 2005) and 0–40%, 40–70% and >70% in FCCS (Prichard et al., 2013) while four ranges are considered in Australian Bushfire: <10%, 10–30%, 30–70% and >70% (Cruz et al., 2018). Different shrub cover thresholds have been proposed for Mediterranean shrublands: >60% (Ria˜ no et al., 2002), >40% (Jakubowksi et al., 2013) and >50% (Viedma et al., 2020). By contrast, Marino et al. (2016) considered a threshold of 50% for a Macaraonesian ecosystem. Further research on this topic is necessary, as indicated by the widely differing criteria applied to the parameter quantifying fuel horizontal continuity, which is a critical property for fire spread, and the lack of any reasoning justifying the ranges used. In the present case, the shrublands typically form closed stands and the combined effect of height and cover on R was accounted for through shrub-weighted height. 3.2. Customized fire behaviour models for bracken-dominated communities As already mentioned, four different situations were considered, as the morphological and physiological characteristics of brackendominated communities vary throughout the year: 1) conditions of full vegetative development (similar to those at the time of the inventory); 2) autumn conditions during which bracken is dead but standing at the same height as when the samples were taken (total load is considered dead fuel); 3) winter conditions during which all bracken is dead and moderately compacted (total load is considered dead fuel and mean height is adjusted to 40 cm based on field observations) and 4) full vegetative development conditions but assuming that only dead fuel burns, due to high moisture contents of live fine fuels, with values generally above 150%, and considering a mean height for this fuel of 40 cm based on field observations. Practical knowledge about bracken-dominated communities and associated fire behaviour in the study area initially indicated that two customized fuel models should be developed. This initial idea was ratified by the results of the determination of the optimal number of clusters, which showed the same result regardless of the morphological and structural state of the bracken communities. According to the criterion of maximization of the mean value of the silhouette variable, the optimum number of clusters is two, while the value of the sum of squares within the cluster decreases considerably for two clusters although it continues with a smoother decreasing trend even up to eight clusters. Fig. 6 shows the pattern of these two selection criteria considering full vegetative development (first condition) of the bracken-dominated communities analyzed in this study. The two clusters obtained for each of the four morphological and structural conditions of development of these communities were then compared. The mean distance from one sample plot to the other sample within the same cluster was never greater than the mean distance from the plots in the other cluster, i.e. the silhouette values for all sample plots and development conditions were positive. A total of 87.76% of the sample plots were classified similarly, regardless of the development condition considered. It is important to note that discrepancies only occurred when considering the second development condition (autumn, when bracken is dry, and the total load is therefore considered dead fuel, and standing at the same height as when the samples were taken), while the classifications and the medoids were the same for the other three conditions. These discrepancies probably occur because the second condition results in the highest (by far) rates of fire spread due to the large dead fuel load with low moisture and low compaction while standing. In any case, the medoids obtained for both clusters in this second development condition are similar to those obtained in the classifications of the other three developmental conditions; moreover, the results of the discriminant analysis indicated similar positions in the canonical space between the two groups of medoids. The results shown below were obtained in the cluster analysis of the sample plots in which the first morphological and structural condition of the bracken-dominated communities (full vegetative development) was considered. The discriminant analysis indicated that only one of 49 sample plots (2.04%) was misclassified and the structural variables with the highest weight in the canonical variable formulation were W Shr_G1 and W Shr , while R was the fire behaviour variable with the highest weight in its formulation. The first customized fuel model for bracken-dominated communities (Bracken-1; Fig S-7) is that with the lowest heights and fuel loads and includes 27 sample plots (55.1%). The second customized fuel model (Bracken-2; Fig S-8) includes 22 sample plots (44.9%) and corresponds to communities with higher mean shrub heights and fuel loads. Table 6 shows the mean value and standard deviation of the main structural variables of the bracken-dominated community in the sample plots included in each cluster as well as the values that correspond to the medoid of each cluster which represent the custom fuel model proposed. As shown in Table 6, any of the structural variables can be used to assign the custom fuel model to a new plot as there are significant differences between means for all of them. Fig. 7 shows, for example, the box plots with the distribution of hShr and W Shr for both customized fuel models. Finally, to facilitate assignment of one of the customized fuel models to the bracken-dominated formations, both in the field inventory and with remote sensing approaches, a classification was made based only on hShr, establishing a threshold of 1 m to differentiate between the Bracken-1 and Bracken-2 models (see Fig. 7). Using this threshold, 83.67% of the 49 sampling plots were correctly classified, 14.29% were Fig. 6. Graphs for determining the optimal number of clusters for brackendominated communities by maximizing the mean silhouette value and minimizing the within cluster sum of square of distances considering full vegetative development of the sample plots (first condition). J.A. Vega et al.