A Comparison of Four Spatial Interpolation Methods for Modeling Fine-Scale Surface Fuel Load in a Mixed Conifer Forest with Complex Terrain

: Patterns of spatial heterogeneity in forests and other ﬁre-prone ecosystems are increasingly recognized as critical for predicting ﬁre behavior and subsequent ﬁre effects. Given the difﬁculty in sampling continuous spatial patterns across scales, statistical approaches are common to scale from plot to landscapes. This study compared the performance of four spatial interpolation methods (SIM) for mapping ﬁne-scale fuel loads: classiﬁcation (CL), multiple linear regression (LR), ordinary kriging (OK), and regression kriging (RK). These methods represent commonly used SIMs and demonstrate a diversity of non-geostatistical, geostatistical, and hybrid approaches. Models were developed for a 17.6-hectare site using a combination of metrics derived from spatially mapped trees, surface fuels sampled with an intensive network of photoload plots, and topographic variables. The results of this comparison indicate that all estimates produced unbiased spatial predictions. Regression kriging outperformed the other approaches that either relied solely on interpolation from point observations or regression-based approaches using auxiliary information for developing ﬁne-scale surface fuel maps. While our analysis found that surface fuel loading was correlated with species composition, forest structure, and topography, the relationships were relatively weak, indicating that other variables and spatial interactions could signiﬁcantly improve surface fuel mapping.


Introduction
Forest fuel inventory and monitoring provide the basis for fuel management activities, including assessing wildfire hazards and risk, prescribed fire planning, designing silvicultural treatments, and predicting fire behavior and effects at various scales. The most commonly assessed physical fuel attribute is the load (kg m −2 ). Fuel load is a required input to nearly all fire behavior and effects models and is coupled to terrestrial carbon inventories and wildlife habitat assessments [1][2][3][4]. Fuel inventory approaches have traditionally assumed that spatial variability in fuel load is of little consequence for management decisions and thus focus on providing estimates of the spatially averaged fuel load for a given area based on a limited set of sampled locations. Yet, recent studies highlight that fine-scale variability in the fuel complex, as exists in virtually all wildland fuel beds, exerts considerable influence on many ecologically relevant fire behavior and effects metrics [5][6][7][8][9][10][11][12]. However, directly mapping fine-scale fuel variability is costly and time-consuming, especially across large areas [11,13]. To overcome these limitations and generate spatially continuous fuel load maps to support forest and fire management activities, spatial interpolation methods (SIMs) can be used to estimate the fuel load at unsampled locations.
Developing continuous surface fuel maps using spatial interpolation requires three sequential steps. The first step involves collecting fuel inventory data at a subset of points within an area of interest. The next step involves fitting a model of the fuel load using a selected SIM, and the final step involves developing a map from the fitted model. Two critical decisions in this process influence the map accuracy, namely creating a sampling scheme at a spatial sampling frequency sufficient to capture spatial autocorrelation and selecting an appropriate SIM. Although decisions in both stages are essential, this manuscript focuses on the effect of SIM selection on map accuracy.
The general approach of SIMs is to use weighted averaging to estimate a value of interest at an unmeasured location,ẑ(x 0 ), based on data from sample points, z(x i ), and a weighting factor, λ i, as per Equation (1): The major difference among SIMs stems from using different approaches to calculate the weighting function, such as including spatial autocorrelation and/or auxiliary explanatory variables [14]. Methods incorporating spatial autocorrelation include geostatistical techniques, such as ordinary kriging, and non-geostatistical methods, such as inverse distance weighting. These approaches exploit the correlation between sample variability and the distances between samples to predict values at unmeasured locations [15]. Alternatively, approaches such as regression modeling and classification rely on a relationship between the primary and auxiliary variables such as terrain, forest type, or vegetation characteristics to predict values at unsampled locations [4,16]. Though these approaches exploit relationships between the variable of interest and covariates, they do not include any information regarding the spatial structure of the data. Hybrid methods such as regression kriging that utilize both correlations with auxiliary data and spatial autocorrelation have gained popularity for spatial interpolation with environmental data [17][18][19]. Although several researchers have compared the accuracy of various SIMs for mapping environmental data [19,20], no universally optimal method for interpolation exists. The selection of the most accurate SIM depends on how well the method reflects the processes governing the spatial distribution of the variable of interest, including the variable's intrinsic spatial autocorrelation and correlation with any auxiliary variables [21,22].
One of the most significant challenges in developing accurate surface fuel maps is capturing the spatial variability in the fuel load within and between different fuel components. This variability arises due to interactions between the physical environment (e.g., climate, soils, and topography) and ecological processes (e.g., productivity, deposition, decomposition, and disturbances) that determine surface fuel accumulation [23][24][25]. One of the most commonly used approaches to capture the variability in surface fuels is to classify an area into unique groups using auxiliary, often remotely sensed, data (e.g., vegetation type, topographic data, or land-use classes) and then assign a fuel load to all areas of a given category based on the sampled data (e.g., 4). A drawback of a classification approach is that it reduces variability in the data to a few unique values, often with unrealistically stark class boundaries. A national-scale example of this is the fuel characteristics classification system (FCCS) maps of fuel load produced by LANDFIRE [26] at a 30 m spatial resolution. Although Prichard et al. [27] applied an aspatial method to consider the inherent variability and associated uncertainties in the fuel components that comprise wildland fuel beds, their intention was to inform larger-scale emissions estimates. For physical fire behavior modelling at finer scales, however, spatially explicit estimates of component fuel loads are required. Other researchers have used regression-based methods where the relationship between auxiliary variables and the surface fuel load is modeled and used to predict values at unmeasured locations. Several studies have found significant correlations between the surface fuel load and topographic and forest structural metrics [11,28,29]. However, these studies often report weak correlations, indicating that surface fuel maps developed based on these relationships may have limited accuracy. In other cases, studies have found no discernable correlation between surface fuel load and overstory forest structure (e.g., [30]), suggesting that including these variables in the development of surface fuel maps will not improve map accuracy. Other studies have found that surface fuels exhibit strong fine-scale (1-20 m) spatial autocorrelation, which, if taken into account, may increase fuel map accuracy [10,25,[31][32][33]. Although previous research has acknowledged spatial autocorrelation in fuel loads (e.g., 11,16), only Pierce et al. [33] have explicitly assessed the degree to which this improves spatial interpolation. Their results indicate that including spatial autocorrelation did not significantly improve the fine-scale predictive accuracy compared to linear regression approaches. However, the scale of the analysis (30 m) utilized in their study was greater than the spatial autocorrelation scale of surface fuels [10,32], which may have limited improvements in local scale predictive accuracy. More generally, we posit that the differences in the spatial scale of the sample data used in these studies might explain the seemingly disparate results and conclusions.
In this paper, we compare the performance of four SIMs for estimating and mapping fine-scale fuel load (1 m × 1 m resolution) in a mixed conifer forest in Colorado, USA. The four approaches were: classification, multiple linear regression, ordinary kriging, and regression kriging. We chose these methods because they are commonly used in ecological studies and cover a range of SIMs, including non-geostatistical, geostatistical, and hybrid approaches.

Study Area and Data Collection
We conducted this study on the 17.  Figure 1B), a dry continental climate with 660.7 mm of rain per year, and a mean daily temperature ranging from −4.7 • C in January to 14.0 • C in August (prism.oregonstate.edu; accessed on 6 June 2022). Topographically, the PFDP is shaped by two significant ridges, one oriented west-east in the northern portion of the plot and another oriented northwest-southeast in the southwestern portion ( Figure 1B). The dominant overstory vegetation for the study site included ponderosa pine (Pinus ponderosa Lawson and C. Lawson) and quaking aspen (Populus tremuloides Michx.) on the southern aspects and mixtures of Engelmann spruce (Picea engelmannii Parry ex Engelm.), blue spruce (Picea pungens Engelm.), and Douglasfir (Pseudotsuga menziesii (Mirb.) Franco) on the northern aspects. The site's average density, basal area, and quadratic mean diameter were 1040 trees ha −1 , 29.5 m 2 ha −1 , and 19.0 cm, respectively.
All trees in PFDP at least 1.37 m tall were spatially mapped and had their species, diameter at breast height, height, and crown base height recorded. Trees were mapped using tapes within a surveyed grid of 20 m × 20 m cells, and independent validation of the tree locations showed a mean error of less than 0.1 m. To characterize the surface fuel load across the site, we visually estimated the 1, 10, and 100 h dead, down, and woody fuel load on 1 m 2 irregularly located plots (n = 429) using the photoload method ( Figure 1C; [34]). The irregular design resulted in sample spacing ranging from 1 to 580 m covering the reported spatial autocorrelation scales for western U.S. forests [10,32]. Because photoload estimates are often consistent but biased, we randomly chose 50 additional 1 m 2 plots adjacent to the PFDP for the double sampling correction of the bias [35]. On each of the 50 plots, we first visually estimated the fuel load for each time-lag size class and then extracted the fuel for oven-drying following the method by Matthews [36]. We then developed a correction equation for each fuel component by linearly regressing the collected fuel load against the visually estimated fuel load [35]. We estimated the litter and duff fuel load using the depth-to-load method [37]. We first calculated means for the litter and duff depth in each plot from nine randomly selected points and then multiplied the average depth by a locally derived bulk density estimate. We estimated the locally derived bulk density from 20 randomly selected off-site 0.09 m 2 plots. On each plot, nine sample pins were inserted to the top of the litter layer, we extracted all the litter to estimate the dry weight, multiplied the average height of the exposed pins by the sample area to calculate the volume, and divided the dry weight by the volume to determine the bulk density with the same process performed for the duff layer. Finally, to calculate each plot's total fuel load, we summed the individual components' fuel load (1, 10, and 100 h dead, down, and woody fuel, and litter and duff). All trees in PFDP at least 1.37 m tall were spatially mapped and had their species, diameter at breast height, height, and crown base height recorded. Trees were mapped using tapes within a surveyed grid of 20 m × 20 m cells, and independent validation of the tree locations showed a mean error of less than 0.1 m. To characterize the surface fuel load across the site, we visually estimated the 1, 10, and 100 h dead, down, and woody fuel load on 1 m 2 irregularly located plots (n = 429) using the photoload method ( Figure 1C; [34]). The irregular design resulted in sample spacing ranging from 1 to 580 m covering the reported spatial autocorrelation scales for western U.S. forests [10,32]. Because photoload estimates are often consistent but biased, we randomly chose 50 additional 1 m 2 plots adjacent to the PFDP for the double sampling correction of the bias [35]. On each of the 50 plots, we first visually estimated the fuel load for each time-lag size class and then extracted the fuel for oven-drying following the method by Ma hews [36]. We then developed a correction equation for each fuel component by linearly regressing the collected fuel load against the visually estimated fuel load [35]. We estimated the li er and duff fuel load using the depth-to-load method [37]. We first calculated means for the li er and duff depth in each plot from nine randomly selected points and then multiplied the average depth by a locally derived bulk density estimate. We estimated the locally derived bulk Using the stem-mapped data and a 10 m digital elevation model (DEM; available at nationalmap.gov; accessed on 6 June 2021), we estimated several forest structural auxiliary variables across multiple scales ranging from 1 to 20 square meters using a rectilinear grid around the center of each sampling point. Auxiliary variables associated with forest structure included cover type, basal area, distance to nearest tree, distance to nearest large tree, trees per hectare (TPH), and TPH of large trees. For this analysis, we considered all trees with a diameter at breast height greater than 12.5 cm to be a large tree. Cover type was classified by applying a Gaussian kernel smoother [38] to the stem-mapped data to estimate the probability of occurrence by species and then assigning the cover type based on the species with the highest probability of occurrence at each location. Topographic variables were calculated from the DEM including slope, aspect, and flow direction [39], topographic wetness [40], and curvature [41].

Spatial Interpolation Methods
We compared the performance of four SIMs for mapping fine-scale (1 m 2 ) fuel loads: classification (CL), multiple linear regression (LR), ordinary kriging (OK), and regression kriging (RK). These methods represent commonly used SIMs and demonstrate a diversity of non-geostatistical, geostatistical, and hybrid approaches. We assessed the SIMs for the total surface fuel load, litter, duff, and 1, 10, and 100 h dead, down, and woody fuels as these are standard inputs to fire behavior and effects models [10]. A description of each SIM is provided below.

Classification
Classification is an aspatial SIM that utilizes the relationships between a categorical factor and sampled data. Classification approaches consist of two main steps. First, secondary information is used to identify several classes, k, within the region of interest such that all points exclusively belong to one and only one class. For each class, a value is estimated from the sample data and assigned to each point following Equation (2): whereẑ ik is the predicted value at location x ik belonging to class k, µ is the global mean of samples, α k is the difference between µ and the mean of class k, and ε ik is a random error term. In keeping with the general formula presented in Equation (1), the weighting function for classification approaches is estimated per Equation (3), where n k represents the number of observations in class k: The use of classified means as an interpolation method assumes that in the absence of direct observation, the best prediction at a location is the within-class mean. CL approaches assume that the within-class values are random and independently distributed, are not spatially autocorrelated, and that the variance is homogeneous. A more thorough discussion of CL approaches in spatial interpolation can be found in the study by Webster and Oliver [42].
As described above, we developed cover type maps for scales from 1 to 20 square meters. Our final CL model for each fuel component was based on the scale that provided the lowest Akaike information criterion (AIC).

Multiple Linear Regression
Multiple linear regression is an aspatial interpolation approach that models the relationship between auxiliary information and the primary variable of interest to predict values at unsampled locations. This relationship is frequently estimated following the classic ordinary least-squares approach, as per Equation (4): where x i are the auxiliary explanatory variables, β i is the linear slope coefficient that corresponds to each x i , α is the intercept, and ε is the residual error term. When using LR approaches, the weighting function, λ i , is set to 1 such that the prediction of the variable of interest at any point is purely a function of the auxiliary information at that location. Estimates based on LR assume that the data are random, not spatially autocorrelated, and have homogeneous variance. Furthermore, LR assumes a linear relationship between primary and auxiliary data and nominal multicollinearity among the auxiliary data. Additional details of LR can be found in the study by Eberly [43].
We modeled each fuel component as a function of the forest structure and topographic auxiliary variables for a single scale using a backward model selection process to remove non-significant (α > 0.05) predictor variables. We evaluated normality and multicollinearity through a visual inspection of the Q-Q plot and VIF, respectively. For each fuel component, we selected the regression model that resulted in the lowest AIC for further assessment.

Ordinary Kriging
Ordinary kriging is one of the most common geostatistical spatial interpolation methods in the environmental sciences. In OK, values at unsampled locations are estimated based on distance-dependent spatial autocorrelation modeled from observed data. The spatial autocorrelation term is identified through a sample semivariogram, which estimates the semivariance, λ(h), as a function of the lag distance, h, the number of data pairs at a given lag distance, N(h), and the measured value of the primary variable of interest, z, at locations x i and x i + h, as per Equation (5): A statistical model is fit to the sample semivariogram to estimate the coefficients for kriging. Three coefficients are typically estimated based on the fitted model: the range, sill, and nugget. The sill represents the maximum variation in the data, the range parameter is the lag distance at which the sill is reached and provides an estimate of the limit of spatial dependence, and the nugget represents the sum of the unresolved error occurring in the measurement or autocorrelation at a scale below the lag distances sampled [14]. In this study, we estimated the experimental variance by fitting a Matern function [44] to our sample variogram using an ordinary least-squares approach, following Equation (6): where c is the semivariance or sill, α is the range parameter, ν is a shape parameter, Γ is the gamma function, and K v is the modified Bessel function. Unlike CL and LR, OK does not assume sample independence and results in exact predictions at sampled locations (i.e., predictions equal observations at sampled locations). OK leverages spatial dependencies within the observed data for the prediction at unsampled locations. However, OK does not consider any auxiliary information in the spatial interpolation process. As the autocorrelation of the observed data decreases, the fitted models approach a pure nugget model, resulting in the prediction of the global mean at all locations. A more detailed description of ordinary kriging can be found in the study by Li and Heap [14].

Regression Kriging
Regression kriging is a hybrid spatial interpolation technique that combines the linear relationship between auxiliary information and the primary variable of interest estimated through LR with OK for spatial interpolation. RK links these two interpolative approaches together as consecutive steps. First, linear regression estimates the expected value as a function of the auxiliary information described above for LR. Then, semivariograms are fitted to the regressed residuals. The predicted value of the variable of interest at a given location is estimated following Equation (7): whereẑ(s) is the predicted value of the variable of interest at locations s, β o is the estimated intercept, β k are the estimated regression model coefficients, X k (s) are the values of the independent variables at locations s, λ i are the kriging weights, and (s i ) are the observed regression residuals at measured location s i . Because the regression estimation residuals have a mean of 0, RK utilizes simple kriging rather than ordinary kriging to estimate the weights, λ i . Depending upon the strength of the relationship of primary data to auxiliary data and the strength of the unresolved spatial autocorrelation of the residuals, RK will produce results similar to pure kriging or pure regression [45]. This has led some to suggest that pure kriging approaches and pure regression approaches can be considered as a special case of RK [46]. Given that RK can exploit both the spatial structure of the data and the extra information provided by the auxiliary information, RK theoretically should outperform either OK or LR. For a more detailed discussion of regression kriging approaches, see the study by Hengl [45]. Our implementation of RK used our LR models with the minimized AIC for each fuel component as the starting model structure.

Comparison of Interpolation Methods
We assessed the comparative performance of each SIM with a k-folds cross-validation approach [47]; in this manuscript, we set k to ten. The data were split into ten subsets, one subset was withheld, and the SIM was fit using the remaining data (i.e., k − 1). The withheld data were then utilized to assess the predicted values. We calculated the average error statistics for each of these ten iterations of cross-validation. We evaluated the performance of each SIM using the following statistics: mean error (ME); mean absolute error (MAE); mean absolute percent error (MAPE); and Pearson's correlation coefficient.

Results
The fuel sampling locations were well-distributed across the vegetation types, slopes, and aspects found on the PFDP ( Figure 1C). Thirty-two percent of the sample locations were in areas dominated by Engelmann spruce, 26% were in areas dominated by quaking aspen, 22% were in areas dominated by ponderosa pine, 17% were in areas dominated by Douglas-fir, and 4% were in areas dominated by blue spruce. The sample locations had an average slope of 9%, with a range from 1.6% to 22%, and were distributed across aspects.
The total surface fuel load ranged from 0.25 to 19.9 kg m −2 with a mean of 4.0 kg m −2 ( Table 1). The total fuel load was evenly divided between duff (48%) and dead down and woody fuels (45%), with litter making up the remaining 7% (Table 1). The dead down and woody fuel load increased with the time-lag size class, with the mean 100 h dead down and woody fuel load being three times greater than the 1 h dead down and woody fuel load. The mean fuel load was greater than the median for all the fuel components, indicating that the fuel load distributions were right-skewed. The fuel load was generally lowest in areas dominated by ponderosa pine and quaking aspen ( Table 2) and greatest in areas dominated by Engelmann spruce and blue spruce. The areas dominated by Douglas-fir tended to have intermediate fuel loads compared to the other vegetation types, except for the litter load, which was greater than the other cover types. Cover type was a significant predictor for all the components and total fuel load (p < 0.01 for all components; Supplementary Table S1). The best-fitting CL models had scales from 3 to 10 m with coefficients of determination (r 2 ) ranging from 0.04 to 0.29 (Table 3). Estimating the total fuel load and fuel load by component using CL resulted in prediction maps with stark transitions between the cover types ( Figure 2).    The forest structure and topographic variables were significant predictors in the linear regression models for all the fuel components and the total fuel load (Table 3). Cover type remained a significant predictor for the total fuel load and all components expect for the 10 h dead down woody fuel load. Local basal area or TPH were significant predictors for the total fuel load and all components except for the 100 h dead down woody fuels. At least one topographic variable was selected in the final regression for all the fuel components and the total fuel load. However, there was not consistency in the variable selected across the fuel components nor for the total fuel load, and the scale of the best-fitting models had a wider range than the CL models. When retained in a model, increasing the local basal area or TPH increased the fuel loading, the fuel load decreased as the distance to the nearest tree increased, and the fuel loading increased in areas of greater topographic sheltering (i.e., north aspects or greater topographic wetness index values; Supplementary Table S1). The inclusion of additional predictor variables beyond cover type resulted in predicted maps that generally captured the same patterns as those using CL but with variability within the vegetation classes (Figure 2).
Positive nugget-to-sill ratios for OK and RK indicated the presence of spatial correlation in all the fuel load estimates ( Table 4). The scale of autocorrelation (i.e., the range) varied from 6.08 to 39.70 m for OK and from 6.08 to 24.10 m for RK. The greatest spatial autocorrelation scales occurred for 1 h and 10 h down woody fuels and duff for both OK and RK (Table 4). RK had similar or lower ranges than OK, suggesting that the auxiliary variables accounted for some of the autocorrelation within the data. The differences in the predicted nugget values between OK and RK were minor, while the predicted sill values tended to be slightly lower for RK compared to OK. Where ranges exceeded 10 m, the predictive surface fuel maps based on OK captured the general patterns in the fuel load associated with vegetation class boundaries but with a smeared gradation along the transition between the cover types and greater variability in the fuel load than the maps based on CL and LR. However, the predictive fuel load maps based on OK for litter, 10 h down woody fuel load, and total fuel load did not indicate any discernible pattern associated with any overstory or topography. Instead, they illustrated a speckling pattern associated with the sampling locations and small range values. The surface fuel maps based on RK appeared as a hybrid of the LR and OK surface fuel maps with stark contrasts in the fuel load along the vegetation boundaries and less variation within a vegetation type than those based on OK (Figure 2). All of the SIM methods resulted in small mean errors, indicating that the predictions were unbiased ( Table 5). The MAPE varied from over 100% to around 40% depending on the specific fuel component and spatial interpolation method ( Table 5). Regardless of the SIM, the greatest MAPEs were associated with the dead, down, and woody fuel loads, followed by duff, litter, and total fuel load. RK resulted in the lowest MAPE for the 1 h and 100 h dead down woody fuels components, litter, duff, and total fuel load, with reductions ranging from 1 to 3% compared to the next best SIM and 3-12% for the worst SIM. For the 10 h dead down woody fuel load, OK resulted in the lowest MAPE, though there was only a 3% improvement over RK, which had the greatest MAPE of all the SIMs. RK had the greatest r obs,pred for the 1 h and 100 h dead down woody fuel components, litter, duff, and total fuel loads and was tied with LR for the greatest r obs,pred for the 10 h dead down woody fuel load. The r obs,pred for RK was from 1.17 to 2.54 times greater than that of the SIM with the lowest correlation coefficient. However, relative to the next greatest r obs,pred , RK had from 0 to 1.12 times greater correlation coefficients. Table 5. The mean error (ME), mean absolute error (MAE), mean absolute percent error (MAPE), and correlation coefficient between observed and predicted (r obs,pred ) from cross-validation.

Discussion
Fuel maps are commonly utilized by fire and land managers across a range of spatial scales to assist with the planning and locating of fire suppression opportunities, evaluating fire hazard and risk, designing wildland fuel treatments, and simulating fire behavior and effects [56][57][58]. Fine-scale variation in fuel characteristics is increasingly recognized as an important driver of fire behavior and effects [9,59], and as such, there is a growing demand for high-resolution maps of wildland fuels. Our results indicated that all four SIMs produced fine-scale (1 m 2 ) surface fuel maps with unbiased estimates. Although there were relatively minor differences in the predictive capability between the four spatial interpolation methods, regression kriging provided a lower MAPE and explained a greater proportion of the variance relative to the other approaches. Riech et al. [16] also found that the inclusion of a spatial autocorrelation term in a regression-based model resulted in improved predictions of surface fuel loads. These findings indicate that regression kriging may be a better choice for developing fine-scale surface fuel load maps than approaches relying solely on interpolation from point observations or regression-based approaches using auxiliary information.
Though there have been few comparisons of various SIMs for predicting fuel loads, there have been several previous studies that have implemented regression-based approaches to predict surface fuel loads from field-based and active and passive remote sensing data [13,33,60]. While previous studies have found moderate success with the prediction of coarse woody debris using linear and machine learning regression approaches, results for fine woody and litter surface fuels have been mixed. Bright et al. [61] reported r 2 values of 16 and 21 percent with a %RMSE of 39 and 80%, respectively, for litter/duff and 1-100 h surface fuel loads. Hudak et al. [62] explained 32% of the variation in surface fuel loading (combined duff, litter, woody, and herbaceous) in frequently burned longleaf pine (Pinus palustris Mill.) forests. Arellano-Pérez [63] explained 12% of the surface fuel load variability in Spain with a %RMSE of 35%. Lydersen et al. [11] found that the overstory structure explained 16-29% of the surface fuel variation in Californian forests. Our RK models resulted in r 2 and MAPE values that were well within the range of previously published studies despite differences in the scale, extent, and characteristics of the study's ecosystems and the sampling design and data collection methods used. The similar predictive capabilities across the various forest types, productivity gradients, and disturbance regimes found among these studies likely stem from the high spatial variability in the surface fuel load, as noted by Keane et al. [10] and Vakili et al. [32] among others. Given the large variations in and small ranges associated with surface fuel load, it might be possible to increase the overall accuracy through additional sampling. Though more studies are needed to better understand the effect of the sample size on fuel mapping accuracy and SIM selection, the large number of samples required would likely be impractical for fire management applications. Therefore, future mapping efforts might consider spending more effort on obtaining high-quality auxiliary information to improve the correlation in regression rather than increasing the sample size.
Conceptually, the spatial distribution of surface fuel load is controlled by a complex suite of biophysical and plant functional traits that influence productivity, litterfall, and decomposition rates [64]. Although we did not directly measure litterfall and decomposition on the site, we did find significant correlations with several indirect measures of decomposition and litterfall, including vegetation structure, composition, and topography. Overall, we found that the dominant species explained the most variation in our RK model, followed by metrics representing forest density and topography. These results were consistent with previous studies that have related litterfall and decomposition rates with species traits, forest structure, and environmental variables [11,[64][65][66][67]. These findings indicate that the statistical modeling of fuel accumulation based on relevant forest structure, composition, and biophysical variables can be effective at capturing the high spatial variability associated with the long-undisturbed mixed conifer ecosystems of the Rocky Mountains. However, the low correlation coefficients and large MAPE indicate that there is still room for improvement in fuel mapping.
There are several possibilities for increasing fuel mapping accuracy through the inclusion of other biophysical, forest structural, and disturbance history predictor variables describing the spatial and temporal dynamics of fuel accumulation and decay processes. One possible approach is to utilize repeat measurements of the forest structure through the incorporation of linked remote sensing technologies and forest inventory data. The use of repeat measures can supply estimates of both forest productivity and disturbances, both of which influence the spatial and temporal distribution of surface fuels. Though several studies (e.g., 61) have used repeated measures of relatively coarse-scale (<900 m 2 ) remote sensing data to develop fuel maps, there have been fewer studies that have utilized repeated measures to develop fine-scale fuel maps. At these finer scales, several studies have shown the promise of terrestrial and aerial LiDAR and uncrewed aerial system (UAS) structures for motion approaches for capturing forest structure parameters for mapping fuels and for capturing vegetative structural changes (e.g., [68][69][70][71][72][73]), but we are unaware of any studies that have utilized long-term repeated remote sensing measures for fuels mapping. With these three-dimensional remote sensing techniques providing spatially explicit estimates of individual trees, these data could provide a better correlation of how a forest structure impacts the input of fuels and microclimatic influences on fuel decomposition to improve surface fuel mapping. Further studies are needed to assess how repeated measures from three-dimensional remote sensing technologies can be used to improve the precision of fuel mapping and the characterization of fuel development relative to other mapping methodologies. Additionally, it is important to recognize that other approaches, including mechanistic-based simulation modeling of fuel accumulation (e.g., [74]), also need to be more broadly assessed and compared. Finally, hybrid fuel mapping approaches that link multiple data types (e.g., satellite-based remote sensing, terrestrial LiDAR, and destructive sampling) with mechanistic modeling may have the greatest potential to capture the spatial and temporal feedbacks that drive fuel accumulation, resulting in more accurate multi-scale fuel maps.
While there is considerable room to improve fine-scale surface fuel mapping methods, there is also a need to understand how improvements in surface fuel mapping effect fire behavior and effects predictions. Fine-scale surface fuel maps such as those we created and tested in this study are primarily used for fire behavior and effects predictions from next-generation three-dimensional process-based models. Further studies that assess the uncertainty and sensitivity of model predictions associated with surface fuel maps are needed to advance these models. Ideally, such evaluations would take place across a range of ecosystems, environmental conditions, and for both prescribed and wildland fires and would ultimately be tied into new experimental research aimed at validating the predictions from such models [75].

Conclusions
Surface fuel maps are used by land managers to predict fire behavior and effects during prescribed fire planning, assessing wildfire hazards and risk, and locating and designing effective silvicultural fuel-hazard-reduction treatments. Increasing the resolution of fine-scale fuel variability will enable the application of next-generation fire behavior modeling tools for fire management planning [58,75]. In this study, we compared the ability of four different spatial interpolation methods to capture the fine-scale (1 m 2 ) variability in the surface fuel load in a mixed conifer forest with terrain-driven variation in the vegetation and microclimate. We found that regression kriging outperformed the other approaches that relied solely on interpolation from point observations or regression-based approaches using auxiliary information for developing fine-scale surface fuel maps. Although our analysis found that surface fuel loading was correlated with species composition, forest structure, and topography, the relationships were relatively weak, indicating that additional predictors and spatial associations could significantly improve the surface fuel mapping. In addition to developing more accurate approaches to map surface fuels, additional research is needed to understand how uncertainty in surface fuel maps influences fire behavior and effects predictions across a range of model platforms, environmental conditions, and ultimately land management decisions associated with wild and prescribed fires.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/fire6060216/s1, Table S1: Anova summaries of classification and linear regression methods predicting total fuel load and fuel load by component. Table S2: Regression tables of classification and linear regression methods predicting total fuel load and fuel load by component. PIEN is the dummy-coded cover type. Were PIEN = Engelmann spruce, PIPO = ponderosa pine, PIPU = blue spruce, POTR = quaking aspen, and PSME = Douglas-fir.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to ongoing research efforts.