Prediction of Forest Canopy and Surface Fuels from Lidar and Satellite Time Series Data in a Bark Beetle-Affected Forest

Wildfire behavior depends on the type, quantity, and condition of fuels, and the effect that bark beetle outbreaks have on fuels is a topic of current research and debate. Remote sensing can provide estimates of fuels across landscapes, although few studies have estimated surface fuels from remote sensing data. Here we predicted and mapped field-measured canopy and surface fuels from light detection and ranging (lidar) and Landsat time series explanatory variables via random forest (RF) modeling across a coniferous montane forest in Colorado, USA, which was affected by mountain pine beetles (Dendroctonus ponderosae Hopkins) approximately six years prior. We examined relationships between mapped fuels and the severity of tree mortality with correlation tests. RF models explained 59%, 48%, 35%, and 70% of the variation in available canopy fuel, canopy bulk density, canopy base height, and canopy height, respectively (percent root-mean-square error (%RMSE) = 12–54%). Surface fuels were predicted less accurately, with models explaining 24%, 28%, 32%, and 30% of the variation in litter and duff, 1 to 100-h, 1000-h, and total surface fuels, respectively (%RMSE = 37–98%). Fuel metrics were negatively correlated with the severity of tree mortality, except canopy base height, which increased with greater tree mortality. Our results showed how bark beetle-caused tree mortality significantly reduced canopy fuels in our study area. We demonstrated that lidar and Landsat time series data contain substantial information about canopy and surface fuels and can be used for large-scale efforts to monitor and map fuel loads for fire behavior modeling at a landscape scale.


Introduction
Wildfires and other forest disturbances such as insect outbreaks and drought cause widespread forest mortality across the western United States [1][2][3]. Forest disturbance has increased in recent years, and is projected to continue to increase in coming decades, in part due to climate change [4][5][6][7][8].
Many aspects of wildfire behavior and effects are greatly influenced by forest fuels [9,10], which are shaped by prior disturbances such as insect epidemics, drought, fire, and disease. Spatially-explicit estimates of forest fuels and the effect of insect-caused tree mortality on fuels can serve as valuable information for forest managers and researchers in their efforts to predict and mitigate the effects of future fire.
Landsat time series data using random forest (RF) modeling; (2) applied RF models to produce fuel maps; and (3) compared fuel maps to an independent assessment of forest mortality to assess relationships between fuels and mortality severity. Our results showed that canopy and surface fuels can be predicted from lidar and Landsat time series data with moderate accuracy. Such predictions could be used to improve fuel maps, fire behavior modeling, and programs like LANDFIRE. We also showed that canopy fuels decreased 2-6 years after a severe bark beetle outbreak, and that many aspects of fuels did not appear to be strongly affected by variation in the severity of tree mortality within our study area.

Study Area
The study area covers 108,551 ha in eastern Grand County in north-central Colorado, USA ( Figure 1). An outbreak of mountain pine beetles (Dendroctonus ponderosae Hopkins) began in the study area around 2003 and was subsiding by 2010, with most tree mortality occurring from 2004-2008 ( Figure 2) [68]. Lodgepole pine (Pinus contorta Dougl. ex Loud.), a principal host of mountain pine beetles, is the predominant tree species across the study area; Engelmann spruce (Picea engelmannii Parry ex Engelm) and subalpine fir (Abies lasiocarpa (Hook.) Nutt.) occur at higher elevations and aspen (Populus tremuloides Michx.) is also present throughout. Elevation averages 2800 m and ranges from 2370-3755 m above sea level. Average annual mean temperature is 1.6 °C and average annual precipitation is 51.4 cm (1961-1990) [69].    [68] for 30-m forested pixels across the study area.

Field Observations
In the summer of 2010 tree and surface fuel observations were gathered for 119 plots across the study area ( Figure 1) [70]. The location of each plot was recorded with a Trimble ProXH global positioning system (GPS) and post-processed to differentially correct the locations. Horizontal accuracy of the post-processed GPS locations ranged between 0.1 and 0.7 m and averaged 0.3 m. Tree measurements included diameter at breast height (DBH, 1.37 m), tree height, tree crown base height, crown canopy dominance (dominant, codominant, intermediate, suppressed), health, position relative to plot center, and species of each tree within an 8.0-m radius of plot center (plot area = 200 m 2 ). For each beetle-killed tree, the number of years since death was estimated based on the proportion of needles and branches remaining [63,65]. Height and canopy base height were not measured on some trees. If not measured, height was estimated from DBH using the allometric equations of [71], and canopy base height was assumed to be 50% of tree height [72].
Surface fuels were measured along three 20-m transects [72,73] and plot-level fuel loads were summarized by size class (1-, 10-, 100-, and 1000-h fuels corresponding to 0-0.6, 0.6-2.5; 2.5-8, and ≥8 cm diameters, respectively) [74]. In addition to downed dead wood, fuel loads for duff, litter, and herbaceous and shrub components were estimated within two 1-m 2 radius plots along each of the 3 transects and then summarized to the plot level [72]. In our analyses, response variables for surface fuels included loads of litter and duff, 1-100 h downed dead wood, and 1000-h downed dead wood.

Canopy Fuel Estimation
Canopy fuels for each plot were calculated in R [75] following the methodology of Scott and Reinhardt [11] as implemented in the Fire and Fuels Extension of the Forest Vegetation Simulator [76] and FuelCalc [12,77] software packages. For both live and dead trees, total crown weight and crown weight contained in foliage and small branches (0-6 mm in diameter) were estimated from DBH measurements using the species-specific allometric equations of Brown [78]. Following Lutes et al. [12], crown weights were adjusted based on crown canopy dominance, so that crown weights of intermediate and suppressed trees were reduced.
For live trees, foliage and 50% of small branch weight was considered available canopy fuel (ACF) [11]. Because the beetle outbreak had been ongoing in the study area for seven years and beetlekilled trees lose needles 2-4 years after attack [52,53], little or no foliage remained on the majority of beetle-killed trees measured in the field; Landsat imagery confirmed this, showing that the forest was largely in the gray phase of attack, in which killed trees have lost most red needles. Trees killed one, two, three, and four years previously were assumed to have 100%, 75%, 25%, and 0% of needles, respectively, and 100% of fine branches available as canopy fuel [61,63,65,79]. Trees killed five or more years previously were assumed to have no ACF. Tree-level ACF was aggregated to the plot level to be used as a response variable.

Field Observations
In the summer of 2010 tree and surface fuel observations were gathered for 119 plots across the study area ( Figure 1) [70]. The location of each plot was recorded with a Trimble ProXH global positioning system (GPS) and post-processed to differentially correct the locations. Horizontal accuracy of the post-processed GPS locations ranged between 0.1 and 0.7 m and averaged 0.3 m. Tree measurements included diameter at breast height (DBH, 1.37 m), tree height, tree crown base height, crown canopy dominance (dominant, codominant, intermediate, suppressed), health, position relative to plot center, and species of each tree within an 8.0-m radius of plot center (plot area = 200 m 2 ). For each beetle-killed tree, the number of years since death was estimated based on the proportion of needles and branches remaining [63,65]. Height and canopy base height were not measured on some trees. If not measured, height was estimated from DBH using the allometric equations of [71], and canopy base height was assumed to be 50% of tree height [72].
Surface fuels were measured along three 20-m transects [72,73] and plot-level fuel loads were summarized by size class (1-, 10-, 100-, and 1000-h fuels corresponding to 0-0.6, 0.6-2.5; 2.5-8, and ≥8 cm diameters, respectively) [74]. In addition to downed dead wood, fuel loads for duff, litter, and herbaceous and shrub components were estimated within two 1-m 2 radius plots along each of the 3 transects and then summarized to the plot level [72]. In our analyses, response variables for surface fuels included loads of litter and duff, 1-100 h downed dead wood, and 1000-h downed dead wood.

Canopy Fuel Estimation
Canopy fuels for each plot were calculated in R [75] following the methodology of Scott and Reinhardt [11] as implemented in the Fire and Fuels Extension of the Forest Vegetation Simulator [76] and FuelCalc [12,77] software packages. For both live and dead trees, total crown weight and crown weight contained in foliage and small branches (0-6 mm in diameter) were estimated from DBH measurements using the species-specific allometric equations of Brown [78]. Following Lutes et al. [12], crown weights were adjusted based on crown canopy dominance, so that crown weights of intermediate and suppressed trees were reduced.
For live trees, foliage and 50% of small branch weight was considered available canopy fuel (ACF) [11]. Because the beetle outbreak had been ongoing in the study area for seven years and beetle-killed trees lose needles 2-4 years after attack [52,53], little or no foliage remained on the majority of beetle-killed trees measured in the field; Landsat imagery confirmed this, showing that the forest was largely in the gray phase of attack, in which killed trees have lost most red needles. Trees killed one, two, three, and four years previously were assumed to have 100%, 75%, 25%, and 0% of needles, respectively, and 100% of fine branches available as canopy fuel [61,63,65,79]. Trees killed five or more years previously were assumed to have no ACF. Tree-level ACF was aggregated to the plot level to be used as a response variable.
For each tree, we assumed that ACF was distributed throughout 0.3-m height strata along the length of the tree crown, defined as the distance between tree crown base height and tree height. For most species, ACF was assumed to be evenly distributed throughout the canopy. However, for lodgepole pine, the majority of field-measured trees, ACF was distributed unevenly using an equation from Lutes et al. [12]: where W is the proportion of ACF weight for height percentile H. Tree-level AFC in each 0.3048-m height strata was summed by plot, converted to density units of kg m −3 , and smoothed using a 1.5-m running mean to create a canopy density profile for each plot ( Figure 3). Canopy bulk density (CBD) was defined as the maximum of the canopy density profile; plot canopy base height (CBH) and plot canopy height (CH) were defined as the minimum and maximum heights, respectively, at which the canopy density profile dropped below 0.012 kg m −3 [11]. CBD, CBH, and CH were used as response variables in analyses. For each tree, we assumed that ACF was distributed throughout 0.3-m height strata along the length of the tree crown, defined as the distance between tree crown base height and tree height. For most species, ACF was assumed to be evenly distributed throughout the canopy. However, for lodgepole pine, the majority of field-measured trees, ACF was distributed unevenly using an equation from Lutes et al. [12]: where W is the proportion of ACF weight for height percentile H. Tree-level AFC in each 0.3048-m height strata was summed by plot, converted to density units of kg m −3 , and smoothed using a 1.5-m running mean to create a canopy density profile for each plot ( Figure 3). Canopy bulk density (CBD) was defined as the maximum of the canopy density profile; plot canopy base height (CBH) and plot canopy height (CH) were defined as the minimum and maximum heights, respectively, at which the canopy density profile dropped below 0.012 kg m −3 [11]. CBD, CBH, and CH were used as response variables in analyses.

Lidar
Lidar data were acquired across the study area in August 2010 by Aerometric, Inc. using an Optech 3100 AE lidar system aboard a fixed-winged aircraft. Flying height above ground and flight speed of the aircraft during the lidar acquisition were 1250 m and 160 kts, respectively. Side lap of the acquisition was 50%, with a scan angle of ±22° from nadir. Nominal pulse spacing was 0.7 m, or about 2 pulses m −2 . Aerometric, Inc. delivered point cloud data in the LAS file format.
Delivered lidar data were processed with LAStools software [80]. Returns with scan angles >18° and <−18° from nadir were thinned to ensure that only high-quality returns were used for analysis. The remaining returns were classified as ground or nonground, and nonground return heights were normalized to heights above ground. Lidar metrics coincident with 8-m-radius plot extents and lidar metric grids at 30-m resolution were created with the 'lascanopy' function using all returns (Table 1). Metric grids were created with 30-m resolution to be comparable to Landsat data products. A 30-m meter digital terrain model (DTM) was created from ground-classified returns; this DTM was input

Lidar
Lidar data were acquired across the study area in August 2010 by Aerometric, Inc. using an Optech 3100 AE lidar system aboard a fixed-winged aircraft. Flying height above ground and flight speed of the aircraft during the lidar acquisition were 1250 m and 160 kts, respectively. Side lap of the acquisition was 50%, with a scan angle of ±22 • from nadir. Nominal pulse spacing was 0.7 m, or about 2 pulses m −2 . Aerometric, Inc. delivered point cloud data in the LAS file format.
Delivered lidar data were processed with LAStools software [80]. Returns with scan angles >18 • and <−18 • from nadir were thinned to ensure that only high-quality returns were used for analysis. The remaining returns were classified as ground or nonground, and nonground return heights were normalized to heights above ground. Lidar metrics coincident with 8-m-radius plot extents and lidar metric grids at 30-m resolution were created with the 'lascanopy' function using all returns (Table 1). Metric grids were created with 30-m resolution to be comparable to Landsat data products. A 30-m meter digital terrain model (DTM) was created from ground-classified returns; this DTM was input into a DTM derivative tool for ERDAS IMAGINE [81] to create 30-m resolution topographic variable grids. Values of the DTM and derivative grids were extracted at the 8-m-radius plots to be used as explanatory variables along with lidar metrics. Lidar metric grids and DTM derivative grids were created for mapping. Nonground returns >0.15 and <0.5 m in height D02 1 Nonground returns >0.5 and <1.37 m in height D03 1 Nonground returns >1.37 and <5 m in height D04 1 Nonground returns >5 and <10 m in height D05 1 Nonground returns >10 and <20 m in height D06 1 Nonground returns >20 and <30 m in height 1 Each density strata included five variables: average return height (AVG), standard deviation of return heights (STD), skewness of return heights (SKE), and kurtosis of return heights (KUR) within that strata, as well as the proportion of returns within that strata (PROP).

Landsat Time Series
Landsat surface reflectance imagery of Path 34, Row 32 for years 1984-2010 was processed using the Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) set of algorithms [17]. LandTrendr fits straight-line segments to annual spectral trajectories through time, pixel-by-pixel. Generalized pixel trajectories are then described by a set of variables and output as grids (Table 2). For our analysis, the normalized difference ratio of Landsat bands 4 and 7 (band 4 − band 7)/(band 4 + band 7), also known as the normalized burn ratio (NBR) [82] was used to create generalized pixel trajectories and associated variables. A decrease in NBR indicated disturbance; an increase post-disturbance indicated recovery. Generalized trajectories were also imposed upon Landsat tasseled-cap bands of brightness (TCB), greenness (TCG), and wetness (TCW) [83] and output as additional grids. Values of LandTrendr output grids coincident with 8-m-radius plot extents were extracted for use as explanatory variables.
The Landsat time series product of Meddens et al. [84], which quantifies annual bark beetle-caused forest mortality across our study area for years 1996-2011, was also used to create explanatory variables. Meddens et al. [84] describes the creation of this product in detail. Briefly here, Landsat top-of-atmosphere images of Path 34 Row 32 from late summer for years 1996-2011 were downloaded and radiometrically normalized. For each annual image, vegetated pixels were separated from non-vegetated pixels with Landsat band 2 (band 2 < mode of band 2 indicated vegetation); vegetated pixels were further separated into non-forest and forest pixels with Landsat TCG (TCG <~1.1 indicated forest); forest pixels were further separated into undisturbed and beetle-disturbed pixels based on anomalies from the undisturbed mean of the Landsat band 5/band 4 ratio (B5/B4). This method was further expanded in Meddens and Hicke [68] with a multiple linear model that predicted percent beetle-killed tree cover from Landsat band 1, the normalized difference vegetation index (NDVI; (band 4 -band 3)/(band 4 + band 3)), and the B5/B4 (band 5/band 4) anomaly. This model explained 77% of the variability in percent red tree cover and was applied to pixels classified as beetle-disturbed, resulting in annual grids of percent red tree cover across our study area. For the purpose of our study, we excluded the 2011 mortality grid and summarized annual grids into three explanatory grids: cumulative percent tree mortality cover in 2010 (MORT%), bark beetle mortality duration (number of years new red tree mortality cover was detected, MORTdur), and bark beetle mortality rate (MORT%/MORTdur; MORTrate) ( Table 2). We used the values of these grids coincident with each of our 8-m-radius plot extents as explanatory variables in analyses.

Aerial Detection Survey
The U.S. Forest Service conducts annual aerial detection surveys (ADS) of insect-caused forest disturbance [85]. Observers aboard fixed-wing aircraft record the extent, severity, and cause of forest damage within digitized polygon boundaries. Meddens et al. [3] converted this polygon ADS data to a 1-km gridded dataset of bark beetle-caused forest mortality cover for years 1997-2010 across the United States; middle estimate grids were later added to original lower and upper estimate grids. We summed 1-km grids of mortality cover (middle estimates) for years 1997-2010 coincident with our study area to create a 1-km grid of cumulative mortality cover in 2010 across our study area [86] (see Figure 2 of [86]). This 1-km ADS grid served as an independent measure of mortality with which to compare our mapped fuel predictions.

Random Forest Modeling
Random forest (RF) modeling has become a widely-used methodology for predicting response variables from explanatory variables [87], including the prediction of forest attributes, e.g., [88][89][90]. Advantages of RF modeling include: (1) its nonparametric nature, making response variable transformation to a normal distribution and subsequent back transformation unnecessary; (2) the ability to perform both classification and regression tasks; (3) the internal partitioning of data into training and evaluation datasets; (4) relatively simple model parameterization; and (5) the generation of variable importance scores that allows for the ranking and selection of explanatory variables so that parsimonious models can be created.
Canopy and surface fuel response variables were predicted from lidar and Landsat time series explanatory variables via RF modeling implemented in the 'randomForest' package of R [87,91]. Out-of-bag sampling inherent in the RF algorithm was used to develop models and evaluate errors. RF models considering only lidar explanatory variables and considering both lidar and Landsat time series explanatory variables were developed. For each response variable and set of explanatory variables (lidar only and lidar plus Landsat), variable and model selection were performed with the 'rf.modelSel' function of the 'rfUtilities' R package [92]. In this approach, all candidate explanatory variables are assigned model improvement ratio (MIR) scores, a standardized measure of the mean decrease in accuracy when a variable is withheld. Several models including variables above iterative, increasing MIR thresholds are produced (10 models in our analyses), and the model that maximizes percent variance explained in the response variable using the most important and fewest number of explanatory variables is then identified as the final model, e.g., [93]. If a pair of highly correlated (Pearson's r > 0.9) explanatory variables was selected, the variable with a smaller importance score was discarded and the model selection tool was rerun. Spearman's correlations between selected explanatory variables and response variables were also calculated. Final RF models were applied to explanatory variable grids to map the response variables at a resolution of 30-m with the 'raster' R package [94] and the 'AsciiGridPredict' function of the 'yaImpute' R package [95]. Mapped response variables were then aggregated, taking the mean, to the resolution of 1-km ADS grid cells that quantified forest mortality cover. Relationships between the values derived from the aggregated fuel maps and the values derived from the ADS estimates of forest mortality were assessed via Spearman's correlation tests.

Results
Across the 119 8-m-radius plots, available canopy fuel averaged 7.3 Mg ha −1 and canopy bulk density averaged 0.09 kg m −3 (Table 3). Canopy base height and canopy height averaged 2.5 and  RF models explained 28-70% of the variation in canopy fuel variables (Table 4); including Landsat time series variables in the models increased the percent of variance explained by 2-7%. Surface fuel variables were predicted less accurately, RF models explaining 16-32% of variance. The addition of Landsat time series variables to the models increased the percent of variance explained in surface fuels by 2-8%. The greater predictive capacity of lidar variables relative to satellite variables is not surprising, as lidar directly senses the three-dimensional structure of the canopy and to a lesser extent, the subcanopy and forest floor. In addition, the lidar-generated variables have a higher spatial resolution (generated from multiple returns per m 2 ) compared to the multispectral Landsat data (30-m spatial resolution), adding to the increased predictive power of the lidar data. Surface fuel variables were predicted with less accuracy than canopy variables because both airborne lidar and satellite image sensors are more sensitive to fuel structure in the forest overstory canopy than in the forest floor. However, lidar and satellite sensors did provide some information about the forest floor, collectively explaining 24-32% of variation in surface fuel variables.
Of the lidar-derived metrics, canopy density variables, especially the percentage of returns >1.37 m in height (DNS) and the proportion of returns >10 and <20 m in height (D05 PROP), were important predictors of canopy fuel variables ( Table 5). The skewness of canopy returns heights (SKE) was the most important predictor of CBH, and maximum canopy return height (MAX) was the most important predictor of CH. SKE and MAX are logical selections for describing CBH and CH, respectively; SKE is a measure of canopy profile shape from which CBH is derived (Figure 2), and MAX closely approximates CH. Elevation above sea level (ELEV) helped explain variation in ACF and CH; both were positively correlated with ELEV. Only three LandTrendr variables were selected as important predictors of canopy fuel variables; pre-disturbance Landsat tasseled cap wetness (GDpreTCW), a measure of vegetation moisture content, was an important predictor of both ACF and CBD. Percent cover of bark beetle-caused tree mortality (MORT%) was selected as an important predictor of-and negatively correlated with-ACF. Both MORT% and the rate of bark beetle-caused tree mortality (MORTrate) helped explain variation in CBH.
Similar to canopy fuel variables, lidar density variables, both lower strata variables and upper strata variables were important predictors of surface fuels ( Table 5). The selection of both lower and upper density strata explanatory variables in the prediction of surface fuels suggests that lidar was sensitive to surface fuels both directly and indirectly. In particular, the fact that the skewness of returns >0.15 and <0.5 m in height (D01 SKE) was chosen as an important predictor of surface fuels is evidence that lidar directly interacted with surface fuels, whereas the selection of an upper canopy variable, the 75th percentile of canopy return height (P75), as the most important predictor of 1-100-h and total surface fuels shows that this variable described surface fuels indirectly through association with canopy conditions. Landsat time series variables were chosen as important predictors of surface fuels more frequently than they were chosen for canopy fuel variables, likely because lower performance of lidar variables left substantial variation for the satellite variables to explain. The most important predictor of litter and duff was post-disturbance canopy cover (LDpost.val), which presumably represents an indirect measure of litter and duff as they accumulate below dead trees. Pre-disturbance Landsat tasseled cap greenness (LDpreTCG), an indicator of vegetation density [96], was a similar indirect measure of fine surface fuels. Pre-disturbance Landsat tasseled cap brightness (LDpreTCB), sensitive to soil characteristics, might have been a more direct indicator of fine surface fuels. Pre-disturbance Landsat tasseled cap wetness (GDpreTCW) and the duration of the greatest disturbance (GDdur) were important predictors and positively correlated with 1000-h and total surface fuel loads. GDdur described both the recent bark beetle disturbance and previous disturbance, which LandTrendr detected across much of the study area ( Figure 4). MORT% was an important predictor of 1000-h surface fuels, the relationship between 1000-h surface fuels and MORT% being negative (ρ = −0.4). Table 5. Selected explanatory variables, variable importance scores (0-1), and Spearman's correlation values (ρ, given in parentheses) of random forest models predicting canopy and surface fuels from lidar and Landsat variables. Variable importance scores of one indicate most important, zero least important. Spearman's correlation values range from −1 to 1; negative and positive values indicate negative and positive correlation between response and explanatory variables, respectively. Explanatory variable definitions: maximum canopy height (MAX); standard deviation of canopy return heights (STD); skewness of canopy return heights (SKE); kurtosis of canopy return heights (KUR); 25th percentile of canopy return heights (P25); 75th percentile of canopy return heights (P75); percentage of all returns >1.37 m in height (DNS); proportion of nonground returns <0.15 m in height (D00 PROP); skewness of nonground returns >0.15 and <0.5 m in height (D01 SKE); proportion of nonground returns >0.5 and <1.37 m in height (D02 PROP); proportion of nonground returns >1.37 and <5 m in height (D03 PROP); proportion of nonground returns >10 and <20 m in height (D05 PROP); average return height of returns >10 and <20 m in height (D05 AVG); skewness of returns >10 and <20 m in height (D05 SKE); proportion of nonground returns >20 and <30 m in height (D06 PROP); average return height of returns >20 and <30 m in height (D06 AVG); standard deviation of returns >20 and <30 m in height (D06 STD); kurtosis of returns >20 and <30 m in height (D06 KUR); elevation above sea level in meters (ELEV); curvature, the second derivative of slope (CURV) [97]; duration of greatest disturbance (GDdur); recovery duration after greatest disturbance (GDpost.dur); Landsat tasseled cap wetness before greatest disturbance (GDpreTCW); cover value after greatest disturbance (LDpost.val); Landsat tasseled cap brightness before longest disturbance (LDpreTCB); Landsat tasseled cap greenness before longest disturbance (LDpreTCG); longest change in Landsat tasseled cap wetness (LD∆TCW); percent cover of bark beetle-caused tree mortality in 2010 (MORT%); rate of bark beetle-caused tree mortality (MORTrate).  Two topographic variables, ELEV and curvature, the second derivative of the slope (CURV) [97], were selected as important predictors of fine surface fuels. ELEV was negatively correlated with 1-100-h surface fuels. The relationship between CURV and litter/duff is unclear; however, CURV describes slope shape, whether upwardly convex (positive CURV) or upwardly concave (negative CURV), which might affect the degree to which litter and duff are subject to erosion.
Maps of fuel response variables show mortality-related patterns (Figures 1 and 5) that were confirmed through Spearman's correlation tests using independent ADS-derived mortality data (Table 6). In forested areas affected by bark beetles, ACF, CBD, and CH decreased with increasing tree mortality severity, while CBH increased with increasing tree mortality severity. Though weak, most relationships were significant.
Relationships between surface fuels and mortality severity were predominantly negative, somewhat contrary to expectations. No significant relationship between litter and duff and ADSderived mortality severity was found. A weak negative relationship between 1-h to 100-h surface fuels and mortality severity existed, as well as a moderate negative relationship between 1000-h surface fuels and mortality severity. Coarse surface fuels were likely greater in areas of less severe bark beetle mortality at the time of our sampling (2-6 years after peak epidemic mortality in the region) because these areas had experienced previous disturbance, as revealed by LandTrendr (Figure 4), resulting in reduced forest density and lower availability of host trees for bark beetles in the 2000s epidemic. Two topographic variables, ELEV and curvature, the second derivative of the slope (CURV) [97], were selected as important predictors of fine surface fuels. ELEV was negatively correlated with 1-100-h surface fuels. The relationship between CURV and litter/duff is unclear; however, CURV describes slope shape, whether upwardly convex (positive CURV) or upwardly concave (negative CURV), which might affect the degree to which litter and duff are subject to erosion.
Maps of fuel response variables show mortality-related patterns (Figures 1 and 5) that were confirmed through Spearman's correlation tests using independent ADS-derived mortality data (Table 6). In forested areas affected by bark beetles, ACF, CBD, and CH decreased with increasing tree mortality severity, while CBH increased with increasing tree mortality severity. Though weak, most relationships were significant.
Relationships between surface fuels and mortality severity were predominantly negative, somewhat contrary to expectations. No significant relationship between litter and duff and ADS-derived mortality severity was found. A weak negative relationship between 1-h to 100-h surface fuels and mortality severity existed, as well as a moderate negative relationship between 1000-h surface fuels and mortality severity. Coarse surface fuels were likely greater in areas of less severe bark beetle mortality at the time of our sampling (2-6 years after peak epidemic mortality in the region) because these areas had experienced previous disturbance, as revealed by LandTrendr (Figure 4), resulting in reduced forest density and lower availability of host trees for bark beetles in the 2000s epidemic.  Areas where canopy cover is <10% are masked in white. Water bodies are black. Maps were created by applying random forest models to explanatory variable grids.

Discussion
We predicted canopy and surface fuel loads with moderate accuracy in a lodgepole pine-dominated forest affected by an extensive bark beetle outbreak. When predicting these fuel metrics, lidar-derived variables explained much more variation than did multispectral Landsat variables, and Landsat variables added minor to moderate improvements to the predictive capacity of lidar-based models. Similar findings have been reported for predicting forest stand structure attributes such as basal area [86]. We found that the main effects of the bark beetle outbreak shortly (2-6 years) after tree mortality in our study area were a decrease in canopy height and canopy fuels, with an increase in CBH. The maps we generated of these metrics could be used by land managers, policy makers, and fire fighters to assess the impacts of bark beetle disturbances on landscapes that are generally perceived as highly flammable and hazardous [98].
We explained less variation in canopy fuels than others have in previous studies that also used lidar, and in some cases multispectral imagery, to estimate continuous canopy fuel variables in undisturbed forests. Andersen et al. [21] explained 86%, 84%, 77%, and 98% of the variation in sqrt(ACF), ln(CBD), CBH, and CH, respectively, in a Douglas-fir forest. For the same forest type, Hermosilla et al. [28] reported R 2 values of 79%, 67%, 78%, and 79% for ACF, CBD, CBH, and CH, respectively. In a ponderosa pine forest, Hall et al. [22] achieved R 2 values of 83% and 80% for CBD and CBH, respectively. Popescu and Zhao [24] and Vauhkonen [26] explained 80% and 84% of the variation in CBH of individual trees. Erdody and Moskal [34] included both lidar and aerial imagery explanatory variables to achieve R 2 values of 90, 87, 83, and 95% for ACF, CBD, CBH, and CH, respectively, in a mixed conifer forest. Skowronski et al. [27] predicted ACF and CBD in a pitch pine forest with R 2 values of 71 and 83%, respectively. Using ICESat and GLAS satellite lidar data, García et al. [32] explained 78% of the variation in CBD. Jakubowski et al. [36] predicted CBD, CBH, and CH with R 2 values of 25, 41, and 87%, respectively, in a mixed conifer forest. In a Pinus radiata forest in Spain, González-Ferreiro et al. [29] predicted ACF, CBD, CBH, and CH with R 2 values of 82%, 44%, 98%, and 98%, respectively, using low-density lidar.
There are a few reasons why we explained less variation in fuels than others have in previous studies. First, previous studies have not predicted canopy fuels in lodgepole pine forests, which, in contrast to other forest types previously studied, tend to be denser, with closed canopies and relatively shorter, uniform tree heights in many stands leading to less variability available to explain in our data. Using higher density lidar that is better able to penetrate the forest canopy or adjusting for the attenuation of the lidar from the canopy to the surface [99] might improve the sensitivity of lidar to variation in fuels, especially surface fuels, in these relatively shorter and denser lodgepole pine forests. Second, we might have explained more variation in fuels if our field plots had been larger; preliminary models that used the field plot data of Meddens et al. [100], which had larger 400-m 2 plot extents, performed comparably to those of previous studies [101,102]. We chose not to use this field data, however, because it had no associated surface fuel measurements and was gathered two years previous to our lidar data in a forest that was undergoing rapid change. Third, we did not include lidar intensity although our preliminary models also showed lidar intensity (amount of energy reflected back to the sensor) variables to be important predictors that improved model performance; however, in this study we were not able to normalize intensities so intensity information was withheld from analysis. Normalization of lidar intensity values is necessary to control for the effect of differing range between the plane and location of returns, and automatic gain control (AGC) of the lidar instrument during acquisition; intensities could have been normalized if we had had access to information on the smoothed best estimate of trajectory (SBET) of the plane and AGC information from the lidar instrument.
Few previous studies have estimated surface fuels from lidar or satellite data. The percentages of variance explained by our models-25-32% for surface fuel variables-are similar to those of previous studies. Pesonen et al. [45] explained 61% of the variation in downed dead wood volume using airborne lidar. Using airborne SAR and optical data, Huang et al. [46] achieved moderate correlation (R = 0.54) between measured and predicted coarse woody debris. Meigs et al. [47] explained 29% of the variation in coarse woody debris with LandTrendr disturbance magnitude. Using lidar and multispectral variables in a mixed conifer forest, Jakubowski et al. [36] attained R 2 values of 32% and 48% for 1000-h and total surface fuels, respectively, and also explained 35% of the variation in fuel bed depth. Hudak et al. [49] explained 32-44% of the variation in surface fuels in an open longleaf pine forest. Price and Gordon [31] explained 24% of the variation in surface fuel loads in dry sclerophyll forests in Australia. Like canopy fuel variables, we found lidar intensity variables to be effective predictors of surface fuel variables in preliminary models. Future attempts to predict surface or canopy fuels from lidar might benefit from the inclusion of normalized lidar intensity information [103].
We found that canopy fuel responses to bark beetle mortality generally conformed to expectations and previous research. Decreases in ACF and CBD caused by the loss of needles and fine branches during the gray stage are well-documented [56,[58][59][60][61]. However, to our knowledge no previous study has demonstrated the negative correlations we found between ACF, CBD and severity of tree mortality across a landscape. Previous studies are less clear about the effect that tree mortality has on CBH, with some reporting a decrease in CBH [58,62] or no change in CBH [51,56,60]. Like Klutsch et al. [59] and Donato et al. [61], we found that CBH increased following mortality, and demonstrated that CBH was positively and significantly correlated with tree mortality severity across the landscape. A post-outbreak increase in CBH is logical, as a decrease in overall CBD would cause the height of the critical lower threshold of 0.012 kg m −3 that defines CBH to increase [61]. For the same reason, the decrease in CH with increasing mortality severity that we found is also logical, as a decrease in ACF and CBD in the upper canopy would decrease the height of the upper critical threshold that defines CH.
Relationships between surface fuels and mortality that we found should be treated with less certainty, as percent variance explained by surface fuel models was fairly low (Table 4). Low model performance could be due not just to lack of satellite and lidar data sensitivity but also by the inherent high variability in fine surface fuel distributions on the ground, making them difficult to sample and predict [104]. Variables selected as important for the prediction of surface fuels suggest that lidar and Landsat data measured surface fuels both directly and indirectly, i.e., surface fuels were indirectly predicted through their association with canopy characteristics [44,50,105]. We found no significant relationship between litter and duff and mortality severity, and a weak negative relationship between 1-100-h surface fuels and mortality severity. Although a positive relationship between fine surface fuels and mortality severity might be expected due to bark beetle-killed trees dropping dead needles and small branches [51,60,[62][63][64], others have also found bark beetle mortality to have no significant effect on fine surface fuels in the gray stage [56,59,61], and have not detected significant relationships between fine surface fuels and satellite-derived mortality measures [47]. While trees killed in a beetle outbreak appear to contribute a large single pulse of needles and fine branches to the forest floor, this phenomenon only occurs over several years in a given area during the course of an epidemic, and may not strongly influence the distribution of fine surface fuels relative to contributions from live trees over longer periods of time.
The negative relationship between coarse surface fuels and mortality severity that we found was somewhat unexpected. Others have reported no significant increase in coarse surface fuels in the gray stage before snag fall [56,59,61,63], and Meigs et al. [47] similarly found no significant relationship between bark beetle outbreak duration and coarse woody debris. The outbreak in our study area began only seven years previous to our field and remote sensing measurements, with most mortality occurring 2-6 years prior to measurement. Snag fall timing is variable, generally beginning 3-5 years post attack, with low initial snag fall rates that increase with time [65][66][67]; however, no studies have yet measured fall rates in the 2000s mountain pine beetle epidemic in Rocky Mountain lodgepole pine forests. In our study area, most beetle-killed trees were still standing and most 1000-h fuels were from prior disturbance. Stands regenerating from or thinned by previous disturbance were likely less vulnerable to bark beetles in the 2000s epidemic and lost fewer trees to beetles, hence the negative relationship that we found between 1000-h fuels and mortality severity. Downed coarse woody debris should accumulate in our study area as beetle-killed snags fall, eventually reversing this negative relationship [51,55,[59][60][61][62][63].
Previous studies disagree on how bark beetle-caused tree mortality influences subsequent fire severity. Studies have reported a decrease [106,107], no change [108][109][110][111][112], an increase [113,114], and a mixed response [115,116] in fire severity following bark beetle-caused tree mortality. The extent and severity of bark-beetle caused tree mortality can vary greatly at the landscape scale, e.g., [117][118][119]; how this variability in severity of mortality influences fire severity is unclear [51], although time since attack is an important factor. The decreases in ACF and CBD and the increase in CBH that we found six years after a bark beetle outbreak in our study area suggest a corresponding decrease in potential crown fire, although additional research on the relationship between beetle-caused tree mortality and fire is still needed, given the complex temporal and spatial variation in these two processes.

Conclusions
We demonstrated the use of lidar and Landsat time series data for predicting canopy and surface fuels in a forest heavily affected by bark beetle mortality. Lidar data, and to a lesser degree Landsat time series data, can add great value to spatially-limited plot-based information in situations where disturbances have affected extensive areas. Where field measurements and lidar are available elsewhere, similar methods could be used to map canopy and surface fuels and describe relationships between fuels and other types of forest disturbance. Where repeat lidar is available, change in fuel loads could be assessed.
Using lidar data alone resulted in models that were nearly as informative as those that included both lidar and Landsat time series data. Including lidar data in fuel mapping has the potential to improve fire behavior modeling and existing methodologies for mapping fuels in extensive spatial databases such as LANDFIRE [120,121] that rely upon information from passive satellite sensors; lidar actively senses three-dimensional vegetation structure thereby providing more direct and finer-scale measurements of vegetation and fuels.