For forests in the western US, we characterized contemporary (1985–2021) fire regime distributions with satellite derived fire severity data and fire perimeters from the Monitoring Trends in Burn Severity (Eidenshrink et al. 2007) program. Historical distributions (36 simulated years approximating 1600–1880) were characterized using the LANDFIRE BioPhysical Setting product (BPS), which contains information regarding the average frequency of fires and proportion of those fires that burned under low, mixed, and high severity (LANDFIRE, 2020). Both the historic and contemporary datasets were 30-m resolution. We calculated the EMD between historical and contemporary distributions for fire frequency and fire severity, creating a distributional fire frequency departure (FFD), and distributional fire severity departure (FSD), which we then combined to obtain a multivariate distributional fire regime departure index (MFRD).
Study Area
We aggregated fire frequency and severity data estimated at 30-m resolution pixels (0.09 ha) to multiple spatial containers including hexagonal grids (hexel), firesheds (Ager et al. 2021), protected areas, among others (see data availability) across the western US. For simplicity, we walk through our analysis and main results using hexels. Our chosen hexel size (150,000 ha) is the median mapped areal coverage among biophysical settings (BPS) according to LANDFIRE (LANDFIRE 2020), which maps pre-euro American settlement vegetation with expert and data derived fire regime simulation models developed over the mapped extent. Hexels were not analyzed if they contained no contemporary forest fires and if they had greater than 50% overlap with the ocean, Canada, or Mexico. We also removed hexels if they contained less than 10% forested area, as fire regime departures are only calculated for forested pixels.
Frequency and severity distributions
For each hexel, we randomly sampled 0.1% of forested 30m pixels and iterated our analysis 100 times to capture variation and acknowledge the inherent spatial autocorrelation and pseudo replication associated with gridded spatial data (workflow of our analysis is shown in Supplemental 1; sensitivity tests for changes in parameters, such as the number of iterations, are in Appendix D). We defined forest as those pixels labeled as “Conifer,” “Hardwood,” “Conifer-Hardwood,” or “Hardwood-Conifer” by LANDFIRE Existing Vegetation Type and BPS. We filtered to forested areas because satellite-derived Composite Burn Index (CBI) is not appropriate for use in non-forest systems (Parks et al. 2019). To produce contemporary fire frequency and severity distributions, we first processed the western US fire perimeters from Monitoring Trends in Burn Severity (MTBS) using the approach described by Parks, et al., (2019) to obtain bias corrected CBI maps of every mapped fire greater than 404 hectares from 1985–2021. Parks et al., (2019) is a predictive Random Forest (Breiman 2001) model for satellite derived CBI using multiple spectral, climatic, and geographic variables in Google Earth Engine (Gorelick et al. 2017). These gridded CBI maps approximate field-based CBI measurements (Key and Benson 2006). Then we mosaicked individual fire maps into yearly fire maps from 1985–2021 and filtered to pixels (30m resolution) with > 0.35 pre-fire NDVI as an additional step to remove non-forest (Parks et al. 2023).
We then counted the number of times each pixel burned (CBI values > 0) from 1985–2021 to calculate fire frequency across the study area, from which we produced the contemporary fire frequency distribution within each hexel. Similarly, we extracted fire severity estimates for each burned pixel, which we then produced distributions for each hexel, with 16 evenly spaced bins with values ranging from 0–3.
To obtain historical fire frequency, we first extracted the mean fire return interval (MFRI) from the gridded BPS layer at each sampled pixel. Then we randomly sampled from a binomial distribution at each pixel with parameters n = 36, P = 1 / MFRI, n represents the number of years in the contemporary period and P is the yearly probability of fire in the pixel given its biophysical setting. We then summarized these historical fire frequencies into a distribution for each hexel.
To infer the severity of each estimated historical fire, we assigned a historical severity class of low, mixed, or high, via weighted random selection with weights given by the percent likelihood for each severity class as described by LANDFIRE BPS. For each assigned fire severity class, we then drew randomly from a uniform distribution within the bounds of CBI and thresholds between severity classes determined in Appendix B. For example, if we selected a low severity fire, we drew from a uniform distribution with a lower bound of zero and an upper bound of 1.75, the value determined from our threshold calculation. This produces an estimated CBI for each estimate fire. With these estimated fire severities, we then created fire severity distributions using the same 16 bins as described for contemporary fires for each hexel.
Fire-regime Departure
To quantify fire-attribute departures for frequency and severity in each hexel, we first rescaled our distributions by the historical mean and standard deviation, producing Z-scores. Then we computed the Earth Mover’s Distance (full description of EMD in appendix A) on each attribute independently, effectively quantifying the difference between the historic and contemporary distributions.
EMD can only be positive, so we modified attribute departures to indicate the net direction of change. We simply added a sign based on the change in percent burned at higher severity and MFRI. Positive values represent lower fire frequency, and higher proportion of high severity fire in the contemporary period compared to historical, respectively. This produces signed distributional fire frequency (FFD) and distributional fire severity departures (FSD).
Lastly, we calculated a multivariate distributional fire regime departure index (MFRD) by adding the fire frequency and severity departures together using Euclidean distance \(\:\sqrt{{\text{F}\text{F}\text{D}}^{2}+{\text{F}\text{S}\text{D}}^{2}}\). Because we simulate fire regime attributes using 100 iterations, the key outputs for each hexel are the median normalized departures for FFD, FSD, and MFRD. See our data availability statement for additional datasets and statistics for each spatial container.
Analysis Summaries
We summarized FFD, FSD, and MFRD by hexel, then grouped by state in the western US to analyze statewide patterns of fire regime departure. Next, we analyzed relationships between MFRD and human influence. For each hexel, we calculated the average human footprint (Venter et al. 2016), the proportion covered by public lands, and the proportion covered by Wilderness and National Park lands (U.S. Geological Survey 2022). Human footprint combines eight human pressures (e.g., structures) into one characteristic score of human influence on landscapes (Venter et al. 2016). Human footprint, proportion of public lands, and proportion covered by Wilderness and National Park all serve as indirect measures of fire exclusion policy and changes to the timing and frequency of anthropogenic ignitions (Balch et al. 2017; Boerigter et al. 2024). Each dataset was binned into five evenly spaced bins. We then calculated pairwise Wilcoxon rank sum tests to assess whether there were significant differences in fire regime departure between bins.
Additionally, to gain an alternative perspective on contemporary human use and its relationship to MFRD, we examined whether firesheds designated as priority landscapes are more departed than non-priority landscapes (US Forest Service, 2023). The fireshed dataset delimits landscapes with similar wildfire risks, land tenure, and planned management (Ager et al. 2021). Priority landscapes are designated by the US federal government as the firesheds (Ager et al. 2021) with highest need for future fire risk mitigation and the ability of communities to enact fire mitigation strategies (US Forest Service, 2023). For this analysis, we analyzed firesheds and then grouped them into priority and non-priority firesheds, determined by those with greater than 50% priority landscape coverage. We then calculated Wilcoxon rank sum tests between priority and non-priority firesheds to assess whether there are differences between these groups based on MFRD.
We compared our distributional metrics to a commonly used suite of departure metrics based on means – the Fire-Regime Condition Class (FRCC; Barrett et al., 2010). In each hexel, we used the fire frequency and severity information generated above to calculate historical and contemporary MFRI and proportion burned at high severity (Barrett et al. 2010; Johnson and Gutsell 1994). These values were then used to calculate FRCC fire frequency departure (FFDcc), fire severity departure (FSDcc) and fire regime departure (FRDcc) following the formulas in Barrett et al., 2010. Because FRCC applies no sign to their statistics, we converted FFD and FSD into absolute values. Next, we transformed these absolute values, along with MFRD and the FRCC departure statistics, into percentiles (denoted by a percentile subscript) relative to the western US, placing them all on the same 0-100 scale. We then subtracted the relevant FRCC departure from the relevant distributional departure (e.g., MFRDpercentile – FRDcc.percentile = Δpercentile) to assess the difference in percentile between our metrics and FRCC. This process was performed on all datasets (supplemental 3)
Finally, we provide two granular examples of our results within the Kalmiopsis wilderness and Olympic National Park. In these examples, we discuss the historical and contemporary MFRI, FFDpercentile, FFDcc.percentile and the underlying fire frequency distributions.
Package Credits
Spatial data was manipulated using the R packages sf and terra (Hijmans et al. 2022; Pebesma 2018; R Core Team 2022). Tabular data was manipulated using tidyverse, data.table (Barrett et al. 2023; Wickham et al. 2019), and overall analysis was performed with foreach, doParallel, and units (Daniel, Corporation, et al., 2022; Daniel, Ooi, et al., 2022; Pebesma et al., 2016). Figures were generated with Rcolorbrewer, ggplot2 (Neuwirth, 2022; Wickham, 2016) and in the Julia language with the packages OptimalTransport, Gadfly, and Distributions (Bezanson et al., 2017; Jones et al., 2021; Lin et al., 2023; Zhang et al., 2022). Calculation of EMD used the transport package in R (Schuhmacher et al., 2023).