Estimating Fuel Loads and Structural Characteristics of Shrub Communities by Using Terrestrial Laser Scanning

: Forest fuel loads and structural characteristics strongly affect fire behavior, regulating the rate of spread, fireline intensity, and flame length. Accurate fuel characterization, including disaggregation of the fuel load by size classes, is therefore essential to obtain reliable predictions from fire behavior simulators and to support decision-making in fuel management and fire hazard prediction. A total of 55 sample plots of four of the main non-tree covered shrub communities in NW Spain were non-destructively sampled to estimate litter depth and shrub cover and height for species. Fuel loads were estimated from species-specific equations. Moreover, a single terrestrial laser scanning (TLS) scan was collected in each sample plot and features related to the vertical and horizontal distribution of the cloud points were calculated. Two alternative approaches for estimating size-disaggregated fuel loads and live/dead fractions from TLS data were compared: (i) a two-steps indirect estimation approach (IE) based on fitting three equations to estimate shrub height and cover and litter depth from TLS data and then use those estimates as inputs of the existing species-specific fuel load equations by size fractions based on these three variables; and (ii) a direct estimation approach (DE), consisting of fitting seven equations, one for each fuel fraction, to relate the fuel load estimates to TLS data. Overall, the direct approach produced more balanced goodness-of-fit statistics for the seven fractions considered jointly, suggesting that it performed better than the indirect approach, with equations explaining more than 80% of the observed variability for all species and fractions, except the litter loads.


Introduction
Fuel load, defined as the fuel dry biomass per unit area, is a particularly important forest fuel characteristic [1,2] as it represents the amount of live and dead biomass available for ignition and combustion. Along with other structural characteristics, fuel load strongly affects fire behavior, regulating the rate of spread, fireline intensity, and flame length [3,4]. Accurate fuel characterization, including the fuel load disaggregated by size classes and vegetative status (live and dead fractions), is essential to obtain reliable predictions from fire behavior simulators in order to support decision-making in fuel management and fire hazard prediction [5,6]. Fuel load is also an essential component in carbon cycle assessment and in models used to predict fire emissions and smoke impacts [7][8][9][10].
In the NW Iberian Peninsula, the fast shrub growth associated with mild annual temperatures and high rainfall promotes large accumulations of surface fuel that favor high intensity surface wildfires [11]. Indeed, the area is one of the most severely affected by fire in Europe [12,13]. Recent studies have reported an increase in fire activity in terms of both frequency and burnt area, mainly corresponding to eastern flows during summer, with a substantial contribution from southeastern advections in winter [14][15][16]. In Galicia, located in the westernmost part of this region, the problem is especially important in shrublands, as 63% of the total area burned in this region in the last 20 years has affected non-tree shrub communities [17], even though these communities represent only about 30% of the total forest area and 20% of the total area of the region [18].
The fuel load of shrub communities, i.e., shrub and litter layers, can be quantified by different methods, including destructive and non-destructive sampling [19]. Destructive sampling is the most accurate method of determining fuel load, but also the most complex and costly in time and money, and it is therefore difficult to apply to large areas [2,20]. Non-destructive sampling enables indirect estimation of the fuel load, usually by using allometric relationships between fuel loads and fuel characteristics that are easily measured in the field [21][22][23].
In the last two decades new models relating fuelbed characteristics to remotely sensed data have been developed by using information obtained from satellite images [24][25][26][27][28] or airborne laser scanner (ALS) data [29][30][31][32][33]. This has enabled the construction of thematic fuel maps for fire simulation and fire risk assessment. However, the accuracy of the model estimates largely depends on the accuracy of the digital terrain model (DTM) generated, which is, in turn, affected by the density of the vegetation and the orography [34]. Thus, in very dense shrubland and rugged terrain, typical of the Galician territory, shrub height and volume tend to be underestimated [35][36][37].
Terrestrial laser scanning (TLS) has evolved as a non-destructive measurement technique that produces data (3D point clouds) of great potential for use in forestry applications. Such uses include the assessment of tree and stand level biomass and carbon [38][39][40][41] and also the accurate quantification of relevant fuelbed characteristics for determining fire behavior, fuel consumed and smoke emissions [42][43][44][45]. In addition, studies in which TLS data are combined or compared with ALS data [46][47][48][49][50], data from unmanned aerial vehicles [51,52] or satellite images [53] have become more frequent in recent years for various forestry applications.
Most TLS-based studies estimating shrub biomass or volume have focused on the reconstruction of the shrub shape from the point cloud, either by piecewise approximation of the volume or by voxelization. For example, Olsoy et al. [43] estimated the biomass of sagebrush (Artemisia tridentata) using both voxel counting and a 3-D convex hull approach, concluding that the 3-D convex hull model estimated total and green biomass more accurately than the voxelization approach. Greaves et al. [54] also compared both approaches for estimating harvested biomass of two dominant Alaskan tundra shrub species, using point clouds obtained from close-range (2 m) and variable-range (<50 m) TLS, and reported better results for voxel counting in close-range data and stronger relationships for the surface difference approach for variable-range data. Hudak et al. [49] and Rowell et al. [55] adapted the "3D field sampling protocol" proposed by Hawley et al. [56] to infer fine-scale fuel mass in 3D, observing that mass was strongly related to the voxel-occupancy volume, porosity, and surface area, which were used to develop equations for fine-scale prediction.
In general, fuel characterization studies with TLS are carried out in tree-covered shrub communities, often characterized by a heterogeneous spatial distribution on the ground, and they usually focus on obtaining detailed fuel information for applying physical-based fire models, such as WRF-SFIRE-CHEM [57,58], WFDSS [59,60], and FIRETEC [61,62], which require spatially explicit, three-dimensional data [63]. However, the use of TLS data should not be restricted to obtaining geometric information on a scale as detailed as that required for the use of physical-based models that operate in 3D, but should be tested for potential use in obtaining the fuel characteristics required for the application of more operational simulators such as Farsite and Flammap [64,65]. In addition, such detailed studies require several scans from different positions, often at least 5 to 9 scans, indicating a very time-consuming field inventory and subsequent processing of the point clouds to obtain a single cloud. Studies that use TLS to predict biomass distribution without reconstructing shapes from the returns are scarce. Accordingly, the aim of this study was to test the potential use of TLS data to provide accurate fuel load estimates of size-disaggregated shrub and litter and live/dead fractions for four of the main non-tree covered shrub communities in Galicia (NW Spain) by using data from a single scan. For this purpose, two alternative approaches will be compared, one direct estimation approach and a two-steps indirect estimation approach.

Shrub Communities and Plots
Shrubland areas in the study area (Galicia, NW Spain) are usually covered by gorse (genus Ulex), heath (genus Erica) and broom (genera Pterospartum, Cytisus and Genista), covering 55%, 22%, and 11% of the shrubland area, respectively [66]. Four treeless shrub communities were selected and characterized by the dominant species (Table 1). The communities selected are widely represented in Galicia [18,66] and are among the most frequently affected by wildfires in the region. All the communities considered are multiannual woody shrubs with stems branching from the base. Subjective sampling was carried out to select a representative sample of each of the 4 shrub communities by taking the main structural fuel attributes observed in the study area into account. The plots were located using information obtained from different sources, particularly the Spanish Forest Map [66,67]. A total of 55 sample sites were selected in areas where the dominant shrub species of each community predominated, and shrub cover was higher than 80%. The spatial distribution of the sampled plots throughout the study area is shown in Figure 1.

Field Inventories and Determination of Fuel Characteristics
Field measurements were made in the shrub and litter layers of each sample plot by randomly locating a 2 × 2 m sampling square with sides parallel to the contour line and to the maximum slope. Each sampling square was delimited with four wooden stakes and all the vegetation in a 0.5-m strip around the square was removed, except the side where the TLS was located, where a 4-m strip of vegetation was removed.
A sampling transect consisting of the plot perimeter plus one diagonal of the square was outlined. The horizontal lengths of the intercepted shrub species (cm) were measured along this linear transect, to estimate the linear shrub cover (CovShr) according to the method proposed by Canfield [68]. Moreover, dead and live base height, total shrub height of the intercepted species (cm) and litter depth were also measured (every 50 cm). As the separation between the litter and duff layers is not always very clear in shrub communities, both layers were sampled together as litter. Shrub height was measured as the vertical distance between the litter layer and the top of vegetation. All the shrub and litter structural measurements were averaged to produce the mean shrub height (ℎ ) and mean litter depth ( ). For each community, the fuel loads were estimated using the species-specific biomass equations systems at stand level developed for treeless shrub communities in Galicia [69], with CovShr, ℎ and as inputs. These systems include equations for each of the following seven fractions: dead fine shrub load (diameter < 0.6 cm; WShr_G1_dead); live fine shrub load (diameter < 0.6 cm; WShr_G1_live); fine shrub load (diameter < 0.6 cm; WShr_G1, dead + live); coarse shrub load (diameter ≥ 0.6 cm; WShr_G23); total shrub load (WShr, fine + coarse); litter load (WLitt); and total shrub; and litter fuel load (WShr+Litt). The basic descriptive statistics of shrub and litter layers for the main structural characteristics of each shrub community are shown in Table 2. Table 2. Description of the main shrub and litter characteristics. Std. dev. is the standard deviation, n is the number of plots, ℎ is the mean shrub height, CovShr is the shrub cover, is the mean litter depth, WShr_G1_dead is the dead fine shrub fuel load, WShr_G1_live is the live fine shrub fuel load, WShr_G1 is the fine shrub fuel load, WShr_G23 is the coarse shrub fuel load, WShr is the total shrub fuel load, WLitt is the litter fuel load, and WShr+Litt is the shrub and litter fuel load.

TLS Data
Each sampling square was scanned once using a FARO Laser Scanner Focus3D X 130 TLS to collect three-dimensional point clouds with a density of ~7.67 mm point spacing at 10 m. The TLS is a class I laser (wavelength = 1550 nm) with a range of 0.6 to 130 m for objects with 90% reflectivity and with an accuracy of ±2 mm at 10 m range. Scanning was conducted by locating the TLS at a distance of 3 m from the contour line side of the square sampling located downhill.
The point clouds from the single scans on each sampling square were cropped and height was normalized. All points outside of the squares were removed, and the height above the terrain of each point (h) was calculated from digital elevation models (DEMs) created from each scan. DEMs were generated by interpolation, according to the lowest laser return within a 2 × 2 cm grid, in order to maintain consistency with the field measurements. The heights were then calculated by subtracting the smoothed DEM from the z values. These values were used to calculate 32 metrics related to the vertical distribution of the point clouds. The metrics and the corresponding descriptions are summarized in Table 3. Table 3. Potential independent variables related to height distribution.

Metric
Description interquartile distance h01, h05, h10, h15, h20,…, h90, h95, h99 (cm) percentiles h25 and h75 (cm) first and third quartiles IQR (cm) Interquartile range (h75-h25) PAhmean ratio of laser returns above hmean to total laser returns PAhmode ratio of laser returns above hmode to total laser returns Due to shrub occlusions, the TLS-derived metrics vary as the depth of the analyzed laser returns cloud measured from the shrub front increases. Therefore, analysis was carried out to determine the optimum depth (measured from the vertical plane defined by the lowest side of the square) from which the metric values are stabilized, by comparing the values for depths varying from 10 to 50 cm, by 1 cm. This procedure is equivalent to constructing a virtual vertical plane parallel to the sampling square front and extracting all the points that are closer to the scanner than that plane. The results indicated that 30 cm depth of the laser returns cloud was adequate for all of the sampled plots, so the laser returns clouds were limited to this depth threshold ( Figure 2).
The theoretical vertical shrub cover was estimated as the ratio between the observed number of laser returns and the maximum number of laser returns lrmax_ver. The maximum number was estimated considering the vertical plane defined by the lowest side of the square as a continuous solid with height equal to hmax and taking the laser scanner density into account. Moreover, the 3-D point clouds for each sample plot (2 × 2 × hmax) were divided into voxels (1 cm 3 ), and the volume of shrub cover was calculated as the ratio between the number of voxels containing points and the total number of voxels. Both of these ratios were transformed by calculating the arcsine-square root of their values to stabilize the variance and improve normality [70], thus producing the Bliss variables TLSvCovShr and TLS3DCovShr. Finally, the vertical point distribution was modelled using the two-parameter Weibull probability density function (pdf). For this purpose, the relative frequency of the number of laser returns in each layer of height 1 cm and depth 30 cm was calculated by dividing the number of returns in this layer by the total number of laser returns. The scale (Wscale) and shape (Wshape) parameters of the Weibull pdf were then estimated using the first and the second moments of the relative distribution of vertical laser returns [71].
Point cloud processing, calculation of the metrics derived from the TLS and estimation of the parameters of the Weibull pdf were programmed in R [72] with the support of the rgl [73] and e1071 [74] libraries.
The Bliss transformed ratios (TLSvCovShr and TLS3DCovShr), the Weibull pdf parameters (Wscale and Wshape) and the 32 metrics related to the vertical distribution of the point clouds (Table 3) were finally considered as potential independent variables in the models for estimating fuel loads.

Models for Estimating Fuel Loads
Two different approaches for estimating fuel loads from community and biomass fractions were compared: indirect estimation (IE) and direct estimation (DE). The IE approach is based on fitting a set of three equations relating the inputs of the biomass equations developed for these shrub communities in Galicia [69] from usual field measurements (CovShr, ℎ and ) to TLS-based metrics and variables and then estimating the fuel loads by using these equations. The DE approach is based on fitting a system of seven equations relating the fuel load of each biomass fraction to TLS-based metrics and variables.
Linear models ( = + ) were used to fit ℎ and equations in the IE approach and a logistic model with a horizontal asymptote equal to 100 ( = 100 (1 + exp ( ⁄ + )) was used for CovShr. The best set of independent variables (Xi) was selected using the stepwise variable selection method by previously log-linearization of the logistic model. Allometric models ( = · ) were tested for the seven biomass fractions fitted in the DE approach. These allometric equations must fulfil the property of additivity, i.e., the sum of biomass predictions from separate fuel fractions must equal the biomass prediction from the total biomass model (e.g., the sum of dead fine and live fine shrub load estimates must equal fine shrub load estimates or the sum of fine and coarse shrub load estimates have to equal total shrub load estimates). A disaggregation method focused first on specifying an equation for total biomass and then constrained the form or parameters of fraction equations to ensure additivity was used [75,76].
In the first step, the equation of each fraction was log-linearized and fitted separately using ordinary least squares (OLS) for linear models. The stepwise approach was used to select the best set of independent variables (Xi), and the complete system of seven equations (one per fraction) was then fitted simultaneously to guarantee additivity. The final system of seven equations to be fitted consisted of the following: 1. An equation for estimating total shrub and litter fuel load 2. Two allometric equations to estimate total shrub load ( ) and litter load ( ): = + ( ) and = + ( ) ; therefore, considering the ratio between and and solving by Equation (1) was disaggregated in two equations, the first one to estimate the total shrub fuel load: where b0 = b0Litt − b0Shr; bi = biLitt − biShr and the second one to estimate the litter fuel load: 3. Two allometric equations to estimate coarse shrub load ( Equation (2) is disaggregated in two equations, one for estimating the coarse dead fuel loads: where c0 = c0g1 − c0g23; ci = cig1 − cig23 and a second one for estimating fine shrub fuel load: 4. Finally, disaggregation of Equation (5)  The equation for estimating the fine dead fuel loads was then as follows: where d0 = d0g1_live − d0g1_dead; di = dig1_live − dig1_dead and the equation for estimating the live fine fuel load was as follows: The condition number was used to evaluate the presence of multicollinearity among variables in the fitted equations. According to Myers [77], condition numbers higher than √1000 indicate problems associated with multicollinearity.
The presence of heteroscedasticity was checked using the White test [78] and by visual inspection of studentized residuals plotted against fitted values. When heteroscedasticity was detected, each observation was weighted by the inverse of its estimated variance ( ), under the assumption that this variance can be modelled as a power function of the independent variables [79], i.e., = ( ) . The value of the exponential term k was optimized to provide the most homogeneous studentized residual plot by using the method proposed by Harvey [80]. The systems were fitted using the nonlinear seemingly unrelated regression (NLSUR) method, which takes the cross-equation correlations into account. The cross-equation error covariance matrix obtained by ordinary least squares was used to initiate the iterative procedure. When necessary, the weighting factor for heteroscedasticity was programmed in the MODEL procedure of SAS/ETS  [81].
Two goodness-of-fit statistics were used to check the accuracy of estimates: the coefficient of determination (R 2 ) and the root mean square error (RMSE).
where Yi, i Yˆ, and are the observed, predicted and mean values of the dependent variable, and n is the total number of observations used to fit the equation.
As the number of sample plots per shrub community is low for fitting the systems for each, the complete database was used to compare the IE and DE approaches. However, differences between shrub communities were evaluated by expanding the parameters a0, b0, c0, and d0 with a linear relationship including associated parameters with dummy variables to differentiate the four communities: where α0 is the expanded parameter (a0, b0, c0 or d0); α00 is a common parameter for all communities; α0i (i = 1, …, 3) is an associated parameter, and Ij (j = Ea, Ug, Pt) is a dummy variable whose value is equal to 1 for the specific community and 0 for the other shrub communities. If the result of the t-test of one specific associated parameter is significant (5%), this indicates differences between communities on the specific equation affected by the associated parameter. Finally, the results of DE and IE approaches will be compared by calculating the R 2 , the RMSE, the bias (mean error ̅ ), the absolute bias (mean error in absolute value | | ), and the percent RMSE (the percentage of RMSE over the mean value of the dependent variable RMSE%) from the observed and estimated fuel loads in the sample plots for each approach.

The Direct Estimation (DE) Approach
The final expressions of the system of equations fitted (Equations (1)- (7)) to the complete database of four shrub communities, including the best set of independent variables, the values of the parameter estimates, the asymptotic standard errors and the goodness-of-fit statistics, are shown in Table 4. The results of the White test and the values of the condition number did not indicate any problems of heteroscedasticity or multicollinearity for any fuel fraction.
The values and signs of all the parameters were biologically consistent and the graphical inspection of studentized residuals showed random patterns of residuals around zero, with homogeneous variance and no discernible trends. Plots of observed versus predicted values are shown in Figure 3.
All fitted equations, with the exception of the equation estimating WLitt, explained more than 80% of the observed variability, with values ranging from 80.33% to 93.92%, even exceeding 90% in the equations estimating WShr (90.18%) and WShr_G1 (93.92%). Differences between shrub communities were observed for most of the fuel fractions, except WLitt and WShr. These results are similar to those obtained by other authors. Thus, for example, the linear equations proposed by Ku et al. [82] for estimating Prosopis glandulosa biomass in rangelands in the southwestern USA yielded R 2 values of 0.77-0.81 for respectively TLS-derived height bins and percentile heights as predictors; the models proposed by Olsoy et al. [43] for estimating total biomass and green biomass of individual plants of Artemisia tridentata in the sagebrush-steppe ecosystem in the western USA explained respectively 92 and 83% of the observed variability; for the same species and study area as in the aforementioned study, Li et al. [47] proposed a model that explained 87% of the observed variability; Greaves et al. [32] fitted models to estimate aboveground total biomass of low-stature (<1.5 m) Arctic shrub species in the northern Alaskan tundra, obtaining R 2 values ranging from 0.91 to 0.94 depending on the scanning conditions (2 m close-range or <50 m variable-range); the linear models developed by Owers et al. [45] for estimating biomass in coastal wetland vegetation in southeast Australia explained between 42 and 99% of the observed variability, depending on the volumetric model used.
Although there are slight differences in the fitted equations, depending on the shrub community considered, the proposed system can be considered general. Because of structural and functional differences in the four communities analysed [11], improvement of the goodness-of-fit statistics is expected when a larger sample of plots to fit community-specific systems is available, especially for coarse shrub fuel load (WShr_G23), in which the inter-species variability is greater.
As previously explained, the shrub fuel loads of each fraction and community were not destructively sampled, but estimated from previously fitted equations, with CovShr, ℎ and as independent variables. Therefore, as expected, the independent variables required as inputs to initiate the disaggregation system fitted ( equation) are TLS-based metrics related to shrub height and cover (h95 and TLSvCovShr). Three other metrics are involved in the disaggregation process: one related to the shrub height (h70) and another two related to the vertical distribution of the fuel (Wscale and IQR). The latter two variables are the only ones that discriminate between coarse and fine fuel (WShr_G23 and WShr_G1) and the categories of live and dead fine fuels (WShr_G1_live and WShr_G1_dead), highlighting the importance of the vertical distribution of fuels in the biomass loads of these fractions. Rowell et al. [44] already highlighted the importance of vertical distribution of fuels in their total load and the potential of metrics derived from the vertical height distribution of the TLS laser returns cloud for that purpose. Furthermore, the equations fitted for these fractions were also those most affected by structural differences between the four shrub communities analysed, resulting in significant community specific-parameters. Table 4. Final expression of the system of equations fitted simultaneously to estimate fuel loads of each fraction, including the best set of independent variables, the parameter estimates, the approximate standard errors and the goodness-of-fit statistics. WShr+Litt = shrub and litter fuel load, WShr = total shrub fuel load, WLitt = litter fuel load, WShr_G23 = coarse shrub fuel load, WShr_G1 = fine shrub fuel load, WShr_G1_dead = dead fine shrub fuel load, and WShr_G1_live = live fine shrub fuel load.

Equation Expression
Par

The Indirect Estimation (IE) Approach
The values of the parameter estimates, the asymptotic standard errors and the goodness-of-fit statistics for the three equations relating CovShr, ℎ and to TLS-based metrics and variables are shown in Table 5. No problems of multicollinearity were detected for any of the three equations fitted. As expected, the best results of the goodness-of-fit statistics were obtained for ℎ , with a model explaining more than 91% of the observed variability; 66.26% for CovShr and the poorest results for with a model explaining a 39.07% of the observed variability. Similar results were obtained in studies characterizing shrub and using TLS. Thus, e.g., Vierling et al. [37] reported values of explained variability in height and shrub cover of respectively 94 and 51%, in Artemisia tridentata dominated ecosystems in eastern Washington, USA. Anderson et al. [83] fitted models to estimate the cover fraction of shrubs (R 2 = 0.77), annual grasses (R 2 = 0.70), perennial grasses (R 2 = 0.36), and forbs (R 2 = 0.52) using the random forest algorithm in desert grasslands and big sagebrush shrublands in southwest Idaho, USA. In a study in the agro-pastoral ecotone of northern China, Xu et al. [84] proposed a model for estimating grass height and canopy cover from TLS data, which explained 70 and 72% of the observed variability respectively. Better results were obtained by Chen et al. [85], who, in a study in South-East Australia, established models for estimating the litter depth and percentage cover of the different vertical layers of the understory from TLS metrics, with R 2 values of 0.9 and of 0.8-0.9, respectively.
The low accuracy of the estimates of the litter layer depth equation is consistent with the low accuracy of the estimates of the litter fuel load equation fitted in the DE approach (R 2 = 0.4080). These results were expected due to the strong relationship between the litter layer depth and the accumulated load per unit area and the difficulty in characterizing this depth from the TLS point clouds because of limited visibility in high-density shrub communities [49]. Therefore, because of the importance of this fraction in both the flame and the smouldering phases, in situations in which this fraction could contribute substantially to energy flux and emissions from surface fires, destructive sampling is required to complement TLS scans.
The estimates of CovShr, h and d obtained with the three models were used as inputs of the biomass equations developed for the studied communities in Galicia [69] to estimate the fuel load of WShr+Litt, WShr, WLitt, WShr_G23, WShr_G1, WShr_G1_dead, and WShr_G1_live. Plots of observed versus predicted values of these seven fuel loads are shown in Figure 3 for comparison with the same plots obtained using the DE approach.  Table 6 compares the values of R 2 , RMSE, bias (mean error ̅ ), absolute bias (mean error in absolute value | | ), and the percent RMSE (the percentage of RMSE over the mean value of the dependent variable RMSE%) from the observed and estimated fuel loads in the sample plots for each approach. Table 6. Goodness-of-fit statistics obtained with the direct estimation (DE) and indirect estimation (IE) approaches for each fuel fraction.

Variable Direct Estimation (DE) Approach
Indirect Estimation (IE) Approach Overall, the best results, in terms of goodness-of-fit statistics, were obtained with the DE approach, although for WLitt and WShr_G23 the IE approach explained more of the observed variability, with a reduction in RMSE of, respectively, 2.6 and 13.8%, relative to the DE approach. However, the plots of observed versus predicted values showed a slight bias in all the estimates in the IE approach, with a slight tendency to underestimate the values, whereas this trend was not appreciable in the plots of the estimates of the DE approach (Figure 3).
The system of compatible equations developed with the DE approach is proposed for estimating the fuel loads of the shrub communities analyzed, with the equations fitted for the IE approach used to estimate three other important variables for complete characterization of fuelbeds (CovShr, ℎ and ). Apart from the important improvement regarding the provision of a compatible system of equations that guarantees additivity, the proposed methodology does not add any apparent advantages in field sampling relative to the equations based on variables obtained from non-destructive sampling, such as those already existing for the species analyzed in the study area [69]. However, as Newnham et al. [86] pointed out, TLS should not be considered just as an extension or automation of traditional field data sampling, but rather as a completely different methodology to characterize the spatial structure at the plot scale with an unprecedented level of detail. TLS point clouds can easily be processed with simple algorithms, such as those used in this study, to provide useful additional information for characterizing fuelbeds that would otherwise be difficult, impractical, or impossible to obtain with non-destructive field sampling. For example, TLS can characterize the vertical distribution of fuel loads with greater precision and at a finer scale than afforded by traditional methods (see Figure 4). Moreover, the differences in water content between green and woody biomass produce differences in reflectivity of the laser energy [87,88] enabling the TLS point cloud to be divided into green and non-green fractions [43]. In addition, previous studies have shown that TLS allows for very accurate estimates of other fundamental variables describing structural characteristics of fuelbeds, such as bulk density [49,89] or packing ratio [55]; both of which strongly influences the flammability and combustibility of a given fuelbed. For these reasons, TLS could provide improved estimates of fuelbed properties driving fire behavior, thereby enhancing the accuracy of fire behavior prediction simulators [64,65].  Several authors have highlighted the usefulness that could have the use of estimation models derived from TLS data as a stepping stone to train remote sensor databases applicable at large scale [47,48,83,90]. Therefore, the development of accurate estimation methodologies, such as those presented in this paper, based on a single TLS scan, with the consequent reduction in costs and time, constitute the basis for the application of TLS in the forest inventory, combining it with inference methods such as model-based or model-assisted approaches.
Finally, the equations proposed in this study could be useful for characterizing shrub communities from the point of view of fire behavior and fire hazard assessment, as well as for designing prescribed burns, for simulation of erosion processes, for estimating carbon stocks or biomass potential as energy source and for characterizing and delimiting wildlife habitats. Although the approach proposed is transferrable to point clouds derived from other shrub communities, new models should be fitted based on destructive sampling taking into account the specific ecosystem conditions, and considerations of timing due to phenology.

Conclusions
Forest managers require accurate models for estimating shrub fuel loads, especially in fire-prone areas such as Galicia, where about 2% of the forest land area is burned annually (period 1978-2019) and some two-thirds of which correspond to shrub vegetation. In this study, two alternative approaches were compared for estimating size-disaggregated shrub and litter fuel loads and live/dead fractions from TLS derived metrics, for four of the main treeless shrub communities in Galicia. The results obtained suggest that the use of metrics derived from a single TLS scan as independent variables produces accurate estimates of shrub fuel loads, mainly by using the proposed systems of equations of the DE approach.
TLS research on surface fuels is limited, and previous studies have focused on the application of physical-based fire models requiring spatially explicit three-dimensional data [49,55,63]. In addition, most of these studies have analyzed the total biomass, at best distinguishing between green and non-green fuels [43]. The proposed method therefore represents an important improvement in relation to characterization of shrub community fuel loads.
In northwest Spain, the abandonment of agricultural land, among other factors, has generated a large increase in the wildland-urban interface, creating complex mosaics in which the fire-prone shrub communities analyzed in this study predominate. Efficient management of these ecosystems increasingly demands more accurate and detailed fuel maps for comparing and prioritizing different fuel reduction strategies and ecosystem restoration. This task will require assembling and combining databases of estimates obtained from both traditional sampling and the use of remote sensors. The study findings suggest that TLS may be valuable for achieving this aim.
Our purpose is to advance towards developing a methodology that ultimately enables characterization of the spatial and vertical distribution of shrub fuels, distinguishing between size and live/dead fractions. The present research is an initial exploratory study of potentially suitable approaches to achieving these objectives. Further sampling that covers all combinations of development stages and site qualities of each shrub community, also including destructive sampling, is required to establish TLS as a robust technology for estimating the main characteristics related to fire behavior and for providing the basis for scaling up from plot to landscape.