Canopy and surface fuels measurement using terrestrial lidar single-scan approach in the Mogollon Highlands of Arizona

Study area

This study was conducted along an elevational gradient near Globe, AZ, in the Mogollon Highlands Ecoregion (Fig. 1). The climate is semi-arid with bimodal precipitation patterns with mean annual precipitation ranging from 400 mm at low elevation (1200 m) to 650 mm at high elevation (2000 m) and mean annual temperature ranging from 16.4 to 12.4°C (Johnson and Pase 1968; PRISM Climate Group 2021).

Fig. 1.

Top: map of study area and sampling plot locations in east central Arizona, US. Bottom: examples of each vegetation type.

We measured three vegetation types (Fig. 1) that are dominant across parts of the Southwestern US designated as high fire risk and high treatment priority (US Forest Service 2022). Species presence is similar across all three types but there are important differences in species abundance and structure. (1) Madrean ponderosa pine–oak forest (PPO) often has a two-tiered structure with a distinct overstory of ponderosa pine (Pinus ponderosa Douglas ex C. Lawson) and a dense understory of evergreen oaks (Quercus arizonica Sarg. and Q. emoryi Torr.) and other woody species such as alligator juniper (Juniperus deppeana Steud.), pinyon pine (Pinus edulis Engelm.) and pointleaf manzanita (Arctostaphylos pungens Kunth). (2) Madrean pinyon–juniper–oak woodland (PJO) is characterized by low-stature trees and tall shrubs such as pinyon pine, Utah juniper (Juniperus osteosperma (Torr.) Little), evergreen oaks and alderleaf mountain mahogany (Cercocarpus montanus Raf.), and at lower elevations includes velvet mesquite (Prosopis velutina Wooton), catclaw acacia (Senegalia greggii A. Gray), desert ceanothus (Ceanothus pauciflorus Sessé and Moc. ex DC), Asparagaceae (Nolina microcarpa S. Watson, Dasylirion wheeleri S. Watson ex Rothr., Yucca baccata Torr., and Agave spp.) and Cactaceae (Opuntia spp. and Cylindropuntia spp.). Our plots covered a range of density and structure in PJO including savannah, persistent woodlands and wooded shrubland (Romme et al. 2009). (3) Mogollon chaparral (CHAP) is dominated by a dense shrub layer 0.5–2 m in height composed of pointleaf manzanita and shrub-form oaks sometimes mixed with sparse trees, typically pinyon pine, Utah juniper, or tree-form evergreen oaks. Across all vegetation types, grass and forb cover was sparse owing to woody vegetation and widespread livestock grazing. There was no field evidence or records of wildfires or fuel treatments within our sampling area in the recent decades preceding our field measurements (Southwest Ecological Restoration Institutes 2025).

Owing to accuracy issues with existing classification maps, we manually delineated patches of PPO and CHAP based on field visits and high-resolution imagery, then classified the remaining matrix as PJO. To ensure that we captured a wide range of PJO density, we delineated three density-based subtypes using normalized difference vegetation index (NDVI) derived from Sentinel-2 satellite images (European Space Agency, https://dataspace.copernicus.eu). NDVI and related indices are widely used in studies that map and classify vegetation density in semi-arid environments (Humagain et al. 2017). We delineated low-density PJO, medium-density PJO and high-density PJO using three approximately equal bins by thresholding on the 30th and 70th percentiles of NDVI across all PJO in the study area. We did not repeat this subtyping process on CHAP or PPO because there was not as much variability in density in those vegetation types within the study area. We visited 20 sampling plots within the chaparral type, 20 plots within the ponderosa pine–oak type and 27 plots within the pinyon–juniper–oak type (approximately nine plots in each NDVI subtype). We selected sampling locations randomly with constraints regarding distance from a road (50–600 m) and maximum slope (30°) due to logistical constraints. All field data were collected between June and August of 2023.

Lidar data collection and processing

We collected fuel measurements using the Leica BLK 360 Gen 1 lidar scanner (Leica Geosystems, Heerbrugg, Switzerland). This simple and lightweight unit uses time-of-flight technology to record single return measurements. The scanner has a horizontal field of view of 360° and a vertical field of view between 45° below the horizon and 90° above the horizon. All scans were performed with the ‘medium’ scanning density, which records approximately 5500 vertically oriented scan lines with 2700 points per scan line, yielding approximately 9600 points/m² on a perpendicular plane at a distance of 10 m.

We collected four scans at each plot to compare single-scan versus multiple-scan modes. For single-scan mode, we positioned the scanner directly over the plot center at a height between 1.6 and 2 m. For multiple-scan mode, we collected the center scan and three additional scans positioned strategically around the plot center to minimize overall occlusion (Fig. 2), similarly to Donager et al. (2021) and Batchelor et al. (2023). Most perimeter scans were 7–12 m from the plot center. We selected this number of scans and scanning configuration to balance low occlusion with low set-up time. We used reference targets and Cyclone Register software (Leica Geosystems 2023) to align and export the point clouds. Subsequent data processing was performed in Python (Python Software Foundation, https://www.python.org) using a package developed for this project that automates lidar processing, canopy bulk density prediction and most of the model training (Fig. 3). Lidar data from each plot were processed independently for single- and multiple-scan modes.

Fig. 2.

Diagram of sampling plot design. Sampling plots included an 11.3 m radius macroplot where all trees and shrubs more than 1 m above ground were measured, four 1 × 1 m subplots where live and dead surface fuels less than 1 m above ground were measured, and four terrestrial lidar scans. Position of the perimeter scans was adjusted based on vegetation arrangement in order minimize overall occlusion and ensure visibility of reference targets (not shown). In single-scan mode, only Scan 1 is considered. In multiple-scan mode, all four scans are considered. In both scanning modes, lidar-derived leaf area density data were spatially clipped to the 11.3 m radius macroplot.

Fig. 3.

Diagram of workflow to pre-process lidar data, then train and test canopy fuel prediction model.

Our primary lidar metric was LAD. Here, ‘leaf area density’ is synonymous with ‘plant area density’ because we did not thoroughly remove woody material (branches) from the analysis. To calculate LAD, we followed Grau et al. (2017) and traced the path of each laser pulse through a 3D grid with a 10 cm resolution (Fig. 4). For each 3D grid cell (i.e. voxel), we counted the number of pulses directed at the voxel, the number of pulses entering the voxel, the number of pulses intercepted within the voxel and the number of pulses transmitted through the voxel. We then used the Beer–Lambert Law to estimate LAD for each voxel based on the ratio between the number of pulses entering and the number of pulses transmitted. As recommended by Grau et al., we adjusted for unequal path length. We assumed a spherical leaf orientation function. We calculated occlusion for each voxel as one minus the number of pulses entering the voxel divided by the number of pulses directed at the voxel. To ensure that all pulses were considered in the path tracing counts, we used the .ptx export format from the Leica Cyclone software (Leica Geosystems 2023), which organizes returns along vertical scan lines and includes records for null-return pulses (e.g. pulses directed at empty sky). To estimate the trajectory of null-return pulses, we calculated spherical coordinates for all returns, then interpolated along each vertical scan line to estimate spherical coordinates for the null-return pulses. To maintain geometric integrity and ensure all pulses were considered in transmission estimates, we did not normalize the point cloud by ground elevation or filter out any returns prior to calculating LAD and occlusion. We classified voxels as non-foliage (opaque objects such as tree trunks and rocks) if the estimated LAD was greater than 6 m²/m³, a threshold we selected through visual analysis (Fig. 4). We classified voxels as occluded if occlusion was greater than 80%. We used a voxel size of 10 × 10 × 10 cm and an occlusion threshold of 80% based on previous studies (Béland et al. 2014; Rowell et al. 2020; Nguyen et al. 2022) and a series of optimization tests to maximize relationships between canopy fuel metrics, single-scan LAD and multiple-scan LAD. When aggregating multiple scans, we used the minimum occlusion value of each voxel, meaning voxels only needed to be visible from one scanning position.

Fig. 4.

Top: 3D isometric view of a single scan in PJO showing voxels with LAD estimates (yellow-green), occlusion mapping (translucent blue) and non-foliage elements (brown). Bottom: profile view of a 1 m wide strip from a single scan in PPO.

After estimating metrics for each voxel, we produced a digital terrain raster aligned with the voxel grid. To estimate terrain elevations, we used an iterative height filtering approach (Caster et al. 2021; Tango et al. 2024) to produce a 0.5 m raster, then downscaled to 10 cm resolution using linear interpolation. We used the terrain raster to estimate height above ground for each voxel. We clipped the voxels to a circular plot with an 11.3 m radius, then calculated the mean LAD within 10 cm bins of height above ground to produce a LAD profile. When calculating mean LAD, we treated occluded and non-foliage voxels as null. With this approach, empty voxels (LAD = 0) are included if they are visible, but occluded voxels and non-foliage elements are removed from the sampling space.

Conventional data collection and processing

We estimated canopy bulk density (CBD) along a vertical profile using conventional stem measurements and allometric equations within circular macroplots with an 11.3 m radius (Fig. 1). We estimated surface fuel load using four 1 × 1 × 1 m subplots within each macroplot (Fig. 2). The macroplot size was chosen because it can be related to 20 × 20 m Sentinel-2 satellite data and is nearly the same size as the 0.1-acre macroplots recommended by the FIREMON protocol (Lutes et al. 2006).

Using a laser hypsometer and tape measures, we measured stem sizes and height for all vegetation taller than 1 m within the macroplot. Following recommendations from Lutes et al. (2006), for all trees greater than 7.5 cm in diameter, we recorded species, diameter at breast height (DBH), total height, height to live crown and snag condition. For remaining shrubs and small trees taller than 1 m, we recorded occurrences by stem cluster. For each stem cluster, we recorded species, quadratic mean diameter of stem cluster, count of stems in cluster and representative cluster height. As an approximation of available canopy fuel, we estimated foliage biomass using allometric equations for ponderosa pine (Kaye et al. 2005), pinyon pine (Grier et al. 1992), juniper species (Grier et al. 1992), oak species (Clary and Tiedemann 1987), mesquite species (McClaran et al. 2013), manzanita species (Hughes et al. 1987) and ceanothus species (Hughes et al. 1987). These species accounted for more than 98% of the total basal area. We used the oak, mesquite, or ceanothus model for any remaining trees or shrubs over 1 m in height, depending on their foliage density.

After estimating available canopy fuel for each tree, we estimated CBD along a vertical profile (Sando and Wick 1972), similarly to tools such as the Fire and Fuels Extension to the Forest Vegetation Simulator (Reinhardt and Crookston 2003). We used height bins of 10 cm to match the lidar voxel size and smoothed the profile using a spline. From the CBD profile, we extracted the maximum canopy bulk density (CBD_max), canopy base height (using a minimum threshold of 0.011 kg/m³) and total canopy fuel load (sum of CBD profile).

Surface fuels were evaluated at each plot using four 1 × 1 × 1 m subplots located at 4 m north, 8 m south, 8 m east and 4 m west of the plot center (Fig. 2). In each subplot, fine woody debris loads were estimated using the photoload method (Keane and Dickinson 2007; Sikkink and Keane 2008). Litter and duff loads were estimated for each subplot by measuring litter cover and depth and multiplying by estimates of bulk density from a national reference dataset (Ottmar et al. 2007a, 2007b). Live surface fuel loads and dead standing herbaceous fuel loads (e.g. dead grass that is still rooted) were estimated by measuring vegetation cover by cover class and applying allometric models developed for ponderosa pine understory (Mitchell et al. 1987). For subsequent analysis, we allocated surface fuels into two classes: (1) standing surface fuels, comprising live graminoids, live forbs, live shrubs and dead standing herbaceous biomass; and (2) downed surface fuels, comprising litter, duff and fine woody debris (<7.6 cm in diameter). We calculated plot-level surface fuel metrics using the mean of the four subplots in each plot. Combining the subplots, which measured biomass less than 1 m above ground, with the microplot, which measured biomass above 1 m, we produced vertically stratified fuel load and foliage bulk density estimates from the ground surface to the top of the canopy. However, we did not measure 1000-h fuels because they were sparse in the study area and beyond the scope of this study.

Relating lidar metrics to conventional estimates of canopy bulk density

We used terrestrial lidar to estimate CBD by multiplying lidar-derived LAD by calibrated estimates of LMA_avg (Eqn 1). In our implementation, LAD, LMA_avg and CBD are vectors representing vertical profiles with a value for each height interval. To produce a predictive model, we created a vertical profile of LMA_avg for each vegetation type. Then, given a new lidar scan in this vegetation type, we can calculate a LAD profile, multiply by the LMA_avg profile and produce a lidar-predicted canopy bulk density profile. Our model assumes that LMA varies by species but is otherwise constant. Therefore, LMA_avg is determined by species composition along the vertical profile and the LMA value for each species (LMA_sp) (Eqn 2). We used the conventional CBD estimates to calculate species composition along the vertical profile for each vegetation type. To estimate LMA_sp, we combined prior estimates of LMA_sp from existing studies, our lidar-derived LAD profiles and our conventional estimates of CBD within a Bayesian regression framework (Eqn 3) to produce calibrated values for LMA_sp (Table 1).

In the Bayesian regression analysis, we adjusted LMA_sp to minimize bias of lidar-derived CBD estimates relative to conventional CBD estimates. Calibration was essential to account for assumptions made during LAD estimation and other biases inherent in both lidar and conventional CBD estimates. For prior coefficient estimates, we used a normal distribution centered on the previously published LMA values (Table 1). We used strongly informative priors to shrink modeled LMA estimates toward the previously published values, enhancing model stability and robustness. We scaled the prior standard deviations relative to the mean and iteratively decreased the standard deviations until predictive performance declined. The selected standard deviations were 5% of the means. A weak prior (σ ~ half-normal(0,1)) was specified for error terms. We estimated posterior distributions using the PyMC package (Abril-Pla et al. 2023) with default sampler settings. We then implemented a second stage of model fitting to address bias observed when summing the lidar-predicted CBD values to calculate total canopy fuel load (CFL). We summed the lidar CBD estimates and conventional CBD estimates to quantify CFL for each plot, then used simple linear regression to produce an adjustment factor. For each species group, we multiplied LMA_sp by this adjustment factor, ensuring that lidar-derived total CFL was unbiased relative to the conventional estimates.

We alternated plots through training and testing sets using k-fold cross-validation (k = 10), creating out-of-bag predictions for each plot. We removed two plots from training and testing that were unexplained outliers in both single- and multiple-scan mode. Using the lidar-predicted CBD profile for each plot, we extracted CBD_max, canopy base height and total CFL using the same methods that were applied to the conventional CBD profiles. We evaluated model performance by comparing lidar estimates and conventional estimates of CBD_max and total CFL. Preliminary testing showed that fitting LMA_sp estimates separately for each vegetation type did not substantially improve model performance. Therefore, we chose to fit one set of calibrated LMA_sp values across all vegetation types (Table 1) for model transferability. The models applied to each vegetation type differ only in species composition.

Relating lidar metrics to conventional estimates of surface fuel load

We also related lidar metrics to surface fuel loads less than 1 m above ground. For standing surface fuels, we related the mean LAD between 10 and 100 cm to the conventional estimates of standing surface fuel load using simple linear regression. This approach is similar to the approach applied to canopy fuels, but we did not have enough data to separate surface fuels by species. For downed surface fuels, we related LAI, defined here as the sum of LAD from 10 cm to the top of the canopy, to the conventional estimates of downed fuel load with a slope for each vegetation type. This approach simply highlights the relationship between accumulated dead fuel and the total density of live vegetation from which it falls. In our relationships between lidar metrics and surface fuels, we used plot-level averages instead of comparing lidar data with each individual subplot. We summed the predicted fuel load for standing and downed surface fuel to estimate the total surface fuel load.

Evaluating effects of scan mode, occlusion and terrain

After predicting fuel metrics, we calculated the residuals relative to the conventional estimates and analyzed the residuals with generalized additive models (GAMs) to evaluate the effects of environmental conditions. We visually analyzed the partial dependence plots to look for effects of occlusion, terrain metrics and height above ground on bias and variance. Terrain metrics included roughness, slope and concavity. We calculated occlusion by height bin and occlusion by plot using the percentage of voxels classified as occluded. Additionally, we compared model results between single- and multiple-scan mode and between vegetation types (Table 2).

Estimating fuel metrics from lidar data

Canopy estimates from single-scan lidar were significantly related to conventional estimates for total canopy fuel load (P < 0.001, R² = 0.744, RMSE = 0.092 kg/m²) (RMSE, root mean square error) and maximum canopy bulk density (P < 0.001, R² = 0.541, RMSE = 0.031 kg/m³) (Fig. 5). Fuel estimates using multiple-scan mode were mostly similar or slightly better (up to 15% lower RMSE) compared with single-scan estimates with a few notable exceptions (Table 2). Multiple-scan mode demonstrated a much stronger relationship with conventional estimates of maximum canopy bulk density (CBD_max) in PPO, with a 30% lower RMSE. In CHAP, multiple-scan mode produced a stronger R² (0.60), but also higher bias, leading to higher RMSE. Single-scan lidar produced better results in PJO, with a 24% lower RMSE in CBD_max and an 8% lower RMSE in total CFL. Relationships within each individual vegetation type were weaker than the relationships considering all vegetation types together (Table 2), particularly in PPO and CHAP. We also found strong agreement in the shape of CBD profiles from single-scan lidar compared with conventional methods (Fig. 6). Comparing the conventional profiles and lidar profiles, we can easily identify similar structural patterns and similar locations of peaks and dips, which describe the bimodal canopy structure in PPO, the shrub-heavy structure of CHAP and the intermediate structure of PJO. However, we found that the lidar-derived CBD profiles tended to show a narrower and more concentrated peak, leading to higher estimates of CBD_max relative to conventional estimates. Average bias in CBD_max across scan modes and vegetation types was 0.012 kg/m³.

Fig. 5.

Comparison of conventional estimates and single-scan lidar predictions for total canopy fuel load (left panel) and maximum canopy bulk density (CBD_max) (right panel). Solid line depicts 1:1 ratio.

Fig 6.

3D canopy structure from single-scan lidar (left), and comparison of canopy bulk density profiles using conventional and single-scan lidar methods (right) for three different vegetation types.

Downed surface fuel loads, composed of leaf litter and woody debris, were significantly related to single-scan LAI (P < 0.001, R² = 0.55, RMSE = 0.24 kg/m²) when considering all vegetation types together. However, each vegetation type had a distinctly different slope (Fig. 7) and relationships were much weaker when considering individual vegetation types (Table 2), particularly in PPO, which did not have a significant relationship between LAI and downed surface fuel load regardless of scan mode. Standing surface fuel loads, mostly comprising shrubs shorter than 1 m, were significantly related to surface LAD (P < 0.001, R² = 0.53, RMSE = 0.03 kg/m²) when considering all vegetation types. The relationship between near-surface LAD and standing vegetation appeared to be consistent across vegetation types. Differences in near-surface LAD were large between vegetation types but small within each vegetation type (Fig. 7), leading to weak relationships between near-surface LAD and standing surface fuel within each vegetation type. Using multiple-scan mode, RMSE was 10% lower for standing surface fuel and 1% lower for downed surface fuel.

Fig. 7.

Relationship between LAI and downed surface fuel load, consisting of litter and fine woody debris (left panel) and the relationship between near surface LAD (mean LAD between 0.1 and 1 m) and standing surface fuel load, consisting of live and dead rooted plants shorter than 1 m (right panel). Lidar metrics are from single-scan mode. In the left panel equation, CHAP, PJO and PPO are binary indicator variables signifying vegetation type.

Evaluating effects of scan mode, occlusion and terrain

Leaf area density profiles from single- and multiple-scan modes were closely related (Fig. 8). There was a strong linear relationship between the two data sources in LAD by height bin (R² = 0.96) and plot-level LAI (R² = 0.94), with LAD from single-scan mode being 19% lower than LAD from multiple-scan mode.

Fig. 8.

Comparison of single-scan versus multiple-scan estimates for LAD of individual height bins (left panel), and plot-level LAI (right panel). Solid line depicts fitted regression. Dashed line depicts 1:1 ratio.

In the single-scan data, on average, 20.6% of voxels were classified as occluded, whereas in the multiple-scan data, just 1.4% of voxels were occluded. Single-scan occlusion was highest in ponderosa pine–oak forest (29.4%) and similar between chaparral (16.4%) and pinyon–juniper–oak woodland (17.5%). Single-scan occlusion was also very high (>30%) near the ground in areas with dense shrub cover. When predicting the canopy bulk density of individual height bins, occlusion was associated with a positive bias and an effect of 0.035 kg/m³ between 0 and 35% occlusion. Height above ground was associated with a negative bias and an effect of −0.04 kg/m³ between 2 and 18 m as well as an increase in variance above 15 m. Increasing CBD was also associated with a negative bias and increasing variance. Terrain concavity (bowl-shaped terrain) may also have a positive bias in some cases. We found weak evidence for a positive bias in total CFL if plot-level occlusion is over 25%. We did not detect any other environmental effects on bias or variance when estimating total CFL, maximum CBD or surface fuel load. Additionally, we did not detect a consistent effect of terrain slope or roughness on bias or variance in any of our output variables. However, extreme terrain features, especially boulders and rock outcroppings, can contribute to occlusion near the ground surface.

Estimating fuel metrics from single-scan terrestrial lidar data

We demonstrate that LAD from single-scan lidar can be an effective tool for monitoring canopy and surface fuel across vegetation types with differing structure. Using single-scan terrestrial lidar and a Bayesian calibration process, our results were considerably stronger than those of Martin-Ducup et al. (2025), who used airborne lidar-derived LAD to produce a CBD profile and estimates of total CFL (R² = 0.18) and maximum CBD (R² = 0.28). Arkin et al. (2023) used mobile lidar scanning and area-based approach to achieve stronger relationships with conventionally derived CBD_max (R² = 0.71), but weaker relationships with total CFL (R² = 0.50). Our study prioritized estimation of total CFL (R² = 0.74) and accepted a weaker relationship with conventional estimates of maximum CBD (R² = 0.54). We prioritized total CFL because conventional estimates of CBD_max may not be reliable as training labels. Ex et al. (2016) and Keyser and Smith (2010) determined that the common-practice approach using a uniform function to distribute fuel along the length of a tree canopy (Sando and Wick 1972) is an inaccurate representation of fuel arrangement, which is better approximated by fitting a Weibull function. In those studies, the non-uniform function led to sharper peaks in the CBD profile and higher estimates of CBD_max. Our lidar analysis agreed with these findings, with lidar-informed CBD distributions showing sharper peaks and higher estimates of CBD_max compared with conventional estimates using a uniform fuel distribution function (Fig. 6). Fuels are often treated as uniform within a tree canopy because more detailed information is not available. With lidar and LAD techniques, there are opportunities to more accurately describe canopy fuel distribution and create more accurate estimates of CBD_max by considering the unique arrangement of fuel within each tree and stand.

Our results also highlight the importance of inter-species differences when relating lidar canopy metrics to CBD. Owing to differences in leaf traits and their interaction with lidar pulses, we estimated that pine and manzanita species contained three times more fuel than other species given the same lidar-derived LAD. Preliminary testing also showed major interspecies differences when using other lidar features such as volume metrics and point density metrics. This is important because ponderosa pine and manzanita are disproportionately distributed in the upper canopy and understory respectively. Accurate inference of CBD with lidar requires understanding where different species occur. We also found notable differences between published LMA values and our calibrated LMA values, implying that there are important interspecies differences in how leaves interact with lidar pulses. Ponderosa pine needles are long and sparsely distributed, which probably leads to overestimation of first-hit distance, overestimation of light transmittance and underestimation of LAD, which may explain why calibrated LMA estimates were much higher than expected (Table 1). Conversely, juniper has short, densely clustered needles that do not allow lidar pulses to pass through fine gaps, leading to underestimation of transmittance and overestimation of LAD, which may explain our unexpectedly low calibrated LMA estimates. Our considerations for species mix by height strata and calibration of LMA_sp may partially explain our stronger results compared with Martin-Ducup et al. (2025), who only considered the dominant species in each strata and did not calibrate LMA_sp. Additional model calibration could further improve LAD, LMA and CBD estimates by utilizing destructive reference samples, improved LAD estimation functions and correction factors for sensor characteristics, foliage characteristics and leaf/wood composition (Pimont et al. 2015, 2019; Soma et al. 2018).

Lastly, our results indicate potentially useful relationships between lidar-derived LAD metrics and surface fuels. Using a similar approach to canopy fuels, we found a promising relationship between single-scan LAD and conventional estimates of standing surface fuels (R² = 0.53, RMSE = 0.028 kg/m²), which were primarily shrubs in our study. Future research could test if these methods could also be applied to forbs and grasses. Future work could also be improved by destructive biomass sampling and more detailed analysis of LMA_sp. We also found relationships between LAI, a measure of total live foliage, and downed surface fuels (leaf litter and woody debris) in CHAP (P = 0.05, R² = 0.194) and PJO (P < 0.001, R² = 0.489). Our study area was undisturbed and quite stable; we expect relationships between surface fuel and canopy will be subject to change and greater variability depending on disturbance and successional state. Interspecies differences in leaf traits and their interaction with lidar pulses were also important to relationships between LAI and downed fuels, leading to model differences between vegetation types (Fig. 7).

Effects of scan mode, occlusion and terrain

The decision to use single-scan mode versus multiple-scan mode will depend on sampling design and vegetation characteristics. For stand-scale or landscape-scale applications, the simplicity and time efficiency of single-scan mode may justify increased error at each plot. Single-scan mode enables more sampling locations within the same timeframe and budget, potentially allowing field crews to evaluate a more extensive study area or install a higher plot density, capturing a wider range of spatial variability. More samples across a broader range of variability are better for training a larger-scale remote sensing model based on airborne lidar or satellite imagery, especially if non-linear machine-learning algorithms are used (Silva et al. 2020). However, the occlusion in single-scan mode also makes within-scan analysis difficult, particularly for techniques that require spatially complete 3D data, such as the object segmentation and classification techniques applied in Krisanski et al. (2021) or the fine-scale 3D fire behavior modeling applied in Marcozzi et al. (2023). For techniques based on LAD profiles, users of single-scan mode should estimate occlusion and evaluate potential bias if dense vegetation or other factors cause plot-level mean occlusion of more than 25%. Additional testing is also needed in tall and dense forest conditions where we found weaker results, partly because we did not have a broad gradient of forest density in our study area. Dense forests may pose a problem owing to increased bias and variance caused by occlusion. Tall trees may be a concern because lidar from the ground inherently results in higher scanning density near the bottom of tree canopies (Yin et al. 2022), which may contribute to the negative bias we observed when estimating CBD in the upper canopy. Aside from a potential positive bias in CBD prediction from concave terrain, we found that our predictions were robust against terrain characteristics. However, we still recommend caution in applying these techniques in extreme terrain because terrain features can contribute to occlusion.

Management implications

The workflow and framework presented here have many potential uses in evaluating potential fire risk, prioritizing treatment locations and evaluating treatment effects. Data from a single scan can provide rapid and repeatable estimates of CFL, bulk density, fuel strata gap (or canopy base height), total canopy height and total surface fuel load with lower potential for human bias than conventional methods. Combined with a user-assigned surface fuel type (Scott and Burgan 2005), our modeling framework provides all of the typical vegetation data required for 1D and 2D fire behavior modeling with tools such as BehavePlus (Heinsch and Andrews 2010) or FlamMap (Finney 2006). We could not compare the lidar and conventional methods of estimating canopy base height because only 5 of 67 plots had a canopy base height greater than 1 m. Beyond wildfire risk, the lidar-derived vegetation structure metrics can also be used to model forest inventory (Newnham et al. 2015; Molina-Valero et al. 2022), carbon storage (Sankey et al. 2013), wildlife habitat (Batchelor et al. 2023; Shokirov et al. 2023) and hydrological variables (Belmonte et al. 2021; O’Donnell et al. 2021).

The lidar processing and modeling techniques described here were designed to be open-source, efficient, repeatable and adaptable. After a model is trained, new predictions of canopy and surface fuels can be obtained quickly by providing filepaths and vegetation type to a single Python script. Processing is highly optimized and processing a new scan takes less than 30 s on a moderately powerful desktop computer (CPU: AMD Ryzen 7 5800X). Training a model for a new area can be performed non-destructively using our methods or by destructively scanning, clipping and weighing individual branches (Pimont et al. 2015). Although the data processing may appear complex, the number of parameters in our model is small and the relationships are linear. Therefore, this approach requires relatively few training plots (e.g. 10–30 per vegetation type) compared with models using an area-based approach, which typically have dozens or hundreds of parameters (Loudermilk et al. 2023) with high potential for non-linear relationships and interactive effects. Additionally, because our model is based on well-defined canopy measurements, a trained model is likely more intuitive, transferrable and reusable across a wider area than an area-based approach (Whelan et al. 2023). Over time, it may be possible for the broader terrestrial lidar community to compile look-up tables for LMA values by tree species, correction factors for various lidar scanners and species composition profiles for various forest types.