Using wildfire as a management strategy to restore resiliency to ponderosa pine forests in the southwestern United States

Abstract

The long‐term outcome from accelerated forest restoration using resource objective wildfire in combination with fuel management on fire‐excluded landscapes is not well studied. We used simulation modeling to examine long‐term trade‐offs and synergies of alternative land management strategies by combining two wildfire management alternatives with three levels of contemporary forest restoration treatments on a 778,000‐ha landscape over 56 years. We found that managing wildfires for resource objectives diminished the likelihood of irregular fire events over time by making wildfire activity more predictable. Overall, adding resource objective wildfire reduced the proportion of high‐severity fire in relation to total area burned, but increased total area burned and the area of high‐severity fire. We also found resource objective wildfire changed the distribution of high‐severity burn patches by increasing their total number and range, their likelihood of containing disjunct core areas, and their edge complexity. The results suggested that alongside the current pace of active forest management, expanding the fire footprint to achieve low‐cost restoration carries the potential for increased high‐severity fire and associated impacts to ecological values including old forest structure and wildlife habitat. Concurrently, adding resource objective wildfire expanded the footprint of conventional restoration treatments by fivefold, and restoration objectives were achieved in 25 years when managing resource objective wildfires alongside restoration treatments five times the current pace. This study demonstrates the first fire suppression model to replicate local decision making by fire managers during simulated fire events to manage risk by considering both fire proximity to values at risk, and daily weather conditions. The study paves the way for further investigations of the synergies between wildfires and conventional forest restoration to improve resiliency in fire‐excluded pine forests.

INTRODUCTION

The restoration of frequent-fire pine forests in the western United States has faced multiple challenges over the past several decades, and there is a growing sense of urgency given recent projections of a warming climate. To expand and accelerate restoration efforts, land managers have been augmenting contemporary forest restoration treatments such as thinning, mastication, and prescribed fire with the practice of using natural ignitions to meet multiple ecological and management objectives (henceforth resource objective [RO] wildfire) (Huffman et al., 2018; Hunter et al., 2014). Resource objective wildfire is rooted as far back as 1968 in US National Parks, and has been extensively practiced in US Forest Service wilderness areas (Van Wagtendonk, 2007). Over its history, RO wildfire has had many identities including the policy of “Let it burn,” and descriptive terms such as prescribed natural fire, wildland fire use, managed fire, or restoration wildfire (Barros et al., 2018; Hunter, 2007). Following 2009, the tactics and strategies used in the application of RO wildfire expanded after changes in national fire policy encouraged the practice (Young et al., 2020), and local policies further supported it (USDA, 2014, 2018). However, the process of transforming fire-excluded landscapes to conditions where RO fire is systematically used in concert with traditional forest management is not well defined. Much of the recent experience with RO wildfire is in frequent-fire ecosystems of the southwestern United States (e.g., Huffman et al., 2018) when fuel conditions allow a wildfire to be carried out without posing a significant risk to vital habitat or infrastructure (Young et al., 2019).

Wider application in fire-adapted ecosystems throughout the western United States presents multiple challenges to land management agencies. Though RO wildfires have the potential to efficiently reduce hazardous fuels and future per hectare suppression costs (Gebert & Black, 2012; Houtman et al., 2013; North et al., 2012), especially where mechanical thinning is not feasible (Doane et al., 2006; Huffman et al., 2020; North et al., 2015), specific conditions are required including discontinuous fuel beds to facilitate fire control (Hunter et al., 2011), and remote locations (Miller & Aplet, 2016) where wildfire is less likely to damage valued resources in the wildland–urban interface (WUI; Scott et al., 2012). Typical western US fire-excluded landscapes with high loads of hazardous fuels are not well suited for either RO or prescribed fire due to the acute risk of negative fire effects, where large patches with high burn severity can catalyze transitions to grass and shrub types (Kane et al., 2013), including non-natives (Kerns et al., 2020; Stevens-Rumann et al., 2016). Near the WUI, forest restoration is more consistently achieved with mechanical thinning accompanied by prescribed fire. This combined technique commonly has a high implementation cost (Hjerpe et al., 2009; Nicholls et al., 2018) but it also reduces wildfire suppression costs in frequent-fire forests (Fitch et al., 2018) while potentially reducing direct and indirect wildfire damage (Combrink et al., 2013; Combrink & Rousse, 2018), for example, from post-wildfire debris flows (Staley et al., 2018).

Land managers in the southwestern United States are increasingly incorporating RO fire management with contemporary restoration treatments, seeking to harness the benefits of both strategies (Young et al., 2020). Despite the ongoing use of RO wildfire throughout the southwestern United States (Huffman et al., 2017, 2018; Hunter et al., 2011, 2014) with numerous successes (e.g., 2008 Marteen Fire, 4366 ha; 2009 Wildhorse Fire, 5581 ha; and 2014 Sitgreaves Fire, 4350 ha), the long-term effects of RO wildfire in combination with conventional restoration are unknown, and escaped fires have caused adverse outcomes (e.g., 2006 Warm Fire, 23,727 ha). Questions remain about how to spatially prioritize landscapes for conventional treatments to optimally prepare them for future RO fire (Barros et al., 2018). For instance, can contemporary restoration treatments pave the way for expanding stochastic RO fire to efficiently create resilient forests with minimal risk to social and ecological values?

In this study, we used a forest landscape model to examine long-term trade-offs and synergies of alternative land management strategies by combining two wildfire management alternatives with three levels of contemporary forest restoration treatments. Principally, we examined the effect of expanded RO wildfire on a frequent-fire landscape in the southwestern United States. The model employed a novel fire suppression submodule that simulated local decision making by fire managers during fire events, considering both fire proximity to values at risk, and daily weather conditions that allow for managing risk from RO fire events. We hypothesized that RO wildfire will improve the resiliency of the ponderosa pine ecosystem of the southwestern United States by (1) reducing the relative occurrence of high-severity wildfire, (2) enhancing wildlife habitat and foraging opportunities for goshawk and Mexican spotted owls by enhancing forest heterogeneity and old-growth forest structure, and (3) reducing the need for contemporary restoration treatments. We measured resiliency with the presence of old-growth forest structures and overall potential for high-severity fire and associated patch dynamics, recognizing that the latter can have a significant deleterious effect on forest regeneration in southwest pine ecosystems and elsewhere (Owen et al., 2017, 2020; Singleton et al., 2021; Stevens-Rumann et al., 2016; Stevens-Rumann & Morgan, 2016). The effects of climate change on fire regimes and forest ecology, which have been found to be limited in northern Arizona when compared to land management actions (O'Donnell et al., 2018), are not examined in this study. We contribute to understanding the long-term implications of RO wildfire, including its interaction with other restoration treatments to affect fire severity, wildlife habitat, and patch dynamics of high-severity fire.

METHODS Study area

The 778,378-ha study area used for analysis sits along the Mogollon Rim on the southwestern extent of the Colorado Plateau in northern Arizona, inclusive of the Coconino Plateau along the northern extent (Figure 1). The larger model area used for wildfire simulations is comprised of 1,488,900 ha of public, private, and tribal lands spanning from deserts, grasslands, and pinyon–juniper vegetation types at lower ecotones, to ponderosa pine forests in the midrange ecotones and mixed conifer in the upper ecotones. Much of the wooded and forested lands in the modeling area are within the study area and are managed by the US Forest Service (USFS). The Forest Service is currently restoring 237,218 ha of ponderosa pine forest in the study area under the Four Forest Restoration Initiative (FFRI; USDA, 2015a) using restoration treatments of mechanical thinning and prescribed fire. The objective is to increase the health and function of the forest by returning old-growth forest structure to the landscape, making it more resilient to wildfire. The study area and wider southwestern United States have a history of using RO wildfire to augment restoration treatments (Young et al., 2019), with increased fuel reduction following multiple wildfires that burned at low to moderate severity (Huffman et al., 2017, 2018; Hunter et al., 2011, 2014). The expanded use of RO wildfire has also shown the potential to promote floristic diversity (Laughlin et al., 2005). Despite local success, the task of restoring the landscape is a vast undertaking complicated by an extensive backlog of fuel reduction needs (Marlon et al., 2012; North et al., 2015) in an area with a recent history of high-severity wildfires (e.g., 2000 Pumpkin Fire, 3545 ha; 2010 Shultz Fire, 6101 ha) that contained large patches of sporadic and inconsistent tree regeneration (e.g., Owen et al., 2017).

FIGURE 1. (a) Model area in northern Arizona, study area on the Kaibab and Coconino National Forests, and active management actions described in Table 1. The model area is used to simulate wildfires, and the study area is used to simulate treatments in 38 delineated planning areas. Planning areas are prioritized based on potential fire mortality as a percentage of total basal area as described in section “Forest and wildfire management scenarios.” DMT, dwarf mistletoe; WUI, wildland–urban interface. (b) Boundaries for wildfire management objectives. Protection objectives are located around WUI, mixed-conifer and spruce–fir ecoregions, national monuments, and steep interior chaparral canyons leading to the wildland–urban interface. Resource objectives are located within ponderosa pine and pinyon–juniper forests

Forest Service lands within the study area where restoration treatments are allowed according to FFRI (henceforth manageable lands) include 237,218 ha across the Kaibab and Coconino National Forests. The study area has 26,944 tree stands over 778,378 ha that have been delineated into 38 planning areas (PAs) used to prioritize landscape treatments (Figure 1). Manageable lands in each PA ranged from 1113 to 12,315 ha of ponderosa pine forest. In total, the PAs are 49% ponderosa pine forests and savannas with a 3% median slope and altitude ranging from 1657 to 2828 m above sea level (asl). Mixed-conifer, aspen, and spruce–fir forests comprise 2% of the PAs with a 17% median slope and altitude ranging from 1576 to 3783 m asl. Pinyon–juniper ecological types comprise 26% with a median slope of 3%, and an altitude ranging from 1019 to 2756 m asl.

Local restoration management objectives were obtained from the FFRI management team and are centered on restoring goshawk and Mexican spotted owl habitat, reducing widespread dwarf mistletoe infestation, restoring ponderosa pine savannas and grasslands to presettlement conditions, and reducing hazardous fuels with prescribed fire alone (Figure 1a, Table 1) (USDA, 2015a). The study area also contains 7425 ha of urban development or WUI in need of protection from wildfire. This classification was determined by WUI designation in the SILVIS Lab's 2010 WUI database (Radeloff et al., 2017) and was then expanded based on WUI data provided by the FFRI management team of the USFS (Figure 1a).

TABLE 1 Treatment type, percentage of study area, and associated active management actions

Treatment type	Acronym	Percentage of study areaa	Active management action
Goshawk habitat restoration	Goshawk	32%	Thin throughout to SDI = 130 followed by prescribed fire
Prescribed fire only	Fire	27%	Prescribed fire
Mexican spotted owl habitat restoration	Owl	14%	Thin throughout to SDI = 115 followed by prescribed fire
Dwarf mistletoe mitigation	DMT	10%	Thin from below to SDI = 135 followed by prescribed fire
Ponderosa pine grassland restoration	Grassland	10%	Thin to 5% canopy closure followed by prescribed fire
Ponderosa pine savanna restoration	Savanna	7%	Thin throughout to SDI = 14 followed by prescribed fire

Note: Each tree stand is assigned an active management action that in addition to the action described could remove no more than half of the ponderosa pine with a diameter at breast height greater than 40.6 cm in order to retain old-growth forest structure (FFRI, 2011). Acronyms are referenced in Figure 1. Abbreviations: DMT, dwarf mistletoe; SDI, stand density index; WUI, wildland–urban interface. a Study area = 237,218 ha. Overview of LSim

LSim is a forest landscape management model that was built by linking the Parallel Processing version of the Forest Vegetation Simulator (FVS; Crookston & Dixon, 2005) with the large fire simulator (FSim; Finney et al., 2011). Both models are extensively described in decades of foundational research, and their coupling into LSim and application to landscape forest and fire management was described in Ager et al. (2020). The modeling components are described in brief below with detailed descriptions of the model included in Appendices S1 through S6. Further modifications to LSim for this specific study are described below (Figures 2 and 3).

FIGURE 2. Diagram of landscape simulator (LSim) model components and modifications to simulate resource objective wildfire (black boxes and circles). FDIST, fire distribution file used to guide the daily occurrence of wildfire ignitions based on the energy release component (ERC); FRISK, fire risk file used to model wildfire occurrence and spread based on fire weather metrics including percentile fuel moisture and ERC; FSim, large fire simulator. Figure modified from Ager et al. (2020)

FIGURE 3. Detailed diagram of the landscape simulator (LSim) program components showing the functionality in the two main submodules, the Forest Vegetation Simulator (FVS) and the large fire simulator (FSim); and the data translator within the parent LSim model. For each simulation cycle, FVS loops through each stand of the landscape to simulate fire mortality, growth, management, and natural mortality. The resulting landscape fuels are translated to a binary raster stack required by FSim. FSim then loops through each day in the fire season and simulates fires based on daily predicted weather. At the end of the fire season, the wildfire intensity grid is read by FVS and wildfire-induced tree mortality is predicted using the functionality of the FVS-Fire and Fuels Extension. ERC, energy release component used in wildfire simulations; FDIST, fire distribution file used to guide the daily occurrence of wildfire ignitions based on ERC; FRISK, fire risk file used to model wildfire occurrence and spread based on fire weather metrics including percentile fuel moisture and ERC; MTT, minimal travel time of wildfire spread. Figure modified from Ager et al. (2020)

Forest vegetation modeling system

The Forest Vegetation Simulator (FVS; Crookston & Dixon, 2005) is used by LSim to model forest succession and management. FVS is a nonspatial, individual tree growth and yield model used to simulate stand development and mortality based on empirical relationships. The simulation unit is a homogenous polygon of trees (i.e., tree stand) described by a record of tree density by species and size, along with elevation, aspect, and ecotype (Crookston & Dixon, 2005). There are over 20 variants that cover forest types in the United States, including the Central Rockies (CR) variant that was used in this study (Keyser & Dixon, 2008).

LSim integrated the FVS-Parallel Process Extension (FVS-PPE; described in Ager et al. [2020] and Seli et al. [2008]) created by Crookston and Stage (1991) to adapt FVS to landscape-scale simulations. FVS-PPE pauses the simulation of forest dynamics (growth, regeneration, and management) at user-defined time intervals (e.g., 1 year) and writes selected stand metrics to an external file. Exogenous programs can perform calculations to reprioritize management activities, apply landscape-scale constraints, and generate input files for wildfire and other simulation programs. Prior research with FVS-PPE and the LSim model include studies on wildlife habitat, carbon dynamics, risk to WUI, the effect of accelerated management scenarios (Ager et al., 2007, 2020; Ager, Vaillant, & Finney, 2010), and the spatial optimization of fuel treatments (Finney, 2007).

Wildfire modeling

Wildfires were simulated with the large fire simulator (FSim; Finney et al., 2011) (Figure 3) on a 90-m grid covering the model area, which is inclusive of a 5- to 35-km buffer around the study area to reduce edge effects of the analysis. Due to the local geography (e.g., Grand Canyon National Park to the north) and sparse vegetation surrounding the study area, edge effects would be limited to neighboring ponderosa pine forests in the southeastern extent of the study area. To simulate wildfires, FSim iterates through each day of the calendar year to model fire ignition and spread (Figure 3). Fire ignitions are determined based on the historic distribution of fire ignitions (FDIST file) as a function of the energy release component (ERC via the fire risk (FRISK) file) in combination with the distribution of multiple large fire days. If an ignition occurs, data in the binary raster stack containing landscape fuels (surface fuel model, canopy cover, canopy height, crown base height, and crown bulk density) and topography data (elevation, slope, and aspect) (Appendix S1: Figure S1) in combination with fire weather conditions (FRISK file) are used to determine fire spread. The fire continues to spread until contained as a function of the wildfire's rate of spread and the occurrence of timber fuels, where slower spreading fires in non-timber fuels are more likely to be contained (Finney et al., 2009). In our case, additional perimeter trimming was used to simulate progressive fire perimeter containment to limit fire growth on portions of the fire's perimeter with low fireline intensity (Scott et al., 2018). Perimeter trimming can be adjusted to assist in the calibration of fire models seeking to emulate historical fire activity. FSim output files used by LSim each cycle included a perimeter shapefile and a fireline intensity grid that was converted to flame length using Byram's equation (Byram, 1959) with the Wilson (1980) modification. Each raster pixel of the model area was burned by one wildfire within each cycle (i.e., FSim switch BurnCellsOnlyOnce was engaged).

Integrating FVS and FSim

FVS and FSim were integrated into LSim using a wrapper that sequences their execution and transfers data each cycle (Figure 2). We chose to use 4-year cycles rather than a shorter time frame to enhance model performance associated with constructing fuel landscape files (i.e., binary raster stack). After each cycle, LSim pauses vegetation modeling (i.e., FVS) to translate forest inventory data into a binary raster stack required by FSim (Appendix S1: Figure S1), which then simulates wildfires for four seasons (Figure 3). A raster grid of wildfire perimeters is then transferred back to FVS where it is overlaid with the stand polygon layer. If more than 50% of a stand's pixels burn during a cycle, the stand is determined to have burned. For each burned stand, the fire line intensity grid from FSim is converted to flame length and used to calculate average scorch height using Byram's equation (Byram, 1959) with the Wilson (1980) modification. The average flame length is then used alongside the tree stand's average crown base height to compute the percent of the stand that burned in a crown fire (Ager et al., 2020, Appendix S1). Together, flame length, scorch height, and percent crown fire serve as parameters in the FLAMEADJust keyword via the FVS-Fire and Fuels Extension (FVS-FFE; Rebain et al., 2010). A wildfire is then simulated with the SIMFIRE keyword, and tree mortality is predicted with the first-order fire effects model (Reinhardt et al., 1997). Tree mortality is then used to adjust the forest inventory before FVS resumes for the next cycle (Figures 2 and 3). Additional details on wildfire-induced tree mortality are described in Appendix S1. This process is similar to methods in Ager et al. (2010, 2020).

Model parameters Forest and fuel management

We developed stand treatments (Table 1) based on management objectives developed by the FFRI management team of the USFS (USDA, 2015a) and implemented them using FVS keywords (Appendix S2: Figure S1). One consistent management objective is to restore low-severity fire regimes by returning old-growth forest structure to a large portion of the ponderosa pine forest (FFRI, 2011; USDA, 2014, 2015a, 2018). Local management objectives are described in Table 1. Each treatment is supported by empirical studies as effective for reducing potential fire behavior (e.g., Fulé et al., 2012; Parks et al., 2018; Stephens et al., 2009). Simulated treatment actions included restoration thinning followed by prescribed fire, precommercial thinning followed by prescribed fire, or prescribed fire alone (Table 1; Appendix S2: Figure S1). LSim determines potential treatment actions based on habitat-specific objectives, thresholds of merchantable timber, trees per hectare, time since fire, and fuel loads. Prescribed fire fuel moistures and wind speed were chosen to replicate typical prescribed burning conditions in both the spring and the fall on the Coconino and Kaibab National Forests (Appendix S2: Table S1) including a blend of fixed (Appendix S2: Table S2) and predicted mortality via FVS-FFE (Rebain et al., 2010).

Forest growth and mortality

The CR variant was built using data from forest inventories, research plots, and silviculture stand assessments to derive growth relationships for 38 tree species or species groups (Keyser & Dixon, 2008). The tree list data populating our stands were provided by the FFRI management team of the USFS (USDA, 2015a) and then supplemented by a forest inventory analysis (FIA)-imputed tree list developed by Riley et al. (2016). Recent mechanical thinning (2009–2019) and wildfire activity (2009–2017) data from the USFS Accomplished Activities database (USDA, 2020a) and Monitory Trends in Burn Severity (MTBS, 2018) were used to adjust stand tree lists to reflect disturbances. The landscape data were then vetted further by land managers of the Coconino and Kaibab National Forests to ensure species were properly distributed.

Non-wildfire-related tree mortality is predicted using species-specific, stand density index (SDI)-based mortality models as determined by the relative SDI of the stand. When SDI is <55% of the maximum, mortality is predicted based on species-specific models, and when SDI is ≥55% of maximum, density-dependent mortality is predicted for the stand (Keyser & Dixon, 2008). Both are applied with greater effect in areas with high basal area and are adjusted for species-specific shade tolerance (Keyser & Dixon, 2008). Wildfire-related tree mortality is described in Appendix S1 using the same methods as Ager et al. (2010, 2020). Refer to Appendix S2 and Keyser and Dixon (2008) for additional details on forest growth modeling and tree mortality.

Forest regeneration

Regeneration models used by LSim for our study area were developed with data from Coconino County of northern Arizona (Young, 2019) and then calibrated based on prior research (e.g., Singleton et al., 2021) and expert opinion (Appendix S3). In this study, we use regeneration models for nine tree species (blue spruce [Picea pungens], corkbark fir [Abies lasiocarpa var. arizonica], Douglas-fir [Pseudotsuga menziesii], Engelmann spruce [Picea engelmannii], Gambel oak [Quercus gambelii], limber pine [Pinus flexilis var. flexilis], ponderosa pine [Pinus ponderosa], quaking aspen [Populus tremuloides], and white fir [Abies concolor]) and three species groups (pinyon pine [e.g., common needle pinyon (Pinus edulis)], sprouting juniper [e.g., alligator juniper (Juniperus deppeana)], and non-sprouting juniper [e.g., Utah juniper (Juniperus osteosperma)]) across three ecological zones: pinyon–juniper woodlands, ponderosa pine transition, and mixed-conifer forests. Data for the pinyon–juniper zone were collected by Huffman et al. (2012), while data for other zones were collected by the FIA program (USDA, 2020b). With these data, we developed regeneration models using zero-inflated (independent modeling of regeneration odds and number of stems) random forest and count modeling procedures as a function of topographic (elevation, aspect, and slope) and ecological factors (species-specific trees per acre and basal area metrics capturing composition and structure). Refer to Appendix S3 for empirical regeneration models. To further account for the probability of random failures of regeneration in undisturbed forests due to seed predation, frost heaving, etc., each tree species is calibrated to randomly engage regeneration models based on expert opinion (e.g., undisturbed stands will engage ponderosa pine regeneration models in 25% of the simulated cycles). Then, to account for failed recruitment to the upper canopy and site index (SI), we calibrated the median-predicted regeneration values for each regeneration model based on existing literature (e.g., Ouzts et al., 2015; Strom & Fulé, 2007) and expert opinion (Randy Fuller, USFS Apache-Sitgreaves National Forest). Regeneration models for numerous sprouting (e.g., Gambel oak, quaking aspen) and mixed-conifer species (e.g., Douglas-fir, white fir) were only initiated during the cycle that followed wildfires or mechanical thin and prescribed fire operations.

After accounting for recruitment, regeneration models were also programmed to respond to SI (FFRI, 2011; USDA, 2020b) and underlying mechanisms of type conversion (e.g., disturbed ponderosa pine forest returning as alligator juniper or Gambel oak woodland) based on revegetation potential (Coop et al., 2020; Winthers et al., 2005) and the presence of sprouting trees species (e.g., Johnson et al., 1962; Kaufmann et al., 2016). For example, if SI was low (<80) the post-wildfire regeneration of ponderosa pine following a high-intensity fire (flame length > 3.66 m) was reduced by 80% for 16 years when revegetation potential was moderate, and by 90% for 40 years when revegetation potential was low, when compared to a moderate-intensity wildfire (Randy Fuller, USFS; Savage & Mast, 2005; Singleton et al., 2021; Appendix S3). Low regeneration potential coincided with high erosion potential occurring on steep slopes. When compared to low-intensity fires, we also enhanced the regeneration of alligator juniper three times the empirical prediction following moderate-intensity wildfire and six times following high-intensity fire to account for its ability to survive wildfire and resprout postfire (Johnson et al., 1962). Regeneration for each species–ecological zone combination was predicted using an FVS written ASCII output file containing current stand information (topographic information and species-specific trees per hectares and basal area). We then wrote regeneration counts to an ASCII file, which was read by FVS and used to update the inventory record for each stand. See Appendix S3 for a full accounting of regeneration parameters for each species within our simulation.

Landscape fuel dynamics

We built a rule set for surface fuel models that deviates from FVS-FFE to account for temporal variability in live fuels, particularly shrubby life-forms (Appendix S4). For every LSim cycle, fuel models are selected based on live woody fuel load, fuel bed depth, plant association (USDA, 1997), canopy attributes, and time since disturbance. Surface fuel models for portions of the study area with tree list data were modeled over time using our custom fuel model logic (Appendix S4) with output data from FVS-FFE. Other FVS-FFE data utilized were canopy cover, canopy base height, canopy bulk density, and stand height. Portions of the study area without tree list data were assigned 2012 surface fuel models and canopy fuels from LANDFIRE data (LANDFIRE, 2013). Other landscape inputs obtained from LANDFIRE for fire modeling included elevation, slope, and aspect (Krasnow et al., 2009; Rollins, 2009). Surface fuel models in our study area use a combination of the original 13 (Anderson, 1982) and new 40 fuel models (Scott & Burgan, 2005).

Wildfire calibration

FSim was calibrated to replicate historical fire occurrence and size in the study area. For this process, we used a 15-year (2000–2014) record of wildfire history from the Kaibab and Coconino NFs (≥4 ha; KCFAST, 2020), and corresponding weather from the Flagstaff Remote Automated Weather Station (RAWS) located 4 km south of Flagstaff, Arizona. With these data, we developed a modified fire distribution file (FDIST file) using a kinked logistic regression with a breaking interaction point at an ERC of 46 (Figure 4, kinked logistic). When compared to a logistic regression derived with FireFamily+ (Figure 4, logistic), the kinked approach allowed us to adequately calibrate the probability of high ERC fire ignitions (Figure 4, empirical data), while retaining the predictive value of numerous RO wildfires with an ERC less than 30. To achieve compatibility with the constraints of the FDIST file (i.e., FDIST file accepts two parameters: intercept and ERC coefficients), we dissected our kinked model and focused on components with ERCs ≥46, which were then extended to low ERC fires resulting in the fire distribution curve in Figure 4. The intercept and ERC coefficients for this curve were then recorded in the wildfire distribution file developed using FireFamily+ (Bradshaw & McCormick, 2000) (Figures 2 and 3). Combined with other calibration parameters in FSim discussed below and in Appendix S5 (e.g., artificially lowing the 80th percentile), this adjustment improved the accuracy of wildfire simulations. For additional information on our kinked logistic model, refer to Appendix S5.

FIGURE 4. Daily large fire probability (FDIST curve, smooth black line) developed using a 15-year (2000–2014) record of wildland fires from the Coconino and Kaibab National Forests and energy release component (ERC) from the Flagstaff, Arizona, Remote Automated Weather Station. The stepped gray line represents fire probabilities classified by ERC. Protection objective alternatives use the portion of the FDIST curve that is to the right of the 80th percentile, where fires to the left of the 80th percentile are suppressed with initial attack success. Resource objective alternatives use the portion of the FDIST curve that is to the right of the 50th percentile, where fires are managed for resource benefit when conditions allow between the 50th and 80th percentile, and fires to the left of the 50th percentile are suppressed with high initial attack success

We developed an ignition probability grid to direct fire ignitions based on the historical occurrence of 9444 wildfires from 1992 to 2015 that ranged in size from 0.004 to 9666 ha (Appendix S5: Figure S1) using data from Short (2017). Much of the northern and southwestern portions of the study area are remote, fuel-limited grass and shrublands with an infrequent fire history. Previous studies conducted in northern Arizona have found climate change to have a limited effect on fire regimes when compared to active land management (O'Donnell et al., 2018). However, future deployments of LSim could investigate this question further using statistically downscaled energy release components to model wildfire using future climate ensembles (Abatzoglou, 2013; Barros et al., 2021) via the FRISK file. Refer to Appendix S5 for more information regarding fire calibration, including seasonal distributions of fire ignitions and area burned.

Forest and wildfire management scenarios

Using a 3 × 2 factorial design (Table 2), we examined potential long-term trade-offs and synergies between restoration treatments (i.e., mechanical thinning followed by prescribed fire) and wildfire management. Specifically, we simulated three restoration treatment alternatives: no treatments (0×), the current pace of treatments (1×), and five times the current pace of treatments (5×). These three restoration treatment alternatives were combined with two wildfire management alternatives: all wildfires managed for protection objectives (PRO), and the option to manage wildfires for resoure objectives (RO) based on ignition location and weather conditions (Table 2, Figure 1b). The business-as-usual scenario is 1×RO. Each of the six scenarios was simulated for 56 years (14 cycles) with 10 replicates to account for stochastic processes within LSim (total = 60 simulations). In total, each scenario simulated 560 wildfire seasons, which adequately reflects the inter-replicate variance of 10,000 seasons (Ager et al., 2021; Appendix S5). The sequence of activities in each cycle, including management, growth and mortality, regeneration, and wildfire effects, is depicted in Figures 2 and 3. In each cycle, 4 years of management and wildfire were simulated.

TABLE 2 Wildfire management alternatives, treatment levels, and scenarios

Wildfire management alternative	Acronym	Treatment level	Scenario
Full suppression to meet protection objectives	PRO	0×	0×PRO
		1×	1×PRO
		5×	5×PRO
Wildfire managed for multiple resource objectives	RO	0×	0×RO
		1×	1×RO (BAUa)
		5×	5×RO

Note: Wildfire management was simulated with two alternatives, full suppression to meet protection objectives (PRO), and the addition of wildfire managed for multiple resource objectives when conditions allow (RO). The pace of treatment was simulated at three levels, which include no restoration treatment (0×), current treatment (1×), and five-time current treatment level (5×) on the Kaibab and Coconino National Forests. The combination of wildfire management alternative and treatment level corresponds to a simulation scenario. a Scenario 1×RO represents business as usual (BAU) on the Kaibab and Coconino National Forests.

The RO fire management alternative was simulated with a new feature in FSim that allows for suppression decisions to be made for each ignition based on location and month of the year. Ignition locations assigned protection objectives (Figure 1b) were either WUI, within 1 km of WUI or contained wet mixed-conifer forests, steep canyons of interior chaparral, or national monuments and other protected sites. The timing for RO management was designed to only allow RO fires during periods of mild weather within the shoulder fire season (first 4 months [January through April] and last 5 months [August–December] of the year). During each month of the shoulder fire season, we adjusted fire management response using FSim calibration parameters. This included modifying perimeter trimming to reduce suppression activity, increasing burn periods to resemble fire growth and burn out operations, and precluding fires that spread beyond 6070 ha to simulate the measured response of quickly suppressing fire during weather conditions likely to lead to extreme fire behavior if managed for an extended period (details in Appendix S5). The upper threshold of 6070 ha was chosen to reflect the largest wildfire used to meet ROs in the study area (e.g., 2009 Wildhorse Fire, 5581 ha).

Planning areas were prioritized for treatment based on potential fire severity and associated basal area losses estimated with the POTFMORT keyword. Active management goals are not always met for each 4-year LSim cycle if there is not enough area to treat based on tree stands exceeding minimum thresholds of merchantable timber, or the minimum passage of time since the last disturbance (Appendix S2: Figure S1). For example, stands were excluded if they burned with low to moderate intensity in the last 10 years including prescribed fire. Stands that burned with high intensity (flame length > 3.66 m) in the last 20 years were also excluded from treatment. Five times the current pace of treatment (5×) was the proposed rate of implementation by the USFS (USDA, 2015a).

Simulation sequence and response metrics

Landscape response to each scenario (Table 2) was assessed by (1) the variability in area burned and fire severity, (2) the combined effects of forest management and wildfire on restoration treatments and the promotion of old-growth forest structure, and (3) landscape heterogeneity as measured by the number and shape of high-severity patches, and their ability to convert forests to woodlands or grasslands. Variability in area burned was assessed with the coefficient of variation (CV). The CV is a relative measure of variability where the wildfire population's standard deviation is divided by the mean. We convert this proportion to percentiles by multiplying the CV by 100. Higher percentile values represent a wildfire population with more variation. The old-growth forest structure was assessed by the scenario's ability to return the landscape's tree loading to its historical range of variation (Reynolds et al., 2013). High-severity patch shape and conversion potential were assessed with a count of disjunct core areas within high-severity patches at least 200 m from the forest edge, and with the average-weighted mean (AM) perimeter-to-area ratio to assess the roundness of the patch perimeter (Helzer & Jelinski, 1999). High-severity patches (defined as those with >75% tree mortality excluding seedlings and saplings with a diameter at breast height <13 cm [Barrett et al., 2010]) were analyzed to understand how restoration and wildfire scenarios might prevent the occurrence of large patches where delayed forest regeneration contributes to type conversion potential (e.g., ponderosa pine to grass or shrub) as reported by Owen et al. (2017, 2020) and Singleton et al. (2021), as well as others. Patch metrics were calculated using functions in FRAGSTATS (McGarigal & Marks, 1995).

RESULTS Area burned and fire severity

The effect of the alternative suppression strategies on the number and size of wildfires was readily apparent over the simulation time frame (Figure 5). The PRO (protection objective) scenarios that managed each fire with a full suppression strategy were dominated by smaller wildfires (Figure 5a), whereas the RO (resource objective) scenarios exhibited larger wildfires on average (Figure 5b), with a greater proportion of the study area burned cumulatively over the course of each replicate (Figure 6a), especially where RO wildfires were allowed (Figure 1b). The cumulative area burned for the PRO scenarios ranged from 13% (1×PRO; 101,000 ha) to 38% (5×PRO; 296,000 ha) of the study area, while the RO scenarios ranged from 44% (1×RO; 342,000 ha) to 76% (5×RO; 592,000 ha) of the study area over the simulation period (Figure 6a). The annual area burned by wildfires in PRO scenarios (Figure 7a) had a greater share of years below the historical mean (2000–2014), while the RO scenarios were more consistent with the mean (Figure 7b). Variability in area burned over time was substantially higher for PRO scenarios compared to RO with the inter-cycle coefficient of variation (CV) in area burned being almost twice as much for the PRO scenarios compared with the RO (CV: 131% vs. 70%). This difference was also observed among the 30 replicates within each scenario (CV: 130% vs. 73%) and among the three treatment levels within each wildfire management scenario (CV: 151% vs. 80%). Thus, the fire management scenarios varied in both the magnitude of area burned and the variability in terms of future predicted fire regimes.

FIGURE 5. Wildfire perimeters for a single replicate over the course of 56 years for the two wildfire management alternatives: (a) protection and (b) resource objective for the 1× treatment level. Wildfire outside the planning areas is not included in our analysis

FIGURE 6. Cumulative area burned by wildfire (a) and cumulative area treated by restoration activities (b) described in Table 1 as a percent of the study area. Bands correspond to 95% confidence intervals. 0×, 1×, and 5× represent treatment levels. PRO, protection objectives; RO, resource objectives. Refer to Table 2 for scenario details

FIGURE 7. Area burned by simulation year showing variability among and within four simulated replicates for the protection (a) and resource (b) objective wildfire management alternatives for the 1× treatment scenario. Also shown is the maximum (27,651 ha) and mean (16,585 ha) annual historical area burned between 2000 and 2014 (KCFAST, 2020). The protection objective scenario (a) has numerous years with limited area burned, while at the same time having numerous years that exceeded the historic maximum. The combined effect is an increase in interannual variability as measured by the coefficient of variation (CV) relative to the resource objective scenario (b) that is dominated by wildfire activity closely resembling the historic mean on many occasions

In terms of fire severity, the RO scenarios had elevated levels of area burned at high severity compared with the PRO scenarios (Figure 8a) with up to 13% of the study area experiencing stand-replacing fire at the end of the simulation period. However, the RO scenarios also had higher levels of low- and moderate-severity fire (Figure 8b), which reversed the trend of high severity as a proportion of the total area burned in each scenario resulting in higher proportions of high-severity fire in PRO scenarios (Figure 8c). The proportion of low- and moderate-severity fire reflected the change in high severity, with higher levels of low and moderate severity for the RO scenarios (Figure 8d). The differences between the fire management alternatives are evident when expressed as annualized outputs (Figure 9) where the annual area burned with high severity for RO scenarios was higher than the PRO scenarios among the three treatment levels (Figure 9b), but the proportion of high-severity fire was consistently lower (Figure 9c). Though RO scenarios had a modest increase in high severity on the landscape compared with PRO scenarios, RO scenarios significantly increased levels of beneficial low- to moderate-severity fire that outweighs the high-severity fire over time.

FIGURE 8. Trends in high- versus low- and moderate-severity wildfire. Cumulative area burned expressed as a percentage of the study area (a, b) versus percent of total fire activity (c, d) for high-severity wildfire (a, c), and low- and moderate-severity wildfire (b, d). Bands correspond to 95% confidence intervals. A larger area burned at high severity for the resource objective alternatives compared with the protection objective alternatives, but proportionately resource objective alternatives burned with lower severity after accounting for total area burned. 0×, 1×, and 5× represent treatment levels. PRO, protection objectives; RO, resource objectives. Refer to Table 2 for scenario details

FIGURE 9. Annual area burned (a) and area burned with high severity as a percent of the study area (b), and area burned with high severity as a percent of total area burned (c). A larger area burned at high severity with the resource objective fire management alternatives compared with the protection objective alternatives (b), but proportionately the resource objective alternatives burned with lower severity after accounting for total area burned (c). 0×, 1×, and 5× represent treatment levels. Refer to Table 2 for scenario details

Combined effects of forest management and fire

The effect of increasing forest management on area burned and fire severity was similar for both wildfire management alternatives (Figure 6a, 8ab and 9a). Area burned decreased between 0× and 1× management, then increased between 1× and 5× management. The increase in area burned between 1× and 5× was amplified in the RO scenario. Increased rates of burning with 5× treatment followed post-treatment increases in trees per hectare due to widespread regeneration across the landscape (Figure 10a) combined with increased rates of fire spread driven by wind penetration into the grass understory (Appendix S6). The 5×RO scenario had the greatest cumulative area burned with high severity across all management types, and 1×PRO had the lowest (Figure 8a).

FIGURE 10. (a) Trees per hectare and (b) basal area in ponderosa pine-dominated stands available for restoration treatments. Solid lines represent ponderosa pine, and dashed lines represent all trees within a stand. Gray-shaded bans represent the interquartile range of the historic range of variation for trees per hectare (56–138 trees) and basal area (9–17 m2/ha) for ponderosa pine forests (Reynolds et al., 2013). Bands correspond to 95% confidence intervals. Note that the 5×RO scenario best maintained all tree species combined and the basal area of ponderosa pine, while also quickly reducing ponderosa pine per acre, which entered the historic range of variation after 25 years. 0×, 1×, and 5× represent treatment levels. All, all tree species; PP, ponderosa pine; PRO, protection objectives; RO, resource objectives. Refer to Table 2 for scenario details

Increasing RO wildfire during 1× treatment did not reduce the area treated (Figure 6b), due to wildfire rarely preempting planned restoration treatments, and increased regeneration following low- to moderate-severity fire that increased the need for future treatments. With 5× treatment levels, RO wildfire reduced cumulative area treated modestly as wildfire began to preempt and maintain treatments more commonly. The average annual area treated for the 1×PRO scenario was 2.5% of the study area, compared with 2.6% for the 1×RO scenario. At higher levels of management, the average annual area treated for 5×PRO was 6.8%, compared with 6.4% for the 5×RO scenario. This equated to RO wildfire increasing treatment needs by 778 ha annually for the 1× treatment level and decreasing treatment needs by 3114 ha annually for the 5× treatment level.

The combined effect of disturbances shows the dynamics between fire and treatments for both fire management alternatives. Over the course of the simulation, area burned at low to moderate severity during RO fire management (Figure 8b, 37%–63% study area) met or exceeded cumulative area treated with 1× restoration treatments (Figure 6b, 38% study area). Conversely, the PRO scenarios had cumulative area burned at low to moderate severity (Figure 8b, 10%–31% study area) that was consistently lower than the area treated with 1× restoration treatments (Figure 6b, 38% study area). Furthermore, area treated with 1× restoration treatments in combination with low- to moderate-severity RO wildfire (75% [= 37% + 38%] to 101% [= 63% + 38%] of the study area combined) met or exceeded area treated with 5× restoration treatments (Figure 6b, 75% of the study area for 5×RO to 83% of the study area for 5×PRO) at the end of the simulation period, which equated the addition of RO management to 4× additional treatment. This key nonadditive effect is further magnified for the 5×RO scenario that quickly restored ponderosa pine density to its historic range of variation in 25 years (Figure 10a) with the added benefit of maintaining a consistent basal area of ponderosa pine and all tree species combined (Figure 10b). In effect, 5×RO best resulted in the retention of old big trees and the promotion of old-growth forest structure by maintaining ponderosa pine basal area and reducing ponderosa pine trees per acre.

Landscape heterogeneity in high-severity fire

Landscape heterogeneity in fire severity as measured by patch metrics showed that the RO scenarios increased the median number and interquartile range of patches that burned at high severity compared with the PRO scenarios (Figure 11a). Apart from fewer high-severity patches when increasing restoration treatments from 0× to 1× concurrent with PRO wildfire management, the addition of restoration treatments and RO wildfire activity increased the prevalence of high-severity patches with 0×RO having the greatest variability and 5×RO having the greatest number. The number of patches between wildfire management alternatives most dramatically departed for 5× management, which is depicted in Figure 12. Consistent with the dramatic increase in high-severity patches for the RO scenarios (Figure 11a, Figure 12b), the number of disjunct core patches at least 200 m from the forest edge was also elevated (Figure 11b) with the increase being dramatic for the 5× treatment level. High-severity patches within RO scenarios with treatment (1× and 5×) were also more complex in shape as indicated with marginally higher AM perimeter-to-area ratios (Figure 11c) (Helzer & Jelinski, 1999). Though disjunct core areas are less likely to contain adequate tree regeneration, RO scenarios with treatment produce high-severity patches that are more complex in shape and more likely to regenerate due to a higher proportion of the patch being adjacent to the forest edge. Complex patches also enhance wildlife habitat and foraging opportunities for goshawk and Mexican spotted owls (Reynolds et al., 1992; USFWS, 2012). Patch size was not discernable among scenarios and is not depicted or discussed.

FIGURE 11. High-severity patch metrics including the number of patches (a), the number of disjunct core areas (b), and area-weighted mean (AM) perimeter-to-area ratio (c). Boxes represent the first and third quartiles across replicates, the internal horizontal line represents the median, and the whiskers represent minimum and maximum values, excluding outliers represented with dots. Resource objective wildfires changed the distribution of high-severity burn patches by increasing their total number and range, their likelihood of containing disjunct core areas, and their edge complexity. 0×, 1×, and 5× represent treatment levels. Refer to Table 2 for scenario details

FIGURE 12. High-severity patches for a single replicate over the course of 56 years for the two wildfire management alternatives: (a) protection and (b) resource objective for the 5× treatment level

DISCUSSION

This study examined the long-term effects of expanding RO wildfire in concert with contemporary restoration treatments. The 778,378-ha study area has a history of fire exclusion and is a national target for accelerated restoration by the US Forest Service. Under current budget levels, restoring the area to a fire-resilient condition will require RO wildfire in addition to mechanical and prescribed fire treatments. This is the first modeling study to examine mixtures of wildfire specifically managed for resource objectives with restoration treatments as practiced on national forests in the region. Overall, adding RO wildfire management resulted in more areas burned at high severity. Though, when taken as a proportion of the total area burned, RO wildfire resulted in a lower proportion of high severity. We also found that managing wildfires for resource objectives diminished the likelihood of irregular fire events over time by making wildfire activity more predictable. Resource objective wildfire management also changed the distribution of high-severity burn patches by increasing their total number and range, their likelihood of containing disjunct core areas, and their edge complexity. These changes diversify wildlife habitat and foraging opportunities for goshawk and Mexican spotted owl and their prey species as old-growth forest structure is enhanced (Reynolds et al., 1992; USFWS, 2012). The RO wildfire management alternative did not significantly reduce the need for restoration treatments by taking their place, whereas 5× treatments had a targeted effect on the landscape that reduced the backlog of restoration needs outlined in FFRI within 12 years (USDA, 2015a).

We did not expect accelerated forest management (5×) to increase area burned relative to business as usual (1×). Other long-term modeling studies reported a steady decline in area burned when contemporary restoration levels were increased (e.g., O'Donnell et al., 2018). Our results suggest restoration treatments followed by prescribed fire trigger asexual regeneration of Gambel oak and Mexican locust (Figure 10a, all trees relative to ponderosa pine) at rates consistent with field data (Abella, 2008; Harrington, 1985; Pavek, 1993). This post-treatment spike in regeneration increased ladder fuels (Figure 10a), while the open canopy simultaneously increased the fine fuel grass component of the understory (Appendix S6: Figure S1) and lowered the wind reduction factor of the canopy, resulting in an increased mid-flame length wind speed (Rebain et al., 2010; Scott & Burgan, 2005). These factors contributed to an increase in the modeled rate of fire spread (Appendix S6: Figure S2) and increased fire intensity predicted by surface fuel models (Scott & Burgan, 2005). As surface fire intensity increased, the minimum wind speed required for torching and crowning decreased (Rebain et al., 2010; Scott & Burgan, 2005). The combined result of these mechanisms was an increase in fire spread and intensity in heavily treated areas that contained widespread ladder fuels, ultimately leading to the increased prevalence of isolated torching and crown fire activity. We believe the model accurately depicted changes in fuels resulting from accelerated forest management (5×). We also observed the occasional occurrence of high-severity patches in both managed and adjacent stands that would have overwise been more amenable to containment or protection, not represented in the model. Regardless, for the current model calibration the 5×RO scenario best restored and maintained manageable lands by steadily reducing ponderosa pine trees per acre (Figure 10a), while stabilizing ponderosa pine and total basal area (Figure 10b). This returned the ponderosa pine landscape to its historical range of variation within 25 years (Figure 10a; Reynolds et al., 2013), which improved wildlife habitat for goshawk and Mexican spotted owl by enhancing old-growth forest structure (Reynolds et al., 1992; USFWS, 2012).

The finding that RO fires resulted in substantially less variability in extreme wildfire events without dramatically increasing the prevalence of large, high-severity patches at risk of type conversion is an important finding in terms of managing wildfires and their impacts to people and ecological conditions. Across each treatment level, the variability in area burned by wildfire in the protection objective (PRO), as measured by CV, was nearly twice that of the RO scenario. This illustrates the wide range of possible wildfire futures for PRO scenarios despite the policy of suppressing all fires as quickly and safely as possible. In fact, under the current treatment level (1×), the PRO policy led to a future with greater variability and occurrences of extreme fire seasons (Figure 7a), which are known to increase high-severity fire resulting in damaged ecosystems (Neary, 2009; Shive et al., 2013), disruptions to communities and human health (Carroll et al., 2011; Liu et al., 2015; Reid et al., 2016), and high suppression costs that increase near WUI (Gorte, 2013; USDA, 2015b), particularly those of high value (Bayham & Yoder, 2020).

Constraints and benefits of using RO wildfire to augment restoration treatments were identified. The increased temporal and spatial variability in restoration treatments due to the random nature of RO wildfire taking the place of some treatments while catalyzing the need for others after the maturation of postfire regeneration (Figures 6b and 10a) constrains RO wildfire's ability to reliably reduce planned treatments on manageable lands. Conversely, the area treated with 1×RO met or exceeded area treated with 5× restoration management, equating RO management to 4× additional treatments, with much of the additional treatments occurring where mechanical treatments are not allowed or feasible.

Our results support the idea that RO wildfire in concert with restoration management on fire-excluded landscapes can potentially enhance fuel and wildfire heterogeneity with positive impacts on long-term forest resilience and old-growth forest structure. When combined with business-as-usual (1×) treatments, the RO alternative only had a slight increase in the count of large disjunct high-severity patches greater than 200 m from the forest edge, and these patches were more complex in shape making them more likely to receive ample seed dispersal from the adjacent forest. Large circular or blocky high-severity patches in ponderosa pine forests have been shown to have seed dispersal limitations (Owen et al., 2017, 2020; Singleton et al., 2021) capable of catalyzing transitions to grass and shrub ecotypes (Tarancón et al., 2014), and facilitating an invasion by non-native plant communities with high fire spread rates (Kerns et al., 2020). Large patches of high-severity fire can have the opposite effect in colder temperate forest systems, where ruderal serotinous conifers (e.g., lodgepole pine) mass regenerate over large areas, thus creating and perpetuating crown fire system in the place of mixed-conifer forests that prior to fire exclusion supported patchy mixed and low-severity fire activity (Churchill et al., 2017). In practice, increases in high-severity fire observed in our RO management scenarios can be avoided beyond our simulations with fine-tuned controls on the day-to-day application of suppression or point zone protection, including techniques designed to lower burn severity such as nighttime burning operations and igniting ridges to allow fires to back-burn downhill.

CONCLUSIONS

This study paves the way for further investigations of the synergies between wildfires and conventional forest restoration to improve the resilience of fire-excluded pine forests. By setting up the landscape for the successful application of RO wildfire when conditions allow (North et al., 2012; Young et al., 2019), land managers can potentially increase the use of natural ignitions to expand and maintain restoration treatments (Huffman et al., 2017, 2018, 2020). However, the expanded use of RO wildfire will require caution to protect spatially diverse ecological and human assets, especially on landscapes fragmented by management restrictions or that have not been previously treated with mechanical thinning and/or prescribed fire. Though RO wildfire had the desired effect of expanding the footprint of low to moderate-severity wildfire, it also increased the risk to protected habitats and unmanaged lands with high fuel loadings. Strategic restoration treatments (i.e., mechanical thinning followed by prescribed fire) in conjunction with RO wildfire policies can enhance forested habitats while protecting fire-sensitive values that are inadvertently exposed to wildfires as they spread from managed to protected areas, which is common on national forests that are fragmented by management restrictions.

Future work with landscape models can help managers envision fire management scenarios to determine how to best minimize extreme events and will become increasingly important as ongoing climate change continues to influence wildfire activity (Abatzoglou et al., 2019; Abatzoglou & Brown, 2012; Abatzoglou & Williams, 2016; Westerling et al., 2006) and tree regeneration (Feddema et al., 2013; Owen et al., 2017, 2020). Future work with the LSim model will leverage climate-FVS to model species envelopes (Crookston, 2014; Crookston et al., 2010; Tarancón et al., 2014), and use statistically downscaled energy release components (Abatzoglou, 2013; Barros et al., 2021) to model wildfire using future climate ensembles. Coupling these climate change model components with the RO fire suppression system reported in this paper will create a novel framework to analyze tipping points for managing RO fires in the future climate.

ACKNOWLEDGMENTS

This research was supported by the US Forest Service, Rocky Mountain Research Station, National Fire Decision Support Center (RMRS-NFDSC), a Joint Venture Agreement between RMRS-NFDSC and Northern Arizona University (17-JV-11221637-123), and the McIntire–Stennis Cooperative Forestry Research Program (G1003488). We thank Luke Sheivlin, GreenWave Professional, LLC, Missoula, MT, for programming support; Robert C. Seli, Rachel Houtman, and Ching-Hsun Huang for assistance during early stages of model development; the land and fire management professionals of northern Arizona; and Michelle Day and two anonymous reviewers for their time and insightful comments.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Simulation and summary data (Young, 2022) used to develop figures are available from the United States Forest Service, Forest Service Research Data Archive: https://doi.org/10.2737/RDS-2022-0022.