Before the fire: predicting burn severity and potential post-fire debris-flow hazards to conservation populations of the Colorado River Cutthroat Trout (Oncorhynchus clarkii pleuriticus)

Study area

Our study covers the putative native range of bCRCT (Metcalf et al. 2012; Bestgen et al. 2019) including parts of Colorado, Wyoming, and Utah (Fig. 1). We predict pre-fire burn severity and potential PFDF hazard across this distribution, but also present more detailed analyses focused on the CP, loosely defined as those displaying less than 10% admixture with non-native trout determined from molecular study (Hirsch et al. 2013). These CP are spatially nested within the historic and current distribution of bCRCT (Fig. 1). We selected watersheds from the Watershed Boundary Dataset WBD (2023) 12-digit Hydrologic Unit Codes (HUC-12) within 1 km (n = 723) of the historic distribution of bCRCT delineated by Hirsch et al. (2013). The extent (Fig. 1) covered roughly 7.5 million ha and includes five vegetation types (Kuchler 1964): (1) Great Basin shrub/steppe (35%); (2) Rocky Mountain conifer forest (34%); (3) Great Basin/Southwest forest (21%); (4) Great Basin shrub/grassland (8%); and (5) alpine meadow (2%). Our analysis involved multiple steps; these are summarised below and in Fig. 2.

Fig. 1.

Study area extent derived from 12-digit Hydrologic Unit Code (HUC-12) watershed boundaries within 1 km of the current and historic distributions of blue lineage Colorado River Cutthroat Trout-Blue Lineage (bCRCT), 4-digit HUC (HUC-4) watershed boundaries, Monitoring Trends in Burn Severity (MTBS) wildfire perimeters (n = 28) from 2016 to 2020, and nested distribution of bCRCT conservation populations. Note: map derived from fig. 1 in Bestgen et al. (2019).

Fig. 2.

Flowchart of analysis steps including data inputs, linkages, algorithms, and hierarchy to model and quantify hazard of potential burn severity and potential post-fire debris flow to conservation populations of blue lineage Colorado River Cutthroat Trout (bCRCT; Oncorhynchus clarkii pleuriticus) in unburned areas.

We identified 28 wildfires (Table 1) with perimeters that intersected the study area spanning 2016–2020 (Fig. 1) from the inter-agency database Monitoring Trends in Burn Severity (Eidenshink et al. 2007; MTBS [Monitoring Trends in Burn Severity program] 2023). We selected fires from 2016 onward to align with imagery available from the Sentinel-2 mission (European Space Agency 2023). We modelled pre-fire burn severity for individual wildfires and cast averaged predictions by associated HUC 4-digit (HUC-4) watershed boundaries that intersected our study area (Fig. 3). This included the Upper Green (8 wildfires), Lower Green (8 wildfires), Yampa (12 wildfires), and Dirty Devil watersheds. We modelled pre-fire burn severity in the Dirty Devil HUC-4 watershed, but because no large wildfires occurred there during 2016–2020, pre-fire burn severity was modelled based on wildfires from the closest adjoining watershed, the Lower Green. To quantify direct effects of potential wildfire, we tabulated predicted pre-fire burn severity within 500-m buffers around streams identified by the Geographic Names Information System (GNIS [Geographic Names Information System] 2023) of bCRCT CP, which were composed of 70% Douglas fir/spruce/fir (Psuedotsuga menziesii, Picea spp., Abies spp.) based on potential natural vegetation in the conterminous United States (Kuchler 1964). We quantified indirect effects of wildfire by modelling potential PFDF at the scale of stream reach associated with bCRCT CP (Fig. 3) and tabulated hazard to individually known localised populations of bCRCT CP by name (node ID).

Fig. 3.

Predicted differenced Normalised Burn Ratio (dNBR) differentiated by 4-digit Hydrologic Unit Code (HUC-4) watersheds (Dirty Devil, Lower Green, Yampa, and Upper Green) for contemporary conditions (September 2022) used to assess hazard of potential post-fire debris flow for conservation populations of blue lineage Colorado River Cutthroat Trout (bCRCT; Oncorhynchus clarkii pleuriticus) in unburned areas.

Predicting burn severity

We developed machine learning model of burn severity based on pre-fire conditions, henceforth referred to as predicted burn severity. Each wildfire was modelled individually using observed post-fire MTBS dNBR pixel values as a response variable. Research shows wildfire depends upon multiple interacting environmental factors, including: (1) physical ecosystem conditions (Hagmann et al. 2021); (2) fuel (Keeley and Syphard 2019); (3) topography (Povak et al. 2018); and (4) weather (Litschert et al. 2012), which formed the basis of our predictor variable development.

We developed predictions of potential burn severity using pre-fire environmental factors as previously demonstrated in Wells et al. (2023). To characterise the pre-fire physical ecosystem, we used Sentinel-2A/B (Drusch et al. 2012) imagery from Copernicus Open Access Hub https://scihub.copernicus.eu and the Sentinel Applications Platform (SNAP ver. 7.0; European Space Agency 2022) to pre-process digital imagery and create predictor variables. We created predictor variables for two periods prior to each wildfire. The first period was from as close to the ignition date as possible and the second was an early season period prior to the onset of the growing season and free of cloud or snow over the wildfire perimeter. Use of two periods provided a more robust range of physical ecosystem conditions and inherent delineation between deciduous and evergreen vegetation cover. Sentinel-2 images were corrected based on user specified inputs (Stand-alone sen2cor ver. 2.9; STEP 2024) for atmosphere, bidirectional reflectance, and topography and converted to reflectance. Individual spectral bands (1–8a; 11–12) were used along with Normalised Difference in Vegetation Index (NDVI, Huang et al. 2021), Soil-Adjusted Vegetation Index (SAVI, Rondeaux et al. 1996), and biophysical parameters of leaf area index, fraction of absorbed photolytically active radiation, chlorophyll content in leaf, fraction of vegetation cover, and canopy water content that were included in past work (Wells et al. 2023).

Fuel and topographic information were acquired from existing data products. To characterise surface fuels, we used yearly estimates of herbaceous cover (HERB) of annual and perennial grasses, and forbs (Rigge et al. 2022). To estimate overstorey fuels, we used LANDFIRE (2016) LF 1.4.0 Forest Canopy Cover (FCC), representing fuels prior to the earliest wildfire. Elevation variables were extracted from the National Elevation Data products at 1/3 arc second (Gesch et al. 2018). Additional topographic variables of slope, terrain position index, terrain roughness index, surface roughness, and cosine and sine transformation of aspect were calculated from elevation.

For weather information, we extracted gridded surface meteorological data (Abatzoglou 2013) and averaged daily values over three time-lags prior to ignition: (1) 7 days; (2) 30 days; and (3) one calendar year. These time-lags represented differing climate, fire weather, and fuel conditions that we expected to influence burn severity. Predictors included 100-h and 1000-h dead fuel moisture, specific humidity, maximum and minimum relative humidity, maximum and minimum air temperature, precipitation, wind speed, vapour pressure deficit, downward surface shortwave radiation, alfalfa and grass reference evapotranspiration, energy release component and burning index.

For each fire, we used a geographic information system (GIS) to sample pixels of MTBS datasets of observed dNBR and associated pre-fire predictor variables. We randomly sampled pixels within each fire perimeter to develop model building (training-testing; 40%), and validation (reserved; 10%) partitions. We used a random forest (Breiman 2001) in Anaconda Python 3.10.0 with the scikit library to fit and validate models of dNBR with the training-testing and validation data partitions, respectively (Pedregosa et al. 2011; Anaconda Software Distribution 2023). We removed highly correlated (R > 0.70) predictor variables with a forward stepwise routine (Sherrouse and Hawbaker 2023) and tested for spatial auto-correlation with regression kriging (Wells et al. 2023). The random forest algorithm used a 10-fold cross-validation (CV) based on root mean square error (RMSE) of the CV test fold, with 500 estimators and a minimum of 10 samples per leaf. To assess variable contributions to model compositions, we tallied variable frequency and normalised variable importance values within each model to calculate the summed proportional variable importance (SPVI) over all models. SPVI is calculated by taking the importance value of each selected variable per individual model, divided by the total variable importance of all variables from that individual model, and summing that particular variable over all models. For each variable we tabulated percentiles of SPVI over all models and averaged SPVI by environmental factor (physical ecosystem, fuel, topography, and weather). Spatial extrapolation of pre-fire burn severity, or model predictions, were based on the variables from each individual wildfire model with data acquired representing a contemporary period (11–19 September 2022), with early season imagery variables matched to May–June 2022. For each wildfire model, predictions were cast over the associated HUC-4 watershed and averaged by arithmetic mean. Averaged predicted dNBR by HUC-4 watershed were classified based on natural breaks (Jenks and Caspall 1971) into four severity classes to mirror MTBS burn severity (Fig. 4). Averaged predictions were merged and scaled for display to show relative dNBR over the entire study area (Fig. 3). To verify model functionality, we used hindcasting of the largest individual wildfires (East Fork) based on the Dollar Ridge Fire model (Fig. 5), which occurred within the same watershed (Lower Green) the year prior (Table 1) and compared predicted severity classes to observed MTBS post-fire burn severity. We tabulated a confusion matrix of predicted pre-fire burn severity classification (binned moderate and high) against observed MTBS thematic burn severity (MTBS [Monitoring Trends in Burn Severity program] 2023) and reported the overall accuracy and Kappa statistic.

Fig. 4.

Predicted burn severity, as of September 11–19 2022, based on classified predicted differenced Normalised Burn Ratio (dNBR Fig. 3) for portions of the (a) Yampa, (b) Lower Green, (c) Upper Green, and (d) Dirty Devil watersheds containing Colorado River Cutthroat Trout (Oncorhynchus clarkii pleuriticus) shown in style of Monitoring Trends in Burn Severity (MTBS) data and as used for input to estimating potential post-fire debris flow.

Fig. 5.

Hindcast of one wildfire model to another wildfire showing (a) pre-fire predicted differenced Normalised Burn Ratio (dNBR) of the East Fork Fire derived from the model of the Dollar Ridge Fire (blue-red = low-high), (b) predicted burn severity (green = unburned or low, red = moderate or high), (c) post-fire observed dNBR of the East Fork Fire from the Monitoring Trends in Burn Severity (MTBS) program (blue-red = low-high), and (d) observed MTBS burn severity (green = unburned or low, red = moderate or high). Note: blank spot in 4A and 4B are due to missing values.

Potential for post-fire debris flow

We used the U.S. Geological Survey (USGS) post wildfire debris flow hazard assessment (U.S. Geological Survey 2023) to estimate potential for PFDFs based on our predictions of burn severity. The hazard assessment uses empirical models to estimate the likelihood (%) and potential volume (m³) of PFDF in response to a series of simulated storms with 15-min peak rainfall intensities (12–40 mm h⁻¹ precipitation by 4 mm h⁻¹ increments), which may be expected within a 1-year recurrence interval (NOAA [National Oceanic and Atmospheric Administration] 2023). We delineated the stream reach drainage network and associated catchments within each HUC-12 watershed based on 10-m elevational data and GIS stream analysis. PFDF model calculations were performed at each catchment outlet and along upstream drainage networks by stream reach. Debris-flow likelihood calculations used logistic regression (Staley et al. 2016) and our inputs of predicted pre-fire burn severity classified as moderate or high (described above), soil properties (K-factor of the universal soil loss equation; Wischmeier and Smith 1965) and rainfall accumulation for a given duration. Potential debris-flow volumes were predicted using multiple linear regression models (Gartner et al. 2014) with simulated rainfall intensity, basin morphology, and predicted burn severity (classified as moderate or high, as described above). Volume estimates were classified by order of magnitude from 0–1000 m³, 1000–10,000 m³, 10,000–100,000 m³, and greater than 100,000 m³. We combined likelihood and volume debris flow results (Table 2) following Cannon et al. (2010).

Conservation population hazard

We quantified hazard to bCRCT CP by assessing both direct effects of potential wildfire burn severity and indirect effects of potential PFDF. We spatially aggregated predicted burn severity and potential PFDF to stream reaches and CP node ID, respectively. We buffered (500 m) stream reaches of CP and quantified averaged predicted burn severity, (dNBR), by stream name following nomenclature from the GNIS since single CP often occupy more than one named stream. Since multiple PFDF stream segments conjoin stream reaches of bCRCT CP, we summarised potential PFDF linear stream segments by averaging to points at the intersections of stream outlets with mainstem segments. The mainstem outlet points were averaged over stream reaches, based on rankings of none (0), low (1), moderate (2), and high classification (3) for increasing simulated storm intensities, by bCRCT CP node ID. This final step (Fig. 2) where we averaged linear stream segments to outlet points on mainstem streams and mainstem outlets by stream reach, created a hazard index value of potential PFDF, per simulated storm intensity, for each bCRCT CP node ID.

Predicting burn severity

Over the study area, the 28 wildfires analysed burned 159,000 ha from 2016 to 2020 (Table 1). Observed mean dNBR from individual wildfires ranged from a low of 0.13 to a high of 0.49 with an average of 0.29 (Table 3). Average RMSE values of the random forest model of pre-fire dNBR, based on validation dataset, were 57% of the averaged observed mean dNBR (Table 3). One model (Fly Canyon) showed large variation of >100% RMSE of observed mean dNBR. Based on validation data, the models had adjusted R² ranging from 0.17 to 0.73, with an overall average adjusted R² = 0.54 (Table 3). There was no significant difference in CV-test RMSE between any of the random forests models and regression kriging models, suggesting model residuals had little to no spatial auto-correlation. Model functionality based on hindcasting the East Fork Fire with the Dollar Ridge model showed strong spatial congruency compared with observed dNBR and burn severity from post-fire MTBS (Fig. 5). Prediction of categorised burn severity of the East Fork Fire based on the Dollar Ridge model was marginally accurate (overall accuracy of 0.64, κ = 0.28). Random forest models predicted 2.6 million ha (34%) of the study area (Fig. 3) as moderate to high burn severity. Categorised burn severity by HUC-4 watershed (Fig. 4) showed the Yampa (Fig. 4a) with the greatest amount (1.1 million ha; 63%) of expected moderate to high burn severity, followed by the Lower Green (Fig. 4b; 689,000 ha; 33%), the Upper Green (Fig. 4c; 645,000 ha; 20%), and the Dirty Devil (Fig. 4d) with the least (144,000 ha; 29%). Predicted dNBR ranged from a low of –0.11 to a high of 0.57 (Fig. 3).

Wildfire models were unique in terms of variable selection. In total, there were 329 instances of variable inclusion across 28 wildfire-specific models, with ecosystem factors accounting for 140 instances of variable inclusion, topography 73, weather 71, and fuel 45 (Fig. 6). Accounting for variable overlap among models resulted in 76 unique variables used, averaging 12 variables per model. Controlling for period and lag resulted in 42 variables accounting for pre-fire predicted burn severity (Fig. 6). The most frequently used (25 of 28 models) predictor among all models was of HERB whereas the predictor with the greatest SPVI was elevation (19/28 models, Fig. 6). The most frequently selected predictors (in 14/28 fires or more) included both variables related to fuels (HERB and FCC), two variables related to topography (elevation and terrain position index), and three remotely sensed variables related to physical ecosystem (Sentinel-2 bands 9 (945 nm) & 12 (2190 nm), and SAVI; Fig. 6). Both Sentinel-2 bands 9 and 12 measure reflectance of short-wave infra-red radiation, useful for detecting water vapour and vegetation moisture stress, respectively (European Space Agency 2023). Based on SPVI, the most influential predictors included elevation and three variables related to ecosystem condition (Sentinel-2 band 12, Sentinel-2 band 9, and leaf area index; Fig. 6). The most frequently included weather variable (7/28 models) was the amount of precipitation with a 7-day lag. When SPVI was averaged by environmental factors the largest contributions were from fuel (1.20), followed by topography (0.83), physical ecosystem (0.81), and weather (0.28).

Fig. 6.

Factors, total predictor variable counts over all models, and percentile of Summed Proportional Variable Importance (SPVI), with period and lags removed, of random forest models used to predict burns severity from 28 wildfires that occurred between 2016 and 2020. Topography abbreviations: TPI, terrain position index; TRI, terrain roughness index; RUF, surface roughness; COS, cosine of aspect; SIN, sine of aspect. Fuel abbreviations: HERB, herbaceous cover (HERB); FFC, forest canopy cover. Ecosystem abbreviations: Sentinel-2 Bands (2 = B2, 3 = B3, 4 = B4, 5 = B5, 6 = B6, 7 = B7, 8 = B8, 8a = B8a, 9 = B9, 11 = B11, 12 = B12); NDVI, normalised difference in vegetation index; SAVI, soil-adjusted vegetation index; LAI, leaf area index; FAPAR, fraction of absorbed photolytically active radiation; CCL, chlorophyll content in leaf; FVC, fraction of vegetation cover; CWC, canopy water content. Weather abbreviations: FM100, 100-h dead fuel moisture; FM1000, 1000-h dead fuel moisture; TMMX, maximum air temperature; TMMN, minimum air temperature; PR, precipitation; SRAD, downward surface shortwave radiation; WS, windspeed; RMAX, maximum relative humidity; RMIN, minimum relative humidity; SPH, specific humidity; ETR, alfalfa reference evapotranspiration; PET, grass reference evapotranspiration; ERC, energy release component; BI, burning index; VPD, vapour pressure deficit.

Potential for post-fire debris flow

The overall number and lengths of stream reaches of bCRCT CP (1903) analysed varied by HUC-4 watershed boundary: 727 stream reaches in the Upper Green (1067 km), 604 in the Yampa (718 km), 496 in the Lower Green (648 km), and 76 in the Dirty Devil (89 km). Excluding null values of predicted debris-flow reduced the number of bCRCT CP stream reaches with classified debris-flow of either low, medium, or high to 1643. By HUC-4 watershed boundary this resulted in 638 of 727 bCRCT stream reaches (88%) in the Upper Green, 532 of 604 in the Yampa (88%), 409 of 496 in the Lower Green (82%), and 64 of 76 in the Dirty Devil (84%) with an elevated hazard classification from PFD.

The combined hazard classifications of potential PFDF of stream reaches currently occupied by bCRCT CP increased with storm intensity (Figs. 7, 8). Overall, 86% of streams reaches with bCRCT CP exhibited elevated hazard classification of potential PFDF in response to increased storm intensity (Table 4, Fig. 7). At maximum storm intensity of 40 mm h⁻¹, potential PFDF increased to the highest hazard classification in roughly a third of all stream reaches (Table 4 and Fig. 7). Above 28 mm h⁻¹, more than two-thirds of stream reaches showed potential PFDF hazards of moderate or high, increasing up to 80% at 40 mm h⁻¹ (Table 4, Fig. 8). Combined hazard classifications of the moderate category were present throughout all storm intensities; above 24 mm h⁻¹ more than half of stream reaches were classified as such (Table 4). Potential PFDF over all stream reaches of bCRCT CP showed elevated moderate and high hazard as peak precipitation intensities increased. Similar patterns were present when tabulated by HUC-4 watershed (Table 4).

Fig. 7.

Potential post-fire debris flow (PFDF) combined (likelihood, %; volume, m³) hazard classifications based on a simulated peak 15-min rainfall intensity of 40 mm h⁻¹ for stream reaches with blue lineage Colorado River Cutthroat Trout (bCRCT; Oncorhynchus clarkii pleuriticus) conservation populations and nodes. Note: inset maps show reference for scale and examples of mainstem outlets used for aggregation to node.

Fig. 8.

Proportion of the combined (likelihood, %; volume, m³) hazard classifications (high, and moderate or high) of potential post-fire debris flow with blue lineage Colorado River Cutthroat Trout (Oncorhynchus clarkii pleuriticus) Conservation Populations (bCRCT CP) and their associated stream reaches by 4-digit Hydrologic Unit Code (HUC-4) watershed boundaries based on predicted burn severity of wildfire in unburned areas and simulated peak 15-min rainfall intensities (mm h⁻¹).

The number of bCRCT CP populations within high PFDF hazard potential varied greatly among basins. The Yampa had the greatest number of stream reaches (343) occupied by bCRCT CP in the highest hazard classification at maximum storm intensity, followed by the Upper Green (135), the Lower Green (131), and the Dirty Devil (14). Proportionally, the Yampa also had the greatest proportion of occupied stream reaches in the highest hazard classification at maximum storm intensity (0.57), followed by the Lower Green (0.26), the Upper Green (0.19), and the Dirty Devil (0.18; Table 4). By combined stream lengths, the Yampa had the greatest distances (480 km) in the highest hazard classification, followed by the Lower Green (245 km), the Upper Green (214 km), and the Dirty Devil (12 km). When both high and moderate hazard classifications were considered, the rankings among watersheds shifted subtly. Specifically, the Upper Green contained the greatest number of stream reaches (584; 80%) occupied by CP classified as moderate or high followed by the Yampa (516; 85%) indicating there were fewer stream reaches, but a higher percentage, in the Yampa classified as high or moderate then in the Upper Green. Analysis of the Lower Green showed 372, or 75%, (Fig. 8) of stream reaches with moderate or high combined hazard classification, and the Dirty Devil with 51, or 67% (Fig. 8).

Conservation population hazard

In total, CP of bCRCT spanned ~2500 km within 422 spatially distinct GNIS streams comprised of 1903 reaches that constituted the 190 nodes (Fig. 7): 66 in the Yampa; 65 in the Upper Green; 37 in the Lower Green; and 22 in the Dirty Devil. (Fig. 7). Within 500-m buffers of streams with bCRCT CP average predicted dNBR (Fig. 3) ranged from 0.17 to 0.46 (see Supplementary Table S1). When tabulated by HUC-4 watershed, the proportion of predicted moderate or high burn severity classifications (Fig. 4) within 500-m buffers was highest in the Yampa watershed (0.76), followed by the Upper Green (0.47), Lower Green (0.43), and the Dirty Devil (0.41). Overall, proportions of moderate or high severity were greater within 500-m buffers around streams occupied by bCRCT CP than across HUC-4 watershed.

The mainstem outlet points averaged over stream reaches of associated bCRCT CP nodes (Fig. 7) produced hazard index values per simulated storm intensity (Table S2). Hazard index values showed a steady increase in hazard to bCRCT CP from PFDF classifications as storm intensity increased (Table S2). At peak storm intensity (40 mm h⁻¹) 184 of 190 bCRCT CP (Fig. 8; 97%) were at risk from moderate (52%) or high (45%) potential PFDF hazards. Reflecting overall trends in PFDF hazard classifications by stream reach (Fig. 8), the bCRCT CP accumulated high hazard classifications at peak storm intensities exceeding 24 mm h⁻¹ (Fig. 8). At the lowest storm intensity (12 mm h⁻¹) 16% of bCRCT CP were at risk from moderate potential PFDF, with the remaining 84% classified as low. At the finest spatial scale of analysis, that of the bCRCT conservation population node (Fig. 7), cumulative combined hazard of potential PFDF was greatest, relative to stream reach and watershed.

Predicting burn severity

Our models of dNBR had an overall average adjusted R² = 0.54, which is within the range (0.26–0.63) of values reported for other models (Birch et al. 2015; Kane et al. 2015; Picotte et al. 2021; Pascolini-Campbell et al. 2022; Wells et al. 2023). Our hindcasting of one model predicting the outcome of a different wildfire resulted in a lower R² value (0.25) but provided confirmation of general spatial patterning. Use of Relativised dNBR in future work may improve model fit as we expect a high degree of heterogeneity in vegetation, and therefore spectral reflectance, of pre-fire factors (Cansler and McKenzie 2012), however this would necessitate recalibration of the PFDF equations. Our study represents an advance in understanding potential burn severity before a wildfire occurs by employing pre-fire conditions of fuel, topography, physical ecosystem, and weather. Much of the unexplained variation between predicted and observed dNBR and burn severity is likely due to the inherently dynamic behaviour of wildfire that can rapidly change with prevailing weather and fuel conditions. Incorporating diurnal fluctuations (Poulos et al. 2021; Pascolini-Campbell et al. 2022) in ambient air temperature, humidity, and wind speed that are specific to the location of the wildfire may improve future predictions.

The threat of wildfire was high across HUC-4 watersheds. The Yampa watershed showed the greatest areal extent of expected moderate to high burn-severity, which influence PFDF. Recent wildfires, such as the East Troublesome Fire within similar ecosystems and fuels conditions (beetle-killed conifer) have exhibited a rapid rate of spread and high burn severity (Meldrum et al. 2022). The spruce/fir vegetation type existing at upper elevations within the Yampa basin has a low fire return interval of ~300 years (Schoennagel et al. 2004) and high predicted burn severity. Similarly, high predicted burn severity was reflected in upper elevations of the Upper Green (Fig. 4a) and Lower Green (Fig. 4b) HUC-4 watersheds, with forested assemblages of Douglas fir/spruce/fir. The lower, more arid elevations of the Upper and Lower Green primarily support sagebrush (Artemisia spp.) and pinyon/juniper (Pinus edulis, Juniperus spp.) vegetation, respectively. Consequently, we expect to have lower predicted burn severity primarily due to the lower amount of woody biomass and the greater spatial heterogeneity of these fuel types. This pattern was mirrored in the arid and southerly Dirty Devil watershed, largely composed of pinyon/juniper but with forested Douglas fir assemblages occurring at higher elevations. These relationships to severity were supported by the frequency and importance of elevation and FCC variables. We expected HERB, which was frequently included in models, to influence burn severity through its influence on rate of fire spread and intensity. Despite its importance, fine fuel is often not well characterised in existing fire models and can be improved upon with future efforts (Wells et al. 2021). Both environmental factors of fuel and topography were readily apparent in the spatial distribution of predicted burn severity. Large swaths of the Upper and Lower Green and Dirty Devil joined the Yampa watersheds with high burn severity classification, comprising roughly 1.1 million ha, or 15% of the study area. Combined with moderate burn severity, those classes used to estimate PFDF, covered roughly 2.6 million ha of the study area. Future research should investigate the differences in variables used to model pre-fire burn severity based on plant community (i.e. lower elevation shrublands vs higher elevation forests) to better understand the drivers of burn severity across different ecotypes.

Potential for post-fire debris flow

We took a novel approach to modelling potential PFDF at the stream reach scale by including predicted burn severity before a wildfire occurs (Cannon et al. 2010; Sedell et al. 2015). Current operational USGS PFDF models (U.S. Geological Survey 2023) rely on burn severity estimated post-fire. Prediction of burn severity using our methodology represents an important step in developing and refining hazard classifications for broad application to unburned but potentially at-risk areas. Additionally, our PFDF estimates reflect landscape factors that likely initiate, propagate, and accumulate debris-flow based on basin morphology, slope, soils, and precipitation. Adjustment of soil parameters including hydrophobicity, the loss of infiltration, and other physical changes due to fire effects (Lopes et al. 2021; Ebel et al. 2023) were not considered in the model but offer opportunities for further study. Likewise, additional research focused on sediment delivery downstream of watershed outlets may improve hazard assessment, as this physical process is not accounted for by the USGS operational model. The application of the debris flow assessment prior to wildfires and across large spatial extents has been historically restricted by the requirement of post-fire dNBR and burn severity estimates. Here, we have presented an opportunity to shift the assessment framework from a post-fire emergency reaction to a pre-fire early warning, enabling wide-spread assessment. This can increase the robustness of the emergency assessment system as well as proactively identifying the risk to streams, aquatic species, or other water resources before a real hazard emerges.

The potential PFDF combined hazard over all stream reaches of bCRCT CP was almost entirely moderate to high for high intensity rainfall events. The hazard in the Yampa was more pronounced based on the greater amount of moderate to high predicted burn severity, but the hazard in the Upper Green was collectively broader based on the greater number of CRCT stream reaches with moderate to high potential PFDF. Hazard was similar between Upper Green and Lower Green in terms of number and length of streams in the highest hazard class. Spatially, the Lower Green had a similar overall length of stream reaches of bCRCT populations as that of the Yampa, but the Yampa was proportionally much greater in terms of overall number and length of streams with the highest potential PFDF. The Dirty Devil watershed showed the lowest potential PFDF by number and length of stream reaches with bCRCT populations. Risk in Dirty Devil should not be disregarded on the basis of lower hazard classifications however, as small population size and structure may increase the vulnerability to localised extinction (Cook et al. 2010).

Past research to assess risk to Cutthroat Trout populations from wildfire (Williams et al. 2009) or PFDF (Sedell et al. 2015) or wildfires and PFDFs (Miller and Bassett 2013) have used a variety of GIS (Williams et al. 2009; Sedell et al. 2015) and statistical tools (Cannon et al. 2010; Miller and Bassett 2013) to identify hazards. This includes analysis based on landscape characteristics derived from elevation models (10-m DEM; Sedell et al. 2015), in conjunction with LANDFIRE data products (Williams et al. 2009) and FlamMap analysis (Miller and Bassett 2013) at the stream reach (Sedell et al. 2015) or basin (Williams et al. 2009) scale. Our framework melds a singular approach for estimating of hazards derived from both wildfire (Wells et al. 2023) and debris flow (U.S. Geological Survey 2023) for both basin and stream reach scales. While this was an important modelling improvement, the emerging effects of climate change increase the need for an early warning system integrated with near-term climate forecasting (Bradford et al. 2018). Similarly, validation efforts to describe model fitting of PFDF events will benefit from further research to detect PFDF based on hydrologic data from widespread monitoring. Our efforts to validate potential PFDF classifications were hindered by the sparsity of data on PFDF events that may have occurred originating from within a study wildfire. Improvements to PFDF predictions can be advanced through further research that detects (Hürlimann et al. 2019) and identifies the impact of PFDF on hydrologic function.

Conservation population hazard

For bCRCT CP, estimates of the combined hazard classification of potential PFDF were cumulatively compiled by stream reaches and outlets to quantify hazard (Table S2), similar to past work assessing vulnerability of CRCT habitat (Conservation Science Partners [Theobald D, Leinwand I, Harrison-Atlas D] 2018). Hazard of PFDF to bCRCT CP was greatest for the Yampa, followed by the Upper Green, Lower Green, and Dirty Devil watersheds. The higher hazard to bCRCT CP in the Yampa basin was a result of the greatest number, length, and proportion of moderate or high classifications of PFDF over all simulated storm intensities. Hazard indices can be included in population viability assessments that predict population persistence in a changing climate (e.g. Zeigler et al. 2019), and more importantly, enable CRCT managers to consider risk from wildfire to individual habitat patches defined by Hirsch et al. (2013) and develop mitigation strategies accordingly. Management actions that could help secure native trout diversity in a more fire-prone future could include founding new CP with translocations of CRCT from high-risk patches or important peripheral populations (Haak et al. 2010b), reducing burn severity by defining and implementing strategic fuel reduction opportunities, reducing the chances of ignition (e.g. public access restrictions or powerline blackouts), and engineered debris-flow mitigation (retention structures e.g. Volkwein et al. 2011). Most importantly, the ability to predict PFDF hazard will help inform pre-fire planning to determine whether emergency fish salvage operations should proceed when a fire threatens or has burned through a CP. These efforts involve risk to personnel as well as resident trout and are therefore generally missions of last resort but have been used with good success to save valuable native trout populations, such as those at Hayden Creek, CO (Rogers et al. 2020).

Our approach provides a quantifiable characterisation of hazard based on the long-term indirect effects of wildfire. These results build upon past work to forecast burn severity and subsequent PFDF before they happen, helping to inform conservation planning and long-term risk management of valued fisheries. By using an explicit spatial evaluation and assessment, we are able to bridge an implicit understanding of burn severity from both direct and indirect effects of wildfire to an at-risk aquatic species in western North America.