Giving Ecological Meaning to Satellite-Derived Fire Severity Metrics across North American Forests

: Satellite-derived spectral indices such as the relativized burn ratio (RBR) allow ﬁre severity maps to be produced in a relatively straightforward manner across multiple ﬁres and broad spatial extents. These indices often have strong relationships with ﬁeld-based measurements of ﬁre severity, thereby justifying their widespread use in management and science. However, satellite-derived spectral for the direct mapping of CBI, which is more interpretable compared to spectral indices. Moreover, because the model and all spectral explanatory variables were produced in Google Earth Engine, predicting and mapping of CBI can realistically be undertaken on hundreds to thousands of ﬁres. We provide all necessary code to execute the model and produce maps of CBI in Earth Engine. This study and its products will be extremely useful to managers and scientists in North America who wish to map ﬁre e ﬀ ects over large landscapes or regions.


Introduction
Fire severity, defined here as fire-induced change to physical ecosystem components such as vegetation and soil [1,2], is a critical fire regime component that influences landscape heterogeneity, soil erosion, nutrient cycling, wildlife habitat, post-fire successional trajectories, and other ecological factors [3][4][5][6][7]. Consequently, there is great interest in mapping fire severity patterns [8,9], opportunities for which have greatly expanded with the availability of multi-spectral satellite imagery [10,11]. Satellite-derived spectral indices allow fire severity maps to be produced in a relatively straightforward manner across multiple fires and broad spatial extents [12] and have proven to be an invaluable resource for both the management and research communities. For example, satellite-derived indices have been used to evaluate temporal and spatial trends in fire severity [13][14][15][16], identify key factors driving fire severity patterns [17][18][19], and develop models to predict the potential ecological impacts of fire prior to a fire occurring [20,21].
Over the last few decades, several satellite-derived spectral indices have been developed to measure fire severity, including the delta normalized burn ratio (dNBR) [1], the relativized delta normalized burn ratio (RdNBR) [22], and the relativized burn ratio (RBR) [23]. Generally, these spectral indices are based on Landsat Thematic Mapper (TM) and Operational Land Imager (OLI) imagery, although several studies have evaluated spectral metrics produced with other sensors such as SPOT, AVIRIS, and Sentinel [24][25][26]. These indices quantify spectral differences between pre-and post-fire imagery and often have strong relationships with field-based measures of fire severity [27,28]. Most often, the spectral indices are justified as measures of fire severity through correlations with the composite burn index (CBI), an aggregated, field-based measure of severity [1], although correlations with specific ecological responses (e.g., canopy mortality, residual litter depth, crown charring) have also been used for validation [29,30]. The robust relationship between spectral fire severity indices and CBI provides strong support for their use in estimating fire effects across landscapes. However, criticisms consistently directed towards spectral fire severity indices are that their non-standardized units make them difficult to interpret in terms of on-the-ground fire effects [31][32][33] and the nonlinear relationship between field and spectral measures of fire severity complicates interpretation [34].
In an effort to provide ecological context to non-standardized fire severity spectral indices, several studies have categorized the indices into severity classes (e.g., low, moderate, and high severity) [29,35]. Although these thresholds are often based on observed relationships to field-based assessments of severity [15,28], this is not always the case [see 32]. Another viable approach to give spectral indices ecological meaning is to produce a statistical model in which a continuous field measurement of fire severity is modeled as a function of a given spectral index [27,[36][37][38]. This approach allows for the direct mapping of fire effects (e.g., CBI, changes to canopy cover) [39][40][41], which is arguably more interpretable than non-standardized spectral indices. The majority of these studies produce models relating a given field-based measure of fire severity to a single spectral index e.g., [20,37].
A small but growing number of case studies suggest that fire effects can be better characterized with machine learning models than by parametric regressions because they can easily incorporate multiple explanatory variables and their interactions [41][42][43][44]. The fundamental reasoning is that different aspects of fire-induced change can be characterized by different indices, and furthermore, machine learning techniques do not require a priori knowledge of variable interactions or functional forms. Incorporating variables sensitive to differences in physiography across broad geographic regions into machine learning models may also improve characterization of fire effects, as suggested by limited evidence that the relationship between field measures of fire severity and spectral indices may differ across biophysical gradients [23,33]. Consequently, models that incorporate diverse spectral indices and biophysical covariates to describe fire effects have the potential to provide more meaningful and accurate representations of fire severity, yet there have been no comprehensive studies across broad biophysical gradients to evaluate the applicability of such models.
Google Earth Engine is a cloud-based satellite image catalog and computing environment that allows users to conduct planetary scale geospatial analyses [45]. Earth Engine is becoming a widely used tool for measuring and quantifying continental-to-global characteristics such as annual gross primary production [46] and changes in the distribution of surface water [47] and land cover [48]. Earth Engine has also been used for mapping satellite-derived fire extent [49] and fire severity [50]. Machine learning models such as Random Forest regression and classification [51] are available in Earth Engine, and as such, models can be built and executed within Earth Engine, including those describing fire severity [44]. Given adequate field data, Earth Engine has the potential to predict and map fire effects (e.g., CBI) at continental to global scales.
We have three specific objectives, the first of which is to build a Random Forest regression model describing CBI (as a continuous variable) and evaluate several spectral, geographic, and climatic variables for inclusion into the model. Second, we aim to evaluate how the interpretation of spectral indices may vary over broad climatic or geographic gradients. Third, we intend to provide a framework and distribute code that allows users to execute the Random Forest model in Earth Engine, thereby enabling rapid production of mapped CBI predictions for fires that have occurred since 1984, with the exception of fires without available pre-or post-fire imagery.

Composite Burn Index Field Data
The composite burn index (CBI) is a post-fire field-based measure of fire severity, usually collected one year after fire [1]. Under the CBI protocol, over 20 individual factors are rated in each of five vertically arranged strata, including (a) substrates, (b) herbs, low shrubs, and trees < 1 m in height, (c) tall shrubs and trees 1-5 m in height, (d) intermediate trees (sub-canopy, pole-sized), and (e) big trees (upper canopy, dominant or co-dominant). There are four or five factors in each strata that include newly exposed soil/rock and duff consumption (stratum a), percent of vegetation that survived the fire and percent black/brown foliage (strata b and c), and canopy mortality and char height (strata d and e). Each factor is based on visual assessments and is rated from 0-3, with zero representing no fire effects and three representing the largest degree of change. Typically, ratings are averaged for each stratum and then across all strata to arrive at an overall CBI for an entire plot [1]. The CBI protocol is sometimes modified to better reflect fire effects in certain ecosystems (e.g., boreal forest; [27]). Figure  S1 shows the standard CBI form and all factors measured in each strata.
We acquired CBI data for fires with readily available fire perimeters that primarily burned forest vegetation (i.e., shrubland, grassland, and tundra fires were not evaluated) ( Figure 1). To do so, we obtained publicly distributed [52][53][54] and agency (e.g., US National Park Service and Parks Canada) CBI datasets; we also requested CBI data from specific individuals who collected such data for specific projects or uses e.g., [39]. Since we opportunistically acquired data collected by other institutions and organizations, we acknowledge that some forest ecosystems may have disproportionately low or high representation in our CBI dataset (Table S1). Nevertheless, the plot data spans broad climatic and geographic gradients of the US and Canada, making our model applicable to a large proportion of North American forests. The majority of the plots, with the exception of the eastern US, can be characterized as conifer forest, though some plots represent mixed forest (conifer and deciduous), deciduous forest, and non-forest. Most plots in the eastern US represent broadleaf deciduous forest with a lesser amount of conifer forest. Due to unique environmental settings or other phenomena, CBI data for several fires that burned forest vegetation were excluded from our analysis. For example, we excluded CBI data for seven fires in Georgia and Florida (US) because the spectral fire severity indices were affected by intermittent or ephemeral wetlands; that is, standing water was likely present in the pre-or post-fire imagery. CBI data for the Glacier Creek Fire (Alaska) were excluded since the pre-fire vegetation comprised beetle-killed spruce trees [55], and were therefore considered non-forested prior to the fire. We excluded CBI data from a handful of other fires (n = 5) whose CBI values were highly inconsistent with post-fire aerial imagery, potentially reflecting problems with the coordinates. We also excluded CBI data from an early season prescribed fire in Canada (Y-Camp). Lastly individual plots were removed if they were duplicates or had obvious issues (e.g., located well outside of the fire perimeter with a non-zero CBI).
We also accounted for two issues with the CBI data pertaining to unburned plots (CBI = 0). First, several fires had very few or no unburned field plots. Such plots are important for characterizing the spectral values of unburned pixels (or those with minimal fire effects), to adequately model CBI, and to better evaluate the models. Our goal was that the plot data for each fire comprised at least 5% of plots represented by CBI ≤ 0.25. Consequently, for those fires with <5% of the plots with CBI ≤ 0.25, we added 'pseudo-plots' representing CBI = 0. We produced pseudo-plots by randomly selecting unburned pixels at a distance of 200 m outside of the fire perimeter and selected the appropriate number of unburned pixels to obtain the 5% goal. Second, several CBI = 0 plots had very high spectral fire severity values; this was not expected and likely due to errors in either the plot coordinates or the CBI field data as entered. We removed CBI = 0 plots if the corresponding RBR value was greater than the 95 th percentile of randomly placed unburned pixels. Our final CBI dataset included 8075 plots covering 263 fires that burned from 1994 to 2017. Table S1 lists all fires used in our study and the source from which we obtained CBI data. For those fires that comprised multiple polygons (i.e., a complex), we considered each polygon as an independent fire for the purpose of this study.

Explanatory Variables: Spectral, Climatic, and Geographic Data
We evaluated several variables for inclusion in the model that can be characterized as spectral, climatic, or geographic (i.e., coordinates) ( Table 1). The spectral indices were derived from combinations of various bands from Landsat TM and OLI sensors and have been previously evaluated as measures of fire severity ( Table 1, Table S2). We produced the spectral indices from Landsat imagery using the Google Earth Engine cloud computing platform [45]. For each spectral variable, we used the mean compositing approach described by Parks et al. [50] to produce pre-and post-fire image composites. This approach selects all valid pixels (i.e., free of clouds, shadows, water, and snow) within a pre-specified time window (i.e., 'image season', Table 2) for one year prior to fire for pre-fire imagery and one year after fire for post-fire imagery. If valid pixels were unavailable during the corresponding image seasons, we extended the time frame to include imagery from two years prior or post-fire, as needed. All pixels meeting our criteria were stacked and averaged on a per-pixel basis. We did not include a dNBR offset [1] as fire perimeter accuracies were variable in Canada, reducing our confidence in offsets from unburned areas. Although most of the spectral variables rely on differences between preand post-fire image composites, two are characterized by post-fire image composites only (Table 1).  (Table S2). Note we did not evaluate dNBR or RdNBR because (1) RBR, dNBR and RdNBR were highly correlated (Pearson's r = 0.99) (Table S3) and (2) RBR performed marginally better than dNBR and RdNBR when validated with the CBI (Table S4). Longitude (rounded to the nearest degree) NA Climatic gradients are characterized by two variables representing the climatic water balance [60], actual evapotranspiration (AET) and the climatic water deficit (CWD). Both AET and CWD were obtained from the TerraClimate dataset [61] (resolution = 4 km) and represent climatic normals spanning 1981-2010. These variables are generally considered strong controls on the distribution of vegetation types [62,63]. Geographic coordinates are simply represented by latitude and longitude, rounded to the nearest degree.

Random Forest Model
Within Earth Engine, we built and evaluated a Random Forest model describing CBI as a function of multiple variables. Our goal was to produce a parsimonious model with a reduced set of variables, with each variable contributing to the prediction of CBI. To do so, we followed the cross-validated stepwise procedure described by Parks et al. [21]. Specifically, we built an initial model with all explanatory variables (Table 1) using a five-fold cross-validation in which data from 80% of the fires were used to build the model; CBI values were predicted for the plots representing the remaining 20% of the fires. This resulted in CBI predictions for all plots in which the plot itself was not included in the model. We then calculated the coefficient of determination (i.e., the R 2 of a linear model between observed vs. predicted CBI) for all plots.
Next, we built subsequent Random Forest models in which each variable was iteratively removed; each of these models was validated using the identical five-fold cross-validation. If the cross-validated R 2 did not decrease by at least 0.005 when any given variable was removed (in some cases, the R 2 increased), it indicated that the variable did not contribute any unique information. As such, any individual variable that resulted in <0.005 decrease in R 2 when excluded from the model was permanently removed from consideration. This process was repeated (i.e., variables were removed) until all remaining variables resulted in a decrease in R 2 of at least 0.005 when removed from the model. This procedure resulted in a parsimonious model in which each variable provides a non-negligible contribution. We report the cross-validated R 2 , root mean square error (RMSE), and the mean absolute error (MAE) of the final model for observed vs. predicted CBI. All predictions were produced using the test datasets from the five-fold cross validation, such that predictions were produced for independent data not used to build the models. We also report these validation statistics for various geopolitical and ecological regions, thereby providing an indication of the spatial variability in model skill.
We initially used the default Earth Engine Random Forest parameters, except for the number of trees, which we set to 500. However, initial explorations showed these Random Forest parameters overfit the relationship between CBI and the explanatory variables ( Figure 2). We reduced overfitting by manipulating the minLeafPopulation parameter, which is the minimum size of a terminal node (default = 1). We set the minLeafPopulation parameter equal to: where num.recs is the number of data samples and num.vars is the number of variables in the model. This equation was based on explorations using a 5-fold cross-validation; this approach substantially reduced overfitting (Figure 2). Initial explorations also indicated that the final model was biased, in that CBI was overpredicted at low CBI values and underpredicted at high CBI values ( Figure 3). This is because the model fits the mean response, and data extremes rarely represent the mean response. The bounded nature of the CBI (range: 0-3) also plays a role. Biased predictions may complicate interpretations of mapped CBI predictions, particularly when defining classes such as low-and high-severity e.g., [17,64]. For the purpose of mapping predicted CBI, we bias corrected the CBI predictions as follows: where CBI bc is the bias-corrected CBI prediction and CBI pred is the predicted CBI from the Random Forest model. The initial subtraction of 1.5 centers predicted CBI at zero, the 1.3 and 1.175 multipliers provide the bias correction, and the final 1.5 addition re-centers the predictions at 1.5. When CBI bc resulted in predictions less than zero or greater than three, they were truncated to zero and three, respectively. We defined the parameters for this bias correction by comparing histograms of the observed CBI vs. CBI bc ( Figure S2). The bias-corrected CBI better characterizes the low and high end of the CBI gradient ( Figure 3). Note that all results shown in the main body of this paper do not use the bias-corrected CBI, with the exception of the mapped predictions.

Model Implementation in Earth Engine
We designed the procedure to run entirely within Google Earth Engine [45]. This includes the production of all spectral explanatory variables, executing the model, and producing mapped CBI predictions. We provide all code and a sample fire history shapefile to produce mapped CBI predictions for 10 sample fires. The code produces two CBI predictions for each fire, one without the bias correction and one with the bias correction. Interested users can easily modify the code to produce mapped CBI predictions for nearly all fires that have burned since 1984 in North American forests, with the exception of fires without available pre-or post-fire imagery. The code is available here: https://tinyurl.com/CBImodel. . We applied a bias correction to address this issue (Eqn. 1). This bias correction only applies to mapped CBI predictions (Figure 9), though the code distributed with this paper produces predictions with and without the bias correction. Except for Figure 9, all tables and figures refer to the predictions without bias correction.

Evaluating the Potential for Spatial Variability in the Relationship between CBI and Spectral Indices
Lastly, we aimed to illustrate how variability in geography and climate influences interpretation of spectral fire severity indices. To do so, we produced two generalized additive models (GAMs), one as a function of RBR and latitude and the other as a function of RBR and CWD. Latitude and CWD were chosen because the final Random Forest model (after the cross-validated variable selection procedure) included these non-spectral variables (see Results). Using each model, we then produced plots holding CBI constant at incremental values (e.g., 1.25 and 2.25) while varying RBR and latitude and RBR and CWD.

Results
The Random Forest model describing CBI as a continuous variable performed well across North America (Figure 4; Table 3). The cross-validated R 2 , RMSE, and MAE were 0.72, 0.47, and 0.37, respectively. The spectral variables, particularly the RBR, provided the largest contribution to the models ( Table 3 and Table S4). The cross-validated stepwise procedure we used for variable selection resulted in a parsimonious Random Forest model describing CBI as a function of RBR, dMIRBI, dNDVI, post.MIRBI, CWD, and latitude. Some spectral variables (i.e., dEVI, dNDMI, post.NBR), AET, and longitude were excluded as a result of the cross-validated model selection procedure.   Tables 4-6). For example, CBI was modeled well (e.g., high R 2 , low RMSE, and low MAE) in Arizona, California, New Mexico, and Alberta ( Figure 5; Table 4), but our model was weaker in states such as Virginia and Alaska. Although model skill for some geopolitical units was fairly poor, small sample sizes make the validation statistics uncertain (e.g., Arkansas and North Carolina). Also, the CBI data in some geopolitical units only captured a fairly small portion of the possible CBI range, adding to the uncertainty in model skill. For example, the CBI data for Missouri spanned from 0.10 to 1.65, thereby omitting about half of the range in CBI. Issues of low sample size and inadequate sampling of the CBI gradient were mostly evident in the southeastern US. When all cross-validated CBI predictions were compared with observed CBI within broad geopolitical units, our model performed best in the southwestern US and Canada ( Figure 6; Table 5) and rather poorly in the SE US. In relation to ecoregions [65], the model performed best in North American Deserts, Northern Forests, and Temperate Sierras (Figure 7; Table 6). Assessments of model skill for the bias-corrected version are provided in the Supplementary Materials ( Figures S3-S5, Tables S5-S7). The relationship between CBI and RBR varied along geographic and climatic gradients in North America. Specifically, the RBR representing specific CBI values (e.g., the threshold between low/moderate severity [CBI = 1.25] and moderate/high severity [CBI = 2.25]) increased with latitude and decreased with climatic water deficit (CWD) (Figure 8).    Table 6. Cross-validated model skill (observed vs. predicted CBI) for ecoregions (see Figure 1) in our study.   . This indicates that spectral fire severity indices such as the relativized burn ratio have slightly different meanings across geographic and climatic gradients. For example, at southern latitudes, a given RBR value corresponds to lower CBI compared to northern latitudes, and in regions that are less moisture limited (low CWD), a given RBR value corresponds to higher CBI compared to moisture limited regions with high CWD.

Discussion
Describing fire effects with satellite imagery is an active area of research across the globe [40,66,67]. Most of this research, however, has focused on describing fire severity with one spectral index (e.g., dNBR). The focus on a single spectral index has two important consequences. First, the de facto strategy when producing fire severity maps is to show the non-standardized index [12]; these non-standardized indices are difficult to interpret in terms of on-the-ground fire effects [36,39]. Second, there is often a failure to recognize that multiple spectral indices or other bioclimatic variables may improve models describing fire severity [33,44]. In this study, we address both limitations by a) directly modelling and mapping CBI, and b) including multiple variables into our models, thereby better describing CBI across forested landscapes in the US and Canada. Consequently, our study provides an essential refinement to mapping fire severity (as it relates to ecological impacts) across North America. Also, because we distribute code that allows users to produce mapped CBI predictions within Google Earth Engine, scientists and practitioners have the ability to map fire effects over hundreds to thousands of fires (Figure 9). Figure 9. Example of predicted CBI for select fires. Note that we do not have CBI data for all of these fires, illustrating that CBI predictions can be produced for fires that have occurred since 1984, provided pre-and post-fire TM and OLI imagery are available.
Our study showed that interpretation of spectral fire severity indices varied to some degree over bioclimatic gradients, in that RBR values corresponding to the thresholds between severity classes vary by latitude and climate (Figure 8). A possible explanation is that latitude and climate are indirect proxies for broad-scale forest type, and that interpretations of spectral indices may vary among forest types e.g., [39]. Nevertheless, our findings agree with those of Harvey et al. [33], who found that identical spectral values across latitudes corresponded with different field-measured fire effects. This finding suggests that evaluations of spectral fire severity across broad region [17,21,68] could be slightly misinterpreted without implicitly accounting for biophysical factors that influence interpretation of spectral indices. The changing interpretation between field-based fire effects and spectral indices provides ample rationale that biophysical variables should be explicitly evaluated for inclusion into models describing fire severity across broad spatial domains, as we did in this study.
The relativized burn ratio (RBR) was clearly the most important variable that we evaluated in our models (Table S4). Although we did not evaluate dNBR and RdNBR for inclusion into the model due to their very high correlation with RBR (Table S3), either would likely have been the most important variable were they used instead of RBR. The addition of other spectral variables (dMIRBI, dNDVI, post.MIRBI), CWD, and latitude provided a modest improvement compared to the model based only on RBR; the full model improved the cross-validated R 2 by 0.04. Although we were somewhat surprised that the improvement of the full (multi-variable) model was not greater, it could be that including additional variables results in model improvement in specific CBI ranges (e.g., low severity), as previously demonstrated [44]. Nevertheless, our study illustrates the utility of using multiple variables to model fire effects and supports the findings of several case studies to this effect [41][42][43][44].
Many studies have demonstrated that the relationship between field-based fire effects and spectral fire severity indices varies across fires and vegetation types [55,69,70]. This has led to questions as to whether or not fire severity indices, or more specifically, the functional relationships between fire effects and a spectral index, are 'transferable' across fires and geography [71]. Although such models may indeed be transferable to some degree [72], it is generally the case that model skill decreases when multiple fires are evaluated simultaneously due to variable functional relationships (among individual fires) between field-measured fire effects and spectral indices [23,39]. We posit that our model describing and mapping CBI increases model transferability among fires and regions. Indeed, our model performed well across North America (cross-validated R 2 = 0.72 for all 8075 plots) even though our study area spanned a broad gradient of the biophysical environment. This is an important outcome of our study given that it is unrealistic to collect field data on all fires and build individual models for each. This said, our models do not perform well everywhere (e.g., SE US; see Tables 4-6), but this is generally due to poor relationships between field and spectral data in certain localities and not a reflection of the model itself.
The main purpose of our study was to build a parsimonious model measuring (or estimating) fire effects (i.e., fire severity). This is a nuanced but important difference from studies that aim to explain fire effects e.g., [17]. The former aims to simply quantify the ecological impact of fire without regard to the causes or drivers, whereas the latter is intended to identify factors that control fire severity such as topography, weather, and pre-fire fuel [38,[73][74][75]. Since our goal was to measure fire effects, we did not include variables such as topography that potentially drive or explain fire effects. Studies that explain fire effects can be particularly useful for identifying factors that can potentially be modified to change fire severity outcomes (e.g., fuel reduction treatments) [76] or making predictions of expected fire severity prior to a fire even occurring [20,21]. We suggest that our model that measures fire effects (i.e., CBI) can be used by other researchers as a dependent variable in models that explain fire effects. This would provide more ecological meaning and an improved interpretation when, for example, evaluating drivers of fire severity or measuring reductions in fire severity from fuel treatments.
Although our models performed well in estimating field-measured CBI across broad North American geographic and climatic gradients, there are some caveats to this approach and the resulting model. First, this model should not be used to predict CBI outside of North America because (a) the geographic and climatic variables in other regions would be outside of the range used to build the model and (b) it is not clear if the relationship between CBI and the explanatory variables are similar on other continents. CBI data from other continents could be incorporated to produce a global model or a model applicable to each continent. Second, because we purposefully omitted CBI data for those fires whose pre-fire vegetation was primarily non-forest, the applicability of our models is limited to forested ecosystems. Third, our model did not perform particularly well in the SE US (R 2 < 0.5), an issue not unique to our study [77]. The poor performance could be due to several factors, including that most fires in the SE US in our CBI dataset were prescribed burns that may have occurred during seasons or under conditions in which fires do not typically burn. The deciduous nature of much of SE US forests could have also posed problems in terms of the spectral signatures of different phenological states (leaf off vs. leaf on). Rapid post-fire regrowth in deciduous forests may also obscure spectral fire severity assessments. Continued study into the proper timing of pre-and post-fire imagery for describing fire severity in deciduous forests would be beneficial.
It is worth noting that the CBI itself has been criticized as an imperfect measure of fire severity [34]. For example, the CBI is a post-fire assessment in which pre-fire vegetation conditions are unknown but must be estimated [78]. The CBI is also qualitative in that each factor is visually estimated and not measured per se [79]. Furthermore, the CBI aggregates several diverse fire effects into a single index, thereby obscuring specific factors of interest [31]. Lastly, in some forest ecosystems, particular factors may be more important in describing fire severity (e.g., organic matter consumption in boreal forests), even though they are given equal weight in the CBI [80]. Nevertheless, because the CBI protocol is fairly rapid to implement (~30 min/plot) [1], two field crews can collect 50+ plots in a week (depending on transportation and hiking time), thereby enabling fairly rapid and inexpensive collection of field data that can be correlated with spectral indices or used to build models such as that described in this study. In spite of the criticisms levied towards the CBI, the consistently strong correlations between CBI and spectral indices reported in the literature are a testament to its utility in fire studies e.g., [23,25,27,38].

Conclusions
This study introduces a publicly available model with spectral, climatic, and geographic explanatory variables that produces gridded estimates (i.e., maps) of the composite burn index (CBI), a frequently used, field-based measure of fire severity. The ability to map across a broad geographic range was achieved through the opportunistic compilation of field-based CBI collected previously for a range of studies. The model developed in our study used data from 263 fires and covered broad gradients in the geography and climate of North America to describe the CBI. Maps of CBI produced with this model are more interpretable in terms of on-the-ground fire effects compared to non-standardized spectral indices. Since our model and the resulting maps give ecological meaning to spectral indices, they will provide managers a better understanding of fire effects. The robust relationships and the fairly high model skill in most regions suggest the resulting CBI maps will be particularly beneficial in remote regions where it is expensive and difficult to acquire field measures of severity (e.g., Alaska and the majority of Canada). This said, additional CBI plot data in underrepresented areas of the continent would likely improve model skill. The full Random Forest model provided a moderate improvement in model skill compared to the RBR-only model. This is because the interpretation of fire severity indices varies geographically [33] and, furthermore, different spectral indices complement and enhance assessments compared to use of one spectral variable alone [44]. More importantly, because we provide code to fully execute the model in Google Earth Engine, we provide managers and researchers with a tool to evaluate fire effects to support management decisions and advance ecological understanding at regional, national, and continental scales of North America. Gridded maps representing estimated CBI offers the potential to assess fire effects that may be applied to a range of user-specified natural resource applications.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/14/1735/s1, Table S1: List of fires with CBI data used in this study, Table S2: Equations for spectral indices used in CBI model prediction, Table S3: Pearson's correlation among the variables we evaluated, Table S4: Cross-validated model skill for Random Forest models describing CBI with each spectral variable, Table S5: Cross-validated model skill (observed vs. predicted CBI [bias-corrected]) for states, provinces, and territories, Table S6: Cross-validated model skill (observed vs. predicted [bias-corrected]) for large geopolitical regions, Table S7: Cross-validated model skill (observed vs. predicted CBI [bias-corrected]) for ecoregions, Figure S1: Composite burn index (CBI) field sheet, Figure S2: Histograms show the distribution of observed CBI, predicted CBI, and bias-corrected CBI, Figure S3: Observed vs. predicted CBI (bias-corrected) for the Random Forest model for each state, province, and territory, Figure S4: Observed vs. predicted CBI (bias-corrected) for the multi-variable Random Forest model for large geopolitical regions, Figure S5: Observed vs. predicted CBI (bias-corrected) for the multi-variable Random Forest model for ecoregions.