Prediction of Fuel Loading Following Mastication Treatments in Forest Stands in North Idaho, USA

: Fuel reduction in forests is a high management priority in the western United States and mechanical mastication treatments are implemented common to achieve that goal. However, quantifying post-treatment fuel loading for use in ﬁre behavior modeling to forecast treatment e ﬀ ectiveness is di ﬃ cult due to the high cost and labor requirements of ﬁeld sampling methods and high variability in resultant fuel loading within stands after treatment. We evaluated whether pre-treatment LiDAR-derived stand forest characteristics at 20 m × 20 m resolution could be used to predict post-treatment surface fuel loading following mastication. Plot-based destructive sampling was performed immediately following mastication at three stands in the Nez Perce Clearwater National Forest, Idaho, USA, to correlate post-treatment surface fuel loads and characteristics with pre-treatment LiDAR-derived forest metrics, speciﬁcally trees per hectare (TPH) and stand density index (SDI). Surface fuel loads measured in the stand post-treatment were consistent with those reported in previous studies. A signiﬁcant relationship was found between the pre-treatment SDI and total resultant fuel loading ( p = 0.0477), though not between TPH and fuel loading ( p = 0.0527). SDI may more accurately predict post-treatment fuel loads by accounting for both tree number per unit area and stem size, while trees per hectare alone does not account for variations of tree size and subsequent volume within a stand. Conditions within treated stands and fuels produced during mastication are highly variable and may explain the lack of relationship between fuel classes and loading. Root-mean-square errors of 36 and 46 percent of the random forest LiDAR models for SDI and TPH may limit the ability to predict the highly variable fuel loads produced from mastication. Use of LiDAR to predict fuel loading after mastication is a useful approach for managers to understand the e ﬃ cacy of fuel reduction treatments by providing information that may be helpful for determining areas where treatments can be most beneﬁcial.


Introduction
Due to the variability of species, management objectives, spatial configuration of management areas, regulatory restrictions, landowner funding availability, fuels characteristics, and other geographic and vegetative factors, developing a one-size-fits-all approach for wide-scale fire management is challenging [1][2][3][4]. However, forests with high fire risk must be actively managed [5]. Over the previous century, forest management practices such as fire exclusion have resulted in historically uncharacteristic stand attributes in many forests in the western United States, including dense, small-diameter stands with increased surface fuel loads [5,6]. Fuel reduction in stands that have lacked prior density management is a high priority in many areas of the western United States, especially on federal lands. Understanding the unique challenges and selecting strategies to best suit the needs of each management area, typically applied through one or more treatments applied at the stand level,

Study Site
The study sites were located in three stands in the Nez Perce-Clearwater National Forest, in north central Idaho, following fuel treatments to gather data pertaining to the resultant surface fuel loads. These treatments were part of the larger Orogrande timber sale and consisted of approximately 38 hectares (95 acres) of mechanical fuel treatment. The management units were predominantly mixed conifer forest type with slopes averaging 35% throughout the units. Mastication was successfully implemented in the three stands as a management alternative to timber harvest. The three stands treated for this study were originally planned for timber harvest. Harvesting was found to be financially infeasible due to the low value of harvested products and long-haul distances to the mill. The prescription developed for the project was intended to release remaining trees to increase timber value for future harvest while simultaneously decreasing stand density and increasing canopy base height to reduce the risk of crown fire using mastication. According to forest personnel, the management approach used in this project was its first application on the Nez Perce-Clearwater National Forest. The machine used to perform the mechanical fuel treatment was a Takeuchi TB290 compact excavator with a Fecon Bull Hog mastication head ( Figure 1). The machine weighs 8685 kg, is 2.2 m wide and 2.9 m long at the undercarriage, has a maximum reach of 7.4 m, and creates only 37.9 kPa of ground pressure when equipped with 450-mm-wide rubber tracks.

Study Site
The study sites were located in three stands in the Nez Perce-Clearwater National Forest, in north central Idaho, following fuel treatments to gather data pertaining to the resultant surface fuel loads. These treatments were part of the larger Orogrande timber sale and consisted of approximately 38 hectares (95 acres) of mechanical fuel treatment. The management units were predominantly mixed conifer forest type with slopes averaging 35% throughout the units. Mastication was successfully implemented in the three stands as a management alternative to timber harvest. The three stands treated for this study were originally planned for timber harvest.
Harvesting was found to be financially infeasible due to the low value of harvested products and long-haul distances to the mill. The prescription developed for the project was intended to release remaining trees to increase timber value for future harvest while simultaneously decreasing stand density and increasing canopy base height to reduce the risk of crown fire using mastication. According to forest personnel, the management approach used in this project was its first application on the Nez Perce-Clearwater National Forest. The machine used to perform the mechanical fuel treatment was a Takeuchi TB290 compact excavator with a Fecon Bull Hog mastication head ( Figure 1). The machine weighs 8685 kg, is 2.2 m wide and 2.9 m long at the undercarriage, has a maximum reach of 7.4 m, and creates only 37.9 kPa of ground pressure when equipped with 450-mm-wide rubber tracks. The treatment prescription for the stands included a target range for post-treatment stocking level. For the units in this study, the operator was instructed to leave 40 to 80 trees per hectare after treatment while removing only stems less than 18 cm (7 inches) in diameter. Further, all dead and down material up to 30 cm (12 inches) were masticated [39]. Post-treatment surface fuel sampling occurred in the masticated portions of three replicate stands: 117 (13 hectares), 120 (15 hectares), and 147 (10 hectares), within the management boundary. Due to many downed trees in stand 147, meeting the prescription specification for dead and downed material was not operationally feasible. Therefore, the mastication intensity for downed trees was reduced after stand 147 was partially treated. This prescription adjustment was used for treating the remainder of stand 147 and for the entirety of stands 117 and 120 [39] (Figure 2). The treatment prescription for the stands included a target range for post-treatment stocking level. For the units in this study, the operator was instructed to leave 40 to 80 trees per hectare after treatment while removing only stems less than 18 cm (7 inches) in diameter. Further, all dead and down material up to 30 cm (12 inches) were masticated [39]. Post-treatment surface fuel sampling occurred in the masticated portions of three replicate stands: 117 (13 hectares), 120 (15 hectares), and 147 (10 hectares), within the management boundary. Due to many downed trees in stand 147, meeting the prescription specification for dead and downed material was not operationally feasible. Therefore, the mastication intensity for downed trees was reduced after stand 147 was partially treated. This prescription adjustment was used for treating the remainder of stand 147 and for the entirety of stands 117 and 120 [39] (Figure 2).

LiDAR Processing and Sample Plot Selection
The Orogrande timber sale and the three stands (117, 120, 147) were within the extent of the 18,450-hectare (45,600-acre) Crooked River LiDAR acquisition flown in 2012 with a pulse density return of ≥4 points per square meter ( Figure 3). Field sampling inventory data from 91 20 × 20-m (1/10 acre) plots were run through the USFS Forest Vegetation Simulator [40] to summarize stand composition and structure. These forest inventory data were part of a previous sampling effort and were collected Sustainability 2020, 12, 7025 5 of 20 using field methods described in Falkowski et al. (2005) [22]. Random forest models [41,42] describing trees per hectare, total volume (m 3 ha −1 ), basal area (m 2 ha −1 ), and stand density index (SDI) were then developed using LiDAR metrics encompassing identical extents to the field sampling plots. These methods are consistent with those described in Becker et al. [43]. All random forest development and metric predictions were performed in the open source statistical analysis program, R, using the randomForest package [42,44].

LiDAR Processing and Sample Plot Selection
The Orogrande timber sale and the three stands (117, 120, 147) were within the extent of the 18,450-hectare (45,600-acre) Crooked River LiDAR acquisition flown in 2012 with a pulse density return of ≥4 points per square meter ( Figure 3). Field sampling inventory data from 91 20 × 20-m (1/10 acre) plots were run through the USFS Forest Vegetation Simulator [40] to summarize stand composition and structure. These forest inventory data were part of a previous sampling effort and were collected using field methods described in Falkowski et al. (2005) [22]. Random forest models [41,42] describing trees per hectare, total volume (m 3 ha −1 ), basal area (m 2 ha −1 ), and stand density index (SDI) were then developed using LiDAR metrics encompassing identical extents to the field sampling plots. These methods are consistent with those described in Becker et al. [43]. All random forest development and metric predictions were performed in the open source statistical analysis program, R, using the randomForest package [42,44].

LiDAR Processing and Sample Plot Selection
The Orogrande timber sale and the three stands (117, 120, 147) were within the extent of the 18,450-hectare (45,600-acre) Crooked River LiDAR acquisition flown in 2012 with a pulse density return of ≥4 points per square meter ( Figure 3). Field sampling inventory data from 91 20 × 20-m (1/10 acre) plots were run through the USFS Forest Vegetation Simulator [40] to summarize stand composition and structure. These forest inventory data were part of a previous sampling effort and were collected using field methods described in Falkowski et al. (2005) [22]. Random forest models [41,42] describing trees per hectare, total volume (m 3 ha −1 ), basal area (m 2 ha −1 ), and stand density index (SDI) were then developed using LiDAR metrics encompassing identical extents to the field sampling plots. These methods are consistent with those described in Becker et al. [43]. All random forest development and metric predictions were performed in the open source statistical analysis program, R, using the randomForest package [42,44].  Random forest ensemble learning algorithms were used because they provide excellent classification results, speed of processing, ability to reduce bias, and correlation and reduce overfitting compared to other classification and regression trees (CART) models, making them a widely used machine learning solution [45][46][47][48]. Random forest produces multiple decision trees using bagging and randomly selected subsets of training samples and variables to provide a majority vote from which a prediction is made. The default number of trees (n tree = 500) and the default number of variables split at each node (m try = square root of the total number of input variables) were used when building random forest models for each forest characteristic. 177 LiDAR metrics were available to develop each unique random forest algorithm, but only a subset of metrics were used in the final models, based on importance, defined using rfUtilities [49]. Random forest models were built using 2/3 of the data and validated using the remaining 1/3 of the data.
The entire Crooked River LiDAR acquisition was processed using the USDA LiDAR processing software FUSION version 3.60 to create an identical post-processed data structure as the initial 91 sampled plots [50]. This enabled the random forest models to be applied directly to the whole area to develop predicted metrics in 20 × 20-m pixels. Stand metrics derived from the LiDAR analysis included trees per hectare (TPH), total cubic foot volume (m 3 ha −1 ), total basal area (m 2 ha −1 ), and stand density index. SDI has been used in even-aged monocultures, and more recently in uneven-aged, mixed species stands to assess stand density as a function of quadratic mean diameter and stem density [51][52][53][54][55]. This metric was selected in addition to TPH to provide a more descriptive indication of stand density. The trees per hectare vector map was then stratified into four classes: 0-247; 248-494; 495-740; 741+. These classes were used to select sample plots within the study stands. Trees per hectare classes were used to stratify the selection of a broad distribution of relative stocking in sampled areas prior to mastication.

Field Sampling Procedures
Twelve plots were sampled within each of the three stands, with three representing each of the four levels of pre-treatment trees per hectare derived from the LiDAR data. The 20 × 20-m pixels chosen for sampling were randomly selected from all available pixels of the trees per hectare class within the stand boundaries. The resulting sampled pixels amounted to 36, with nine plots representing each of the four classes of pre-treatment trees per hectare. Trees per hectare was selected as the stand metric by which to select sample plots, due to the mastication treatment prescription being based on a goal trees per hectare post-treatment. All mastication treatments and sampling of fuel loading occurred during summer 2017.
Center points within the 20 × 20-m pixels were determined via ArcMap, and the resultant coordinates were used to locate the field plot centers. A simple method of plot center relocation was established to address situations where plot centers occur in areas that prohibited the sampling of fuel loading including tree stumps, roadways, rock outcroppings, and exposed mineral soil due to machinery movement. In these instances, plot centers were moved due north 3 m. If needed, plots were moved due west from original plot centers 3 m if the movement due north did not resolve the issue with the obstruction. A variation on destructive plot-based sampling was used to quantify fuel characteristics and fuel loading following mastication treatments within the three stands at each of the 36 plots [10,14,17].
Within each of the 36 sampled plots, fuel size classes were sampled in four quadrats. Once the plot center was located via GPS coordinates, 5-m vectors extending directly north, south, east, and west of the plot center were marked and established as the corner points for the quadrats. For instances in which uncharacteristic site conditions occurred within the quadrats, the frame was reflected over the transect. If this quadrat reflection did not resolve the issue and fuel collection within the quadrat was still not possible, the quadrat was excluded from sampling overall. Situations that would permit quadrat reflection over the transect or exclusion included buried logs, stumps, and rock outcroppings. In addition to the collection of fuels within the 25-cm squared quadrats, the fuel depth of masticated fuels was measured at two locations along the 5-m vectors (2.5 m and 5 m) and at the overall plot center ( Figure 4). To measure fuel depths, a cross-section of the forest floor was cleared using a trowel, and the depths were manually recorded. For the depths of the woody/masticated material, any branch or piece of woody debris above the measurement point was included in the depth measurement. Where site conditions prevented the measuring of fuel depths, the depth was measured 0.5 m from the original measurement point moving away from the plot center along the transect. If the depth was still not measurable at the second location, the measurement was omitted. Each of the 36 plots contained four separate fuel collection quadrats and nine fuel depth measurements.
Sustainability 2020, 12, x FOR PEER REVIEW 7 of 20 forest floor was cleared using a trowel, and the depths were manually recorded. For the depths of the woody/masticated material, any branch or piece of woody debris above the measurement point was included in the depth measurement. Where site conditions prevented the measuring of fuel depths, the depth was measured 0.5 m from the original measurement point moving away from the plot center along the transect. If the depth was still not measurable at the second location, the measurement was omitted. Each of the 36 plots contained four separate fuel collection quadrats and nine fuel depth measurements. Frames made of PVC pipe measuring 25 cm by 25 cm were built to establish the sampling area extent. Within the 25 × 25-cm collection quadrats, all fuel down to bare mineral soil was gathered, and any pieces extending outside the collection frame were trimmed using hand shears, and thus only pieces completely within the frame were collected ( Figure 5). All fuels collected were stored in paper bags labeled by their collection point and were brought back to the lab for detailed fuel composition analysis. Within the sampling quadrats, downed trees and logs were not sampled due to their irregular occurrence and collection difficulty for returning to the lab. Frames made of PVC pipe measuring 25 cm by 25 cm were built to establish the sampling area extent. Within the 25 × 25-cm collection quadrats, all fuel down to bare mineral soil was gathered, and any pieces extending outside the collection frame were trimmed using hand shears, and thus only pieces completely within the frame were collected ( Figure 5). All fuels collected were stored in paper bags labeled by their collection point and were brought back to the lab for detailed fuel composition analysis. Within the sampling quadrats, downed trees and logs were not sampled due to their irregular occurrence and collection difficulty for returning to the lab.
Sustainability 2020, 12, x FOR PEER REVIEW 7 of 20 forest floor was cleared using a trowel, and the depths were manually recorded. For the depths of the woody/masticated material, any branch or piece of woody debris above the measurement point was included in the depth measurement. Where site conditions prevented the measuring of fuel depths, the depth was measured 0.5 m from the original measurement point moving away from the plot center along the transect. If the depth was still not measurable at the second location, the measurement was omitted. Each of the 36 plots contained four separate fuel collection quadrats and nine fuel depth measurements. Frames made of PVC pipe measuring 25 cm by 25 cm were built to establish the sampling area extent. Within the 25 × 25-cm collection quadrats, all fuel down to bare mineral soil was gathered, and any pieces extending outside the collection frame were trimmed using hand shears, and thus only pieces completely within the frame were collected ( Figure 5). All fuels collected were stored in paper bags labeled by their collection point and were brought back to the lab for detailed fuel composition analysis. Within the sampling quadrats, downed trees and logs were not sampled due to their irregular occurrence and collection difficulty for returning to the lab. It was assumed that when locating sample plots in the field, there may be instances where mastication, though planned, does not occur. This was a result of inaccessibility, due to the steep or very uneven terrain where the operator chose not to treat the area for safety reasons. In these instances, three supplemental sample points for each classification level of trees per hectare were randomly selected in each of the three stands. If an originally designated sample plot was found to be within an un-masticated area, a randomly selected supplemental sample plot of the same TPH class was selected for sampling instead. These supplemental plots were randomly generated from the remaining available pixels not included in the initial stratification prior to field sampling, using the same method used in the initial plot selection. Supplemental plots were used once in unit 117, four times in unit 120, and four times in unit 147. The final sampled plots across stands 117, 120, and 147 are shown in Figure 6. It was assumed that when locating sample plots in the field, there may be instances where mastication, though planned, does not occur. This was a result of inaccessibility, due to the steep or very uneven terrain where the operator chose not to treat the area for safety reasons. In these instances, three supplemental sample points for each classification level of trees per hectare were randomly selected in each of the three stands. If an originally designated sample plot was found to be within an un-masticated area, a randomly selected supplemental sample plot of the same TPH class was selected for sampling instead. These supplemental plots were randomly generated from the remaining available pixels not included in the initial stratification prior to field sampling, using the same method used in the initial plot selection. Supplemental plots were used once in unit 117, four times in unit 120, and four times in unit 147. The final sampled plots across stands 117, 120, and 147 are shown in Figure 6.

Lab Measurements and Fuel Characterization
After field sampling was conducted, fuel collection bags were brought back to the lab to be processed by drying and sorting. A total of 133 quadrats of fuel was brought from the field for processing. If all plots had all quadrats collected, there would have been a total of 144. However, some quadrats were excluded from sampling because they were located on rock outcroppings, stumps, or other obstructions. One complete plot of quadrat fuel collections from stand 117, trees per hectare class 1, was misplaced during sampling, which constituted four of the eleven missing quadrat samples. Due to the omitted samples, stand 117 had two plots for trees per hectare class 1, resulting in 35 total plots rather than 36. Each collection sample was oven-dried at 105 degrees Celsius for 48 h, or until the sample weight stabilized, and was then weighted to the nearest gram. All fuels were then sorted, by quadrat, into five time-lag fuel classes: duff/litter and woody/masticated (1- (Figure 7). Sorted fuels were then individually weighed to the nearest gram to determine the proportion of overall mass that each fuel class represented. These proportions for each quadrat were averaged with corresponding plot quadrats to determine the fuel composition proportions by mass for the entire plot. For each of the 35 sample plots, the fuel bed volume was calculated by multiplying the average of the fuel depths at the nine measured locations within each plot by the dimensions of the collection frame. We then determined the bulk density of the fuels in each plot by dividing the average oven-dried weight of the fuel classes in the four collection quadrats by the corresponding volume. Plot level values were calculated using the averages of each quadrat within the plot for fuel loading (Mg ha −1 ) for the whole stand and by fuel class, fuel depth (cm), and bulk density (kg m −3 ).

Lab Measurements and Fuel Characterization
After field sampling was conducted, fuel collection bags were brought back to the lab to be processed by drying and sorting. A total of 133 quadrats of fuel was brought from the field for processing. If all plots had all quadrats collected, there would have been a total of 144. However, some quadrats were excluded from sampling because they were located on rock outcroppings, stumps, or other obstructions. One complete plot of quadrat fuel collections from stand 117, trees per hectare class 1, was misplaced during sampling, which constituted four of the eleven missing quadrat samples. Due to the omitted samples, stand 117 had two plots for trees per hectare class 1, resulting in 35 total plots rather than 36. Each collection sample was oven-dried at 105 degrees Celsius for 48 h, or until the sample weight stabilized, and was then weighted to the nearest gram. All fuels were then sorted, by quadrat, into five time-lag fuel classes: duff/litter and woody/masticated (1- (Figure 7). Sorted fuels were then individually weighed to the nearest gram to determine the proportion of overall mass that each fuel class represented. These proportions for each quadrat were averaged with corresponding plot quadrats to determine the fuel composition proportions by mass for the entire plot. For each of the 35 sample plots, the fuel bed volume was calculated by multiplying the average of the fuel depths at the nine measured locations within each plot by the dimensions of the collection frame. We then determined the bulk density of the fuels in each plot by dividing the average oven-dried weight of the fuel classes in the four collection quadrats by the corresponding volume. Plot level values were calculated using the averages of each quadrat within the plot for fuel loading (Mg ha −1 ) for the whole stand and by fuel class, fuel depth (cm), and bulk density (kg m −3 ). To assess the correlation between pre-treatment LiDAR-derived forest metrics and post-treatment fuel conditions, all statistical analyses were performed using the R statistical programming environment. Pearson correlation coefficients were calculated to evaluate the strength of association between predictors. Additionally, linear mixed effects models were used to model the relationship between fuel loading following mastication treatments and the trees per hectare, stand density index, and basal area of plots prior to mastication, using the nlme R package [57]. The general equation for mixed effect models is described as: where β are fixed effects, u are random effects, X is the model matrix for fixed effects, Z is the model matrix for random effects, ε is the vector of errors, R is the variance-covariance matrix of within-individual measurements, and D is the variance-covariance matrix of random effects [58]. In linear mixed effects, models evaluating predictors, stand, and TPH class were treated as random effects, with the TPH class nested within the stand. Random intercepts were used when fitting models. Due to the inconsistency of the mastication treatment in stand 147, two different mixed effects models were fit for each of the LiDAR predictors. One model contained all three stands, while the second model contained only stand 117 and 120. This was in order to avoid potential influential data artifacts associated with the treatment change in stand 147.

Results
Parameter estimates for the random forest models used in the pretreatment derivation of forest characteristics from LiDAR metrics for 20 × 20-m pixels are shown in Table 1. In this study, an acceptable maximum root-mean-square error (RMSE) of 50% of the prediction means was used based on values derived in previous studies [59,60]. The RMSE was within the acceptable range for the density (TPH), basal area (m 2 ha −1 ), and stand density index models. For the random forest model predicting stand volume (m 3 ha −1 ), an RMSE of 166.93 was about 54% of the predicted mean and just outside of the desired range. Predicted volume was therefore excluded from use in subsequent analyses. Model accuracies for forest metrics were 71.5%, 77.4%, 74.3%, and 79% for trees per hectare, basal area, volume, and stand density index, respectively, which are comparable to those obtained by Falkowski et al. [61] and Hudak et al. [60]. These pretreatment maps of predicted stand characteristics provided the basis for study plot selection and the subsequent regression modeling of post-treatment fuel loading. To assess the correlation between pre-treatment LiDAR-derived forest metrics and post-treatment fuel conditions, all statistical analyses were performed using the R statistical programming environment. Pearson correlation coefficients were calculated to evaluate the strength of association between predictors. Additionally, linear mixed effects models were used to model the relationship between fuel loading following mastication treatments and the trees per hectare, stand density index, and basal area of plots prior to mastication, using the nlme R package [57]. The general equation for mixed effect models is described as: where β are fixed effects, u are random effects, X is the model matrix for fixed effects, Z is the model matrix for random effects, ε is the vector of errors, R is the variance-covariance matrix of within-individual measurements, and D is the variance-covariance matrix of random effects [58].
In linear mixed effects, models evaluating predictors, stand, and TPH class were treated as random effects, with the TPH class nested within the stand. Random intercepts were used when fitting models. Due to the inconsistency of the mastication treatment in stand 147, two different mixed effects models were fit for each of the LiDAR predictors. One model contained all three stands, while the second model contained only stand 117 and 120. This was in order to avoid potential influential data artifacts associated with the treatment change in stand 147.

Results
Parameter estimates for the random forest models used in the pretreatment derivation of forest characteristics from LiDAR metrics for 20 × 20-m pixels are shown in Table 1. In this study, an acceptable maximum root-mean-square error (RMSE) of 50% of the prediction means was used based on values derived in previous studies [59,60]. The RMSE was within the acceptable range for the density (TPH), basal area (m 2 ha −1 ), and stand density index models. For the random forest model predicting stand volume (m 3 ha −1 ), an RMSE of 166.93 was about 54% of the predicted mean and just outside of the desired range. Predicted volume was therefore excluded from use in subsequent analyses. Model accuracies for forest metrics were 71.5%, 77.4%, 74.3%, and 79% for trees per hectare, basal area, volume, and stand density index, respectively, which are comparable to those obtained by Falkowski et al. [61] and Hudak et al. [60]. These pretreatment maps of predicted stand characteristics provided the basis for study plot selection and the subsequent regression modeling of post-treatment fuel loading.  Table 2 shows the summary data of the stands for the fuel collection as averages of the sampled plots and quadrats within each stand. Surface fuel loadings range from 9. , found by performing a Pearson's correlation test. Additionally, there was a positive (0.3577) and significant (p = 0.0349) correlation between SDI and fuel loading. No significant relationship was found between pre-treatment TPH and bulk density of resulting fuels (kg m −3 ); SDI and bulk density of resulting fuels; basal area (m 2 h −1 ) and loading or bulk density; nor between pre-treatment total volume (m 3 ha −1 ) and fuel loading or bulk density (Table 3). Based on the results of the correlation tests, the linear mixed effects model was fitted to evaluate the relationship of pre-treatment TPH and the resulting fuel loads as well as SDI and the resulting fuel loads. The linear mixed effects models predicting the total fuel loading of all time-lag classes from pre-treatment TPH showed no significant relationship between the two factors (test statistic = 0.05318, df = 22, p-value = 0.0527) when accounting for all three stands. In the reduced model, pre-treatment trees per hectare and fuel loading were significant (p-value = 0.0066) ( Table 4). The linear mixed effects models predicting the total fuel loading of all time-lag classes from pre-treatment SDI showed significant relationships between the two factors for all three units (p = 0.0477) and when assessing stand 117 and 120 alone (p = 0.0337) ( Table 4). Figure 8 shows the associated relationships between SDI and the resulting fuel load for all three units. The black line represents the regression line of the complete data set and the individual regression line for each stand. 0.05318, df = 22, p-value = 0.0527) when accounting for all three stands. In the reduced model, pre-treatment trees per hectare and fuel loading were significant (p-value = 0.0066) ( Table 4). The linear mixed effects models predicting the total fuel loading of all time-lag classes from pre-treatment SDI showed significant relationships between the two factors for all three units (p = 0.0477) and when assessing stand 117 and 120 alone (p = 0.0337) ( Table 4). Figure 8 shows the associated relationships between SDI and the resulting fuel load for all three units. The black line represents the regression line of the complete data set and the individual regression line for each stand.  Additional linear mixed effects models were developed to further assess the impact of pre-treatment TPH and SDI on resulting fuels loads for each of the 5 time-lag fuel classes (litter/duff, 1-h, 10-h, 100-h, 1000-h). Only the litter/duff fuel class loading was found to be significantly correlated to pre-treatment trees per hectare for all three stands ( Table 5). The litter/duff fuel class is generally independent of the mastication process, as most fuels in this class were present before treatment. However, when assessing stands 117 and 120, litter/duff (p = 0.0242), 1-h (p = 0.0232), and 100-h (p = 0.0059) were found to be significant. When assessing the relationship between SDI and the resulting fuel loading across all five time-lag cases, the model containing all three units showed that both litter/duff (p = 0.0042) and 100-h (p = 0.0293) were significant, while the model describing units Additional linear mixed effects models were developed to further assess the impact of pre-treatment TPH and SDI on resulting fuels loads for each of the 5 time-lag fuel classes (litter/duff, 1-h, 10-h, 100-h, 1000-h). Only the litter/duff fuel class loading was found to be significantly correlated to pre-treatment trees per hectare for all three stands ( Table 5). The litter/duff fuel class is generally independent of the mastication process, as most fuels in this class were present before treatment. However, when assessing stands 117 and 120, litter/duff (p = 0.0242), 1-h (p = 0.0232), and 100-h (p = 0.0059) were found to be significant. When assessing the relationship between SDI and the resulting fuel loading across all five time-lag cases, the model containing all three units showed that both litter/duff (p = 0.0042) and 100-h (p = 0.0293) were significant, while the model describing units 117 and 120 showed that only 100-h (p = 0.0476) was significant. The data for the SDI model are shown in Figure 9. 117 and 120 showed that only 100-h (p = 0.0476) was significant. The data for the SDI model are shown in Figure 9.   Our second research objective explored whether the in-plot distribution of masticated fuels among the five time-lag fuel classes was impacted by the stand density of the plot prior to treatment. Mixed effects models were developed to assess these questions, with results being found in Table 6. 10-h fuels in the sampled plots were the only fuel class found to change significantly as the pre-treatment TPH changed (test statistic = −0.01581, p = 0.0178, df = 22). SDI was shown to have a significant relationship with post-mastication fuel loading for both the 10-h (p = 0.0302) and 100-h (p = 0.0406) fuel classes (Table 6). Table 6. Mixed effects model summary assessing pre-treatment stand density (TPH) and stand density index (SDI) impact on the percentage of total fuel load (Mg ha −1 ) by time-lag fuel class. The data distribution for the five time-lag classes as a percentage across the range of stand density indices and fuel loading are shown in Figure 10. Woody and mulched fuels (1-h, 10-h, 100-h, and 1000-h) were found to be between 30-85% of the overall fuel loads across all plots, which shows the variability found within the stands. 1-h, 10-h, 100-h, and 1000-h fuels contained 6-27%, 36-94%, 0-57%, and 0-13% of the woody and mulched fuels, respectively. On average, woody fuels made up about 58.2% of the overall surface fuel loading across all sites, with 1-h, 10-h, 100-h, and 1000-h fuel classes averaging 7.5%, 33.2%, 17.3%, and 0.2% of the loading, respectively.
Our second research objective explored whether the in-plot distribution of masticated fuels among the five time-lag fuel classes was impacted by the stand density of the plot prior to treatment. Mixed effects models were developed to assess these questions, with results being found in Table 6. 10-h fuels in the sampled plots were the only fuel class found to change significantly as the pre-treatment TPH changed (test statistic = −0.01581, p = 0.0178, df = 22). SDI was shown to have a significant relationship with post-mastication fuel loading for both the 10-h (p = 0.0302) and 100-h (p = 0.0406) fuel classes (Table 6). The data distribution for the five time-lag classes as a percentage across the range of stand density indices and fuel loading are shown in Figure 10. Woody and mulched fuels (1-h, 10-h, 100-h, and 1000-h) were found to be between 30-85% of the overall fuel loads across all plots, which shows the variability found within the stands. 1-h, 10-h, 100-h, and 1000-h fuels contained 6-27%, 36-94%, 0-57%, and 0-13% of the woody and mulched fuels, respectively. On average, woody fuels made up about 58.2% of the overall surface fuel loading across all sites, with 1-h, 10-h, 100-h, and 1000-h fuel classes averaging 7.5%, 33.2%, 17.3%, and 0.2% of the loading, respectively.

Discussion
Total surface fuel loadings varied widely across our study plots (26.4 to 158.2 Mg ha −1 ), with woody/masticated surface fuels representing a similarly wide range (7.7 to 127.5 Mg ha −1 ). Across all stands and plots, total surface fuel loading averaged 89.8 Mg ha −1 , and woody surface fuels averaged 53.9 Mg ha −1 . These total surface fuel loads were similar to those reported by Stephens and Moghaddas [16] in the Sierra Nevada Mountains (87.1 and 93.8 Mg ha −1 ), Reiner et al. [18] in the Sierra Nevada (81.6 Mg ha −1 ), Kane et al. [14] at study sites in Northern California and southwestern Oregon (83.6, 83.8 and 71.1 Mg ha −1 ), and Hood and Wu [17] in the Northern Rockies (82.0-95.9 Mg ha −1 ), but were higher than those reported by Brewer et al. [62] in mixed conifer Idaho stands (58.4 Mg ha −1 ) and lower than those reported by Battaglia et al. [10] in mixed conifer stands in Colorado (110.4 Mg ha −1 ). This finding was not surprising, as masticated fuel beds and characteristics have a wide variability across sites and regions. The attempt to develop the surface fuel prediction model at a 20 × 20-m resolution in this study was meant to help address this site variability.
Based on results reported by previous mastication studies, the average total fuel depth of 16.7 cm recorded in the mixed conifer stands we studied was significantly higher than those in Battaglia et al. [10], but only slightly higher than those shown by Stephens and Moghaddas [16]: 14.6 and 14.7 cm. All fuel sampling in this study occurred within a month of mastication treatment, so fuels were still green at the time of collection and had not settled to the forest floor, whereas sampling of masticated fuels for other studies occurred 2-6 years post-mastication [10,14]. These temporal changes in masticated fuel beds make the generalization of loadings difficult, especially across broad geographic extents. For example, a recently masticated stand may indicate greater fuel depths than a stand masticated several years ago, due to the decomposition and deterioration of fuel structural integrity. The additional compaction of fuel beds over time as they settle may affect subsequent fire behavior. Therefore, the development of the model estimating surface fuel characteristics directly after masticating would provide a consistent expectation of fuel loads, as was done in this study.

Relationship between Pre-Treatment Stand Characteristics and Fuel Loading
Through the analysis performed across the three stands we studied-117, 120, and 147no significant relationship was found between overall fuel loading following mastication and the pre-treatment tree per hectare we derived from LiDAR. However, SDI was found to be a significant predictor variable for post-mastication fuel loading when accounting for all management units and time-lag classes jointly. Initially, we expected to see an increase in the fuel loading as the pre-treatment TPH increased. It was believed that, given a consistent prescription implementation, greater TPH would result in more fuel, as a greater number of standing trees were mulched to meet treatment objectives. The contrary findings for absolute stand density may have resulted from the inconsistency in the initial treatment of stand 147, which was then corrected. When excluding stand 147, increasing pre-treatment TPH resulted in greater loading, as expected. However, even when including data from stand 147, the p-value of 0.0526 was just outside the level of significance needed to reject the null. The significance of both mixed effects models for SDI (p = 0.0477 and p = 0.0337) indicates fuel loading may be more accurately predicted using a metric that accounts for both tree size and number, as opposed to simply using trees per hectare where only the number of stems is accounted for.
Stand 147 was the first stand treated and was initially treated to prescription specifications. It was found that once treatment began, the original degree to which large downed woody debris was to be treated was operationally infeasible due to the increased treatment time. Further, stand 147 contained a small pocket of lodgepole pine killed by beetle, with a significant portion of downed trees which were, under the original treatment specs, to be masticated heavily. This resulted in a larger amount of masticated fuels in plots with relatively low stand densities.
The decision to retain stand 147 in the analysis was made to maximize the data available for assessment and provide a realistic portrayal of the large variability of mastication treatments. As a relatively new treatment option being deployed over large areas, it is likely that similar inconsistency in operational treatments may occur during implementation and administration of mechanical fuel treatments, particularly as operators familiarize themselves with prescription requirements in new treatment areas. However, when the treatment prescription and execution was consistent for the entire stand, as seen in stands 117 and 120, a clear relationship between pre-treatment trees per hectare and fuel loading was seen. Given the potential for variability of mastication treatments and the heterogeneity of stand conditions in practice, the predictive success of SDI in a "real-world" management scenario is valuable for future modeling efforts.
Given the extent of our results, it remains unclear if a relationship exists between the fuel loading following mastication treatments and the pre-treatment stand density based solely on TPH or stand basal area, but SDI is a useful predictor. In mastication, the conservation of mass must be considered, as materials are not removed from the stand after treatment but rearranged in different physical forms. A larger masticated tree will understandably produce a larger amount of masticated material than a tree of smaller size. For example, two stands may both have similar numbers of trees per unit area, but one stand may have a larger average tree diameter than the other. If both stands are treated to the same prescription and reduced to a defined tree per unit area, it would be expected that the stand with the larger average stand diameter would produce heavier masticated fuel loadings. Accounting for both stem number per unit area and tree size in a single pre-treatment stand metric, SDI addresses this issue. Therefore, alternative approaches to modeling landscape scale fuel loading following mastication based on pre-treatment stand conditions that incorporate both stem numbers and size may offer improved prediction in future research and should be the focus of future study design and implementation.
When assessing the fuel loading for each time-lag class (litter/duff, 1-h, 10-h, 100-h, 1000-h), the litter duff class showed a significant relationship for the TPH and SDI models (Table 5). This may be a result of greater stand density, leading to higher amounts of organic material and litter on the forest floor. In all, minimal 1000-h fuels were collected at the plots, limiting the available data for the particular classes and making predictions difficult. This finding corroborates Kane et al.'s [11], who found that the plot-based method of surface fuel sampling does not assess a large enough area to effectively capture the presence of 1000-h fuels as well as planar intersect methods. This is a result of 1000-h fuels generally occurring less frequently than other fuel classes in fuel beds.

Relationship between Pre-Treatment Stand Density and Fuel Class Distribution
As shown above (Table 6), only the 10-h fuels expressed as a percentage of the overall surface fuel loading were found to change as the TPH increased. It is unclear why the percentage of 10-h fuels would decrease with increasing TPH, but this may be a result of changes in treatment implementation. For example, the operator may spend less time masticating trees to maintain production in a denser stand, resulting in an increase in the proportion of larger fuel classes. It would be expected that, with one fuel class decreasing over increasing TPH, another fuel class would increase. This was seen in the SDI model, where the significant decrease in 10-h fuels (coefficient = −0.04616, p = 0.0302) was matched by a significant increase in 100-h fuels (coefficient = 0.042089, p = 0.0406). With increasing stand density, it is possible the operator attempted to maintain the desired operational production by decreasing the time spent masticating each tree. As a result, trees would be masticated less thoroughly, and there would be a larger percentage of larger fuel particles. The 10-h fuel class accounted for the highest fuel loads across all classes by a considerable amount in our study, which is consistent with other studies [10,14]. Given the variability of the mastication as a whole, and the wide range of fuel loadings across 1-h, 10-h, 100-h, and 1000-h fuel classes, our results show that the distribution of surface fuels among time-lag fuel classes was not clearly modeled as a function of changes in TPH and SDI alone, apart from the 10-h and 100-h classes, in the case of this study.

Study Limitations and Future Work
Mastication is a highly variable operation impacted by many factors, and it is understood that there are some limitations to the scope of our research that should be addressed in future studies. One factor to consider in future applications of this methodology is the pixel size at which the LiDAR metrics were predicted. In the study development, it was believed that maximizing prediction resolution was the best option. Mastication, however, is a variable process, resulting in a scattered distribution of fuels on the forest floor. The directionality, travel distance, and particle size of comminuted materials may be affected by the type of mastication head (disk vs. drum), equipment type (all surface vehicle vs. excavator carrier), equipment horsepower, boom or attachment height, local topography within the stand, or other factors. While the sampling method developed by Hood and Wu [17] attempts to address the variability in stand and site conditions by sampling across multiple quadrats within the same plot, what was not accounted for in our study design was the possibility of fuels from adjoining pixels being distributed inter-pixel. During the observation of the treatment, fuels were clearly distributed more irregularly and further than anticipated. The 20 m × 20-m pixels used in the plot selection may have been too small to limit the influence of surrounding pixels in the resulting fuels found during sampling. For instance, a stand with a high stand density may have resulted in fuels initially in the stand as standing trees being distributed to an adjacent stand of a lower stand density, or to areas within the same stand that were not accounted for in our sampling design. During sampling, it would then appear that the pixel with the lower stand density was responsible for creating greater fuel loads than was possible. By decreasing the resolution and increasing the pixel size, this may be avoided.
In the fuel collection process, future studies should incorporate a hybrid, plot-based, and planar intersect method, as suggested by Kane et al. [14]. Doing so may help to ensure a more accurate representation of the fuel classes, as 1-h and 10-h fuels are more accurately represented in plot-based sampling [14], while planar intersect methods cover greater proportions of the overall masticated area, properly representing the 100-h and 1000-h fuels that may be missed in plot-based approaches [11]. Supplemental planar intersect sampling was not performed in this study due to the small 20 m × 20-m pixel size and the concern that sufficiently long intersect paths would extend too far to sample plot edges and be impacted by the distribution of fuels from adjacent pixels. Increasing the pixel size used in predicting the stand characteristics from the LiDAR, as described above, would enable an easier implementation of supplemental planar intersect sampling. Additionally, due to the variability of fuel distribution across the forest floor, plot-based sampling in future studies should use larger sampling quadrats than the 25 × 25-cm ones used in this study. Alternatively, a larger number of 25 × 25-cm sampling quadrats may also provide a greater representation of overall fuel variability within the sampling plot. The goal in using a smaller sampling quadrat in this study than those described in previous studies [17,63] was to create an efficient and effective sampling procedure. However, larger quadrats will provide a greater representation of overall surface fuel loadings and should be studied.

Conclusions
The ability to quickly, efficiently, and effectively predict surface fuel loads resulting from mastication treatments is a valuable tool, as increased implementation of this management technique occurs. A variety of research and management questions regarding the longevity, fuel bed characteristics, and fire behavior within masticated fuels exist and will increase in relevance as LiDAR data become more widespread, along with the use of mastication to reduce fuels in stands where commercial thinning may be infeasible or more difficult to implement administratively. Existing methods for predicting surface fuel loads rely on intensive, time-consuming sampling following treatment. While existing methods are effective for estimating fuel loading, methods based on remote sensing may help managers to proactively plan and predict post-treatment fire behavior over large areas to optimize treatments in ways that incorporate topography and stand adjacency.
The results from this study showed that pre-treatment stand density metrics that account only for tree number per unit area, such as TPH, were not good predictors of resulting surface fuel loads following mechanical fuel treatments with the sampling design and sample size we evaluated. TPH prior to treatment was not directly related to the distribution of fuel time-lag classes within the fuel bed, although the percentage of 10-h fuels could be predicted from pre-treatment conditions. However, stand density index, which accounts for both the relative stem number and DBH of the stand, effectively predicts post-treatment fuel loading across the whole study area. Further, SDI predicted that as the density of a stand increases, a greater percentage of the overall fuel load consisted of 100-h fuels, while 10-h fuels decreased in percentage, likely a result of operational adjustments. Future modeling efforts to predict post-mastication fuel loading should account for both the stem number and stem size, as stand density alone may not provide the necessary predictive ability. Attempting to predict resulting fuels from the number of trees per unit area alone does not account for variable volumes of materials in trees of different diameters. Stand density measures, such as SDI, provide greater insight into stand composition and overall stand biomass, which is significant when predicting fuel load volumes resulting from the physical conversion of standing biomass to surface-based mulched materials. Two stands with identical TPH may contain varying amounts of biomass as standing trees, whereas it is expected that two stands with identical SDI would have the same amount of overall biomass given similar forest types and species.
We believe that revisiting these methods, while taking into account the sampling considerations mentioned in the discussion, is an important undertaking and could lead to the increased implementation and effectiveness of mastication treatments. The rapid onset of LiDAR-derived models to map individual-tree locations and stem characteristics, coupled with onboard GNSS mapping of spatially, explicit, real-time equipment activities, offer the promise of improved high-resolution fuel bed prediction in the immediate future. Further expanding the scope of the field sampling to multiple, unique forest types, operators, and prescriptions would better capture the variability associated with the masticated surface fuel loads. Future work should address these factors more closely, though the determination of their impacts will likely require sampling at a higher intensity than that performed in this study, or with a sampling design that directly accounts for the spatial resolution at which comminuted material is scattered as a function of localized stand density, treatment prescription, topography, equipment type and size, and the pattern of equipment movements.