Fuel Treatment Longevity in Ponderosa Pine-Dominated Forest 24 Years After Cutting and Prescribed Burning

Fuels reduction treatments to mitigate fire behavior are common in ponderosa pine ecosystems of the western United States. While initial impacts of fuel treatments have been reported, less is known about treatment longevity as live and dead fuels change with time. We analyzed fuel dynamics in ponderosa pine–Douglas-fir forests 21–23 years following experimental fuel reduction designed as two independent studies of cutting combined with burning: one tested a commercial thinning strategy, while a second tested a retention shelterwood strategy to reduce fuels while also restoring ponderosa pine forests. Treated units were harvested in 1992 and half of the units were prescribed burned 1 to 2 years later. After 22 to 23 years post-treatment, few differences in fuel load persist and all treatments have increased ladder fuels in the form of live saplings and seedlings. Canopy fuel loads were lower in treated units compared to untreated control units; however, no other canopy fuel metric differed between treatments. The only persistent difference in surface fuels was in the retention shelterwood, where 1 h fuels were lower in the treated units compared to control units. Crown fire hazard varied greatly, but means were similar between treatments. The increased hazard was driven by increases in live surface fuels from seedlings and saplings in the retention shelterwood, which increased canopy bulk density and reduced canopy base height. The overstory was still dominated by ponderosa pine 22–23 years later for all treatments, but the smaller size classes were primarily Douglas-fir, suggesting that without future disturbance, dominance will shift from pine to Douglas-fir dominated forests. The exception to this was the cut+fall burn treatment in the commercial thinning, where ponderosa pine outnumbered Douglas-fir trees across all size classes. The treatments that included a broadcast prescribed burn killed many existing seedlings and saplings. Our findings support other studies showing fuel reduction and restoration treatments are most successful with a combination of cutting and burning strategies, but also show that fuel treatments in low-elevation dry forests will likely not remain effective for much longer than historical mean fire return intervals.

Fuels reduction treatments to mitigate fire behavior are common in ponderosa pine ecosystems of the western United States. While initial impacts of fuel treatments have been reported, less is known about treatment longevity as live and dead fuels change with time. We analyzed fuel dynamics in ponderosa pine-Douglas-fir forests 21-23 years following experimental fuel reduction designed as two independent studies of cutting combined with burning: one tested a commercial thinning strategy, while a second tested a retention shelterwood strategy to reduce fuels while also restoring ponderosa pine forests. Treated units were harvested in 1992 and half of the units were prescribed burned 1 to 2 years later. After 22 to 23 years post-treatment, few differences in fuel load persist and all treatments have increased ladder fuels in the form of live saplings and seedlings. Canopy fuel loads were lower in treated units compared to untreated control units; however, no other canopy fuel metric differed between treatments. The only persistent difference in surface fuels was in the retention shelterwood, where 1 h fuels were lower in the treated units compared to control units. Crown fire hazard varied greatly, but means were similar between treatments. The increased hazard was driven by increases in live surface fuels from seedlings and saplings in the retention shelterwood, which increased canopy bulk density and reduced canopy base height. The overstory was still dominated by ponderosa pine 22-23 years later for all treatments, but the smaller size classes were primarily Douglas-fir, suggesting that without future disturbance, dominance will shift from pine to Douglas-fir dominated forests. The exception to this was the cut+fall burn treatment in the commercial thinning, where ponderosa pine outnumbered Douglas-fir trees across all size classes. The treatments that included a broadcast prescribed burn killed many existing seedlings and saplings. Our findings support other studies showing fuel reduction and restoration treatments are most successful with a combination of cutting and burning strategies, but also show that fuel treatments in low-elevation dry forests will likely not remain effective for much longer than historical mean fire return intervals.

INTRODUCTION
Fuel reduction projects are designed to reduce wildfire hazard, but the aims of those projects can also include ecological restoration, wildlife habitat enhancement, and forest health improvement (Stephens et al., 2012b;Kalies and Yocom Kent, 2016). In the northern Rocky Mountains of the United States, ponderosa pine (Pinus ponderosa)/Douglas-fir (Pseudotsuga menziesii) forests cover over eight million hectares (Ryker and Losensky, 1983). A century of fire suppression policies has converted many open, seral ponderosa pine forests to dense stands with abundant Douglas-fir regeneration (Smith and Arno, 1999;Hanberry, 2014;North et al., 2015), resulting in reduced vigor, susceptibility to insect infestations, and increased fire hazard (Schoennagel et al., 2004;Hessburg et al., 2015;Hood et al., 2016). For managers interested in moderating fire behavior, the fuel complex remains the factor most amenable to manipulation (Keane, 2015).
Two common approaches to fuels reduction exist-cutting and burning-plus their many permutations. Fuel treatments affect fire behavior by altering both surface and aerial fuel complexes. Surface fuels (both type and quantity), along with weather and topography, influence surface fire intensity (heat output of the fire) and fire severity (ecological impact of the fire). The aerial canopy fuel structure and composition can affect whether a fire transitions from surface fire to crown fire. Fuel treatments affect these two fuel elements both directly-by reducing (or inadvertently increasing) fuel quantities and their contiguity, and indirectly-by altering trajectories of fuel decay, aggradation, and distribution (Jain et al., 2012). Furthermore, because fuel reduction treatments alter the physical environment which in turn alters wind and moisture regimes, future fire behavior may prove more intense (e.g., higher surface fire rates of spread driven by increased wind penetration in open forests) or severe (e.g., higher surface fuel loads from mechanical treatments can increase soil heating) (Agee and Skinner, 2005;Fulé et al., 2012). Seemingly minor differences in fuel treatment prescriptions may cause differences in fuel loads and distribution that become significant with the passage of time.
Fuel treatment longevity, or the persistence of fuel treatment effectiveness, is an important consideration when evaluating fuel treatment alternatives (Reinhardt et al., 2008;Jain et al., 2012;Keane, 2015). While the near-term effects of fuel treatments on fire behavior in the northern Rockies have been well-documented (Smith and Arno, 1999;Fiedler et al., 2010;McIver et al., 2013), the persistence of those effects remains poorly understood (Fulé et al., 2012;Jain et al., 2012). In the short term (<10 years), fuel reduction treatments may reduce surface and canopy fuels, increasing forest resilience to wildfire. Over time, however, vegetation can respond to canopy openings and soil scarification resulting from mechanical treatments, increasing surface and ladder fuels. Greater surface and canopy fuel loads increase forest susceptibility to severe wildfire and decrease fuel treatment longevity (Battaglia et al., 2008;Stephens et al., 2012a;Crotteau et al., 2018). Few studies have analyzed treatment effects for longer than a decade, and there are no studies to date in the Northern Rockies that claim a 20+ year post-treatment response period.
We compared the long-term effects of fuel reduction treatments on surface and canopy fuel loads and hazard in a ponderosa pine-Douglas-fir dominated forest in the northern Rocky Mountains. The experiment was established in 1991 at the Lick Creek Demonstration-Research Forest in western Montana to enable the evaluation of tradeoffs among alternative cutting and burning strategies to reduce fuels and moderate forest fire behavior while restoring historical stand structures and species compositions (Smith and Arno, 1999). The experiment consisted of two independent studies of thinning and retention shelterwood cuttings, with and without prescribed burning treatments.
Our main objectives were to compare (1) surface and canopy fuel loads and (2) crown fire hazard 23 years after treatment initiation. We considered treatments effective if treated units maintained lower surface and canopy fuel levels than control units after 23 years, and if they were resistant to crown fire. We intend this information to inform managers' actions and future planning when considering two particular silvicultural strategies (thinning and retention shelterwood) and prescribed burning for reducing fuel loads and wildfire hazard in dry montane forests of the northern Rocky Mountains.

Study Site
The Lick Creek Demonstration-Research Forest (hereafter Lick Creek) lies on south-facing slopes in the Lick Creek drainage of the Darby Ranger District on the Bitterroot National Forest in southwestern Montana (46 • 5'N, 114 • 15'W) (Figure 1). The elevation varies between 1,300 and 1,500 meters AMSL, with largely 10-30% slopes except at microsites, where slopes range up to 70% (Gruell et al., 1982). Lick Creek has an average winter temperature of −4 degrees C (range −21 to 10 degrees C) and an average summer temperature of 17 degrees C (range 4 to 32 degrees C). The area receives 400 mm of precipitation per year, about half of which occurs as snowfall (RAWS data for Little Rock Creek LRCM8 site near Lick Creek; elevation: 1,678 m, PRISM Climate Group). Soils are typified by Elkner Gravelly Loam, coarse-loamy, mixed, frigid Typic Cryochrepts derived from highly weathered granitic parent material (Gruell et al., 1982;DeLuca and Zouhar, 2000 (Pfister et al., 1977) are Pseudotsuga menziesii/Calamagrostis rubescens, Pinus ponderosa phase, and Pseudotsuga menziesii/Symphoricarpos albus, Calamagrostis rubescens phase. On the lower slopes, the habitat types are Pseudotsuga menziesii/Vaccinium caespitosum, and Pseudotsuga menziesii/Vaccinium globulare, Arctostaphylos uva-ursi phase (Gruell et al., 1982).
Similar to other ponderosa pine/Douglas-fir forests in the northern Rockies (Heyerdahl et al., 2008), the historic, presettlement fire return interval across the Lick Creek drainage averaged 7 years (ranging 3-30 years) (Gruell et al., 1982) and was characterized by low-intensity surface fires (Arno, 1976). In 1906, portions of the Lick Creek area were harvested and additional cuttings have occurred in the 1950s, 1960s, and 1980s (Smith and Arno, 1999.

Treatments
In 1992, ecosystem-based fuels treatments designed to reintroduce low-intensity fire behavior and restore ponderosa pine-dominated forests were implemented, largely based on historical reconstructions of fire frequency in the area (Arno, 1999a). The treatments were also intended to maintain and increase large ponderosa pines in the overstory, improve wildlife habitat for cavity-nesting birds and ungulates, and reduce severe wildfire, insect, and disease hazard. Cutting and burning treatment combinations were tested within two separate, independent studies: a commercial thinning (CT) and a retention shelterwood (SW) (Figure 1). Each study tested four treatments replicated three times for twelve experimental units (1-2 ha each). The treated units were randomly assigned, but the untreated, control units had nonrandom placement due to logistical reasons for the prescribed burns. A permanent plot network was established in 1991 prior to initiation of treatments. Within each unit, a grid of 12 systematically-placed plot centers was installed on 15-40 m spacing depending on the size and shape of the unit, for a total of 144 plots per study. The initial treatment responses were reported in Smith and Arno (1999). Other studies have reported on tree physiology (Sala et al., 2005), reproductive output (Peters and Sala, 2008), and aboveground biomass (Clyatt et al., 2017).
The commercial thinning study was designed to reduce fuels while favoring ponderosa pine in the overstory, boosting tree vigor and growth, and maintaining an even-aged stand structure. The intent of the thinning was to maximize annual increment growth for future yield but not to promote pine regeneration until future cuttings, though regeneration might be achieved. The retention shelterwood study was also designed to reduce fuels while favoring ponderosa pine in the overstory, but with the additional goals of reducing Douglasfir regeneration and initiating development of an uneven-aged stand structure with residual large, old ponderosa pine in the overstory and conditions favorable to pine recruitment in the understory.
We briefly describe aspects of the Lick Creek area that are pertinent to this study, but refer readers to Smith and Arno (1999) for a detailed compilation of the area's rich history of research and management activities. Prior to the treatment implementation, the area designated for thinning had been selectively cut starting in 1907 and partially cut in 1955, 1967, and 1979-1980. In 1991, it supported a 70 year-old second-growth stand of ponderosa pine and Douglas-fir with 370 trees ha −1 , a quadratic mean diameter (QMD) of 27.0 cm, and 21 m 2 ha −1 basal area. Ponderosa pine made up 93% of the trees ha −1 and basal area, and was dominant across all size classes. The area designated for retention shelterwood cutting was an 80-85 yearold second-growth stand of ponderosa pine and Douglas-fir. In 1991, it had 434 trees ha −1 , a QMD of 28.2 cm, and 26.9 m 2 ha −1 basal area. Less than 1% of the overstory was composed of lodgepole pine, with 13% Douglas-fir and 86% ponderosa pine. Saplings were abundant in the understory, with Douglasfir making up 67% (the remainder were ponderosa pine, plus one subalpine fir).
In nine units of each study, harvesting was implemented using chain saw felling. Trees were winched and skidded by crawler tractor along designated trails to road-side landings. Tree tops >15.24 cm in diameter were left on site (Arno, 1999a). Three of the harvested units in each study were not burned (NB), while the remaining six harvested units in each study were subsequently broadcast burned. In the thinning study, two prescribed burn treatments-fall burn (FB) and spring burn (SB)-were conducted post-harvest to examine the influence of burning season on fuel consumption in combination with harvest. Similarly, in the retention shelterwood, two prescribed burn treatments-wet duff burn (WB; 50% moisture, dry weight basis) and dry duff burn (DB; 16% moisture, dry weight basis)-were conducted after harvesting to examine the effect of duff moisture content on fuel consumption and effects in combination with harvest. In the thinning, the FB occurred in fall of 1993 and SB in the spring of 1994; in the shelterwood, the WB and DB prescribed fires both occurred in May of 1993. All prescribed fires were ignited using strip head fires, with distance between strips designed to minimize scorching crowns of the overstory trees. See Harrington (1999a) for additional details on burning conditions. Three unthinned and unburned control units (CO) were included in each study.

Data Collection
Trees (≥15.24 cm DBH) and saplings (>0.10 and <15.24 cm DBH) were measured in 0.04-ha circular plots in 1993-1994 after prescribed burning, in 2005, and in 2015. Each tree's species, diameter, and condition (live/healthy, unhealthy, dead) were recorded. Sapling species and diameter were recorded. In 2015, we also measured height, crown base height, and crown ratio for all trees and a subsample of systematically selected saplings, for a minimum of a 10% sapling sampling, for canopy fuel estimates. Seedlings were tallied by species and size class in 0.004 ha nested circular subplots centered on each plot center in 2005 and 2015. Seedling size classes were 0.5: ≤0.15 m tall, 1: >0.15-0.46 m tall, 2: >0.46-0.76 m tall, 3: >0.76-1.07 m tall, and 4: >1.07-1.37 m tall.
One planar intersect transect (Brown, 1974) per plot was installed in cut-burn treatment units (thinning SB/FB and shelterwood WB/DB) after harvesting and before prescribed burning to quantify woody surface fuels post-harvest (spring 1993 in both studies) and again following burn treatments (spring/summer 1993 in the shelterwood, fall 1993 in the thinning FB and spring 1994 in the thinning SB). These transects were permanently monumented with metal duff spikes at the start (plot center) and end points (15.24 m). Along each transect, 1 h fuels (<0.64 cm diameter) were measured from 0 to 1.8 m, 10 h fuels (≥0.64 and <2.54 cm diameter) were measured from 0 to 1.8 m, 100 h fuels (≥2.54 and <7.62 cm diameter) were measured from 0 to 3.7 m, and sound and rotten 1,000 h (≥7.62 cm diameter) fuels were measured from 0 to 15.24 m. In 2005, transects were remeasured, and transects expanded to all units, with two additional live and dead surface fuels measured at two points (4.60 and 9.10 m) on each transect: (1) litter and duff depth and (2) live/dead herb and shrub height and percent canopy cover. In 2015, all transects were remeasured following 2005 protocols, with the addition of overall average fuel bed depth (m) taken at two points (4.60 and 9.10 m) on each transect.

Data Processing and Analysis
All data were entered into a database using the FEAT/FIREMON Integrated (FFI) software program (Lutes et al., 2009). Archived data is available at https://doi.org/10.2737/RDS-2020-0008 (Lutes et al., 2020). Plot-level surface fuel loads (kg ha −1 ) were output from FFI for statistical analysis, which utilizes formulas from Brown (1974) and Brown et al. (1982). Our analyses compare the long-term effects of fuel reduction treatments on ground and surface fuel load (litter, duff, and downed woody fuels <2 m in height) and canopy fuel load (live fuels regardless of height) and hazard. We categorized surface woody fuels into two groups by particle type: fine woody debris (FWD) consisted of 1, 10, and 100 h fuels, and coarse woody debris (CWD) consisted of sound and rotten 1,000 h fuels. Duff and litter loads are calculated by FFI from the depth measurements and using bulk density values of 0.882559 kg m −3 cm −1 for duff and 0.44126 kg m −3 cm −1 for litter. In the fall of 2014, a natural-ignition wildfire burned a portion of (∼0.5 ha) of one NB treatment unit in the retention shelterwood, eliminating 4 plots from further analysis.
We used 2015 overstory tree characteristics, saplings, live and dead surface fuels, and vegetation plot-level summary output data from FFI as inputs into FuelCalc (Reinhardt et al., 2006) to calculate canopy fuel load, canopy bulk density (CBD), canopy base height (CBH), and crown fire hazard (i.e., Crowning and Torching Indices). FuelCalc defines CBH as the lowest height above ground where CBD is greater than or equal to a calculated threshold CBD value. Threshold CBD is the maximum standlevel CBD × 0.1 for CBD values up to 0.12 kg m −3 . For CBD values >0.12 kg m −3 , threshold CBD is constrained to 0.012 kg m −3 . Canopy bulk density, the mass of canopy fuel load per unit volume (Scott and Reinhardt, 2001), is estimated at the plot-level as the maximum 1.52 m running average in the fuel profile (Lutes, 2014). FuelCalc uses the embedded code of Nexus 1.0 (Scott, 1999) to calculate Crowning and Torching Indices (Scott and Reinhardt, 2001). Nexus creates a custom fire behavior fuel model from the surface and canopy fuel characteristics, but this is not a direct output of FuelCalc. Torching Index is the 6 m wind speed at which crown fire activity can initiate (i.e., surface fire flames begin to burn single and small groups of trees), and Crowning Index is the 6 m wind speed at which active crown fire can occur (i.e., fires can spread from tree crown to tree crown). Torching Index is a function of surface fuel load and composition, surface fuel moisture, canopy base height, slope and wind reduction due to the canopy. Crowning index is a function of canopy bulk density, surface fuel moisture and slope. Maximum wind speed displayed for either index is 160 km h −1 . We used the default "dry" moisture regime.
We obtained average daily wind speeds and gusts using FireFamilyPlus (Bradshaw and McCormick, 2000) to relate the estimated Torching and Crowning Indices to local wind patterns. Daily wind speed and gust data during the primary fire season of July through September were available from 2003 to 2019 from the nearby Remote Automated Weather Station (Little Rock Creek LRCM8).
For the 2015 measurement period, some saplings were missing height values, a required input in FuelCalc to calculate canopy fuels. In the thinning study, about 5% of the values were missing (45 Douglas-fir and 1 ponderosa pine), and in the retention shelterwood study about 20% (1,192 PSME and 6 PIPO) were missing. Therefore, we developed height-diameter equations by study for the two dominant species, ponderosa pine and Douglasfir using simple linear regression. These equations were then used to generate fitted heights for all remaining live sapling records for input into FuelCalc. All live saplings were categorized as "intermediate" crown class and assigned 50% crown ratio. In the thinning, the model for Douglas-fir saplings (Equation 1) explained 91% of the variation while the model for ponderosa pine saplings (Equation 2) explained 85% of the variation. In the shelterwood, the model for Douglas-fir saplings (Equation 3) explained 92% of the variation and the model for ponderosa pine (Equation 4) explained 82% of the variation in height by diameter. Ht = 3.28269 + 1.896637 * DBH /3.28 (1) where "Ht" is calculated height (m) and DBH is measured diameter at breast height (cm). Treatment responses were examined separately by study (thinning, retention shelterwood), as the studies have different site locations and are not designed to be replicates of each other. We evaluated treatment differences in forest structure [tree, sapling, and seedling density, tree basal area, and tree quadratic mean diameter (QMD)] at 1 year (1993/4), 12/13 years (2005), and 22/23 years post-treatment (2015). We also report 2015 diameter distribution of seedlings, saplings, and trees by We used general linear mixed models to compare treatment differences in tree density, QMD, basal area, canopy fuels (load, CBD, CBH, and stand height), surface fuel load by time-lag class (1-1,000 h components, litter, and duff), and crown fire hazard (Crowning and Torching Indices). Treatment units were the experimental unit (n = 3); plots within units were treated as subsamples in the model to account for unitlevel variation. Differences in surface fuel load by component between the 1993/4 and 2015 sampling periods were tested using general linear mixed models, with time as the independent variable and random effects of unit and plot nested within unit. Pairwise differences in treatment means were tested using Tukey's post-hoc test (α = 0.05). Statistical analysis was performed using the GLIMMIX procedure in SAS 9.4 (SAS Institute, Cary, NC, U.S.A.).

Stand Structure
Immediately post-treatment (i.e., time step 1), tree density and basal area were lower in the cut and cut+burn treatments  Figure 3C; F (3,7.961) = 25.18, p = 0.0002]. Sapling density in the NB was higher than the WB (p = 0.0095) but did not differ from the DB (p = 0.1541), and the WB and DB were not different from each other (p = 0.2563; Figure 3C).
In the thinning study, treatment effects persisted through 2005 (time step 12) but by 2015 (time step 22), tree densities in treated units remained largely unchanged, but declined in the CO, such that densities were statistically the same across all four treatments [ Figure 2A; F (3,7.998) = 2.94, p = 0.0988]. Mortality was largely due to mountain pine beetle (Dendroctonus ponderosae) attacks. Ponderosa pine was the dominant species across all size classes (Figure 4, left panel). Basal area increased over the 22-23 years in all treatments (54-62% in treatments and 23% in CO), but was still lower in the cut+burn treatments compared to the CO (Figure 2B; FB: p = 0.0059; SB: p = 0.0046). By time step 22, mean basal area in the NB treatments was similar to the cut+burn treatments (p > 0.9 all comparisons), but because of the variability within the NB units, the NB was also not statistically different than the control (Figure 2B; p = 0.4852). Over the 22-23 years since treatment, QMD increased by 23-26%   p = 0.06276; NB: p = 0.9461; FB: p = 0.0768). Douglas-fir was more prevalent in the sapling strata compared to the tree strata, ranging from a low of 18% in the SB to a high of 56% in the NB (Figure 5, left panel). Seedling density at time step 12 was lowest in the SB compared to the other treatments [ Figure 3B; F (3,132) = 20.7, p < 0.0001]. By time step 22, seedling density was much higher in the SB (p = 0.0005) and FB (p = 0.0001) treatments compared to the CO (1,598 and 1,869 seedlings ha −1 compared to 274 seedlings ha −1 , respectively). Seedling density in the NB was not different than the CO (p = 0.0866) or cut+burn treatments ( Figure 3B; SB: p = 0.3496; FB: 0.187). Similar to the sapling strata, Douglas-fir seedlings were more prominent than in the tree strata, ranging from 28% in the FB to 53% in the NB (Figure 5, left panel).
In the retention shelterwood study, treatment effects on tree density persisted through the 22-23 years post-treatment [ Figure 2D; F (3,8.119) = 6.76, p = 0.0135], with lower densities in the treated units compared to the control. This was despite a 15% decline in control tree density over the time period. Tree densities in the NB and WB increased very slightly, but decreased 10% in the DB from time steps 1 to 22. As in the thinning study, tree density declines in the CO were largely due mountain pine beetle attacks, but mortality in the DB was a result of the prescribed burn (Harrington, 1999b). Ponderosa pine was the dominant species across all size classes and treatments, although the CO had a larger component of Douglas-fir (Figure 4, right panel). Basal area increased over the 22-23 years from 37 to 53% percent in the treatments and 14% in the CO, but was still lower in the NB and cut+burn treatments compared to the CO [ Figure 2E;  Figure 3D; F (3,132) = 1.85, p = 0.1413]. By time step 22, there was a significant treatment effect on seedling density [F (3,128) = 3.07, p < 0.0301], but the only treatment difference was between the CO, which had the highest density with an average of 4,790 seedlings ha −1 and the DB, which had the lowest density at 2,102 seedlings ha −1 (Figure 3D; p = 0.0215). There was a higher percentage of Douglas-fir seedlings compared to ponderosa pine in the CO and NB (92 and 88%) vs. the WB (46%) and DB (69%; Figure 5, right panel).

Canopy Fuels
In the thinning, canopy fuel load was lower in the treated units compared to the CO [F (3,8) = 5.72, p = 0.0217], but there were no differences between the NB and cut+burn treatments ( Figure 6A). In the treated units, canopy fuel load  Smith and Arno (1999) with no standard errors available. Lowercase letters denote statistical treatments differences in duff and litter load within a sampling period at α = 0.05. FWD and CWD loads were not statistically different.
was 28-30% lower than the CO, at about 0.47-0.48 kg m −2 in treated units compared to 0.67 kg m −2 in CO. However, there were no differences in treatment for any other canopy fuel variable (Figures 6B-D). This was likely due to the large variability (Figures 6B-D, compare plot values with treatment means).
In the retention shelterwood, there were also persistent treatment differences in canopy fuel loads [ Figure 6E; F (3,8.036) = 9, p = 0.006]. Fuel loads in the WB and DB were 60 and 56% lower, respectively, than the CO, with average loads of 0.38-0.41 kg m −2 compared to 0.94 kg m −2 . Loads in the NB averaged 0.72 kg m −2 and did not differ between the CO or cut+burn treatments. As in the thinning study, there were no significant differences by treatment for canopy bulk density, canopy base height, or stand height, due to large variability within and between treatment units (Figures 6F-H). Mean CBD was 32% higher in NB units than CO and 56-58% higher than the WB and DB. CBH ranged from 0 to 11 m, with two out of three units within each treatment having median CBH values at or near 0 m, indicating dense ladder fuels.

Surface Fuels
In the thinning, there were no treatment differences for any surface woody fuels component over the course of the study period (Figures 7A,B). It seems plausible that FWD loads were higher in the NB at time steps 0 and 1 compared to the SB and FB, but we do not have these data. Despite 23 years of potential fuel aggradation, surface fuels were generally lower in 2015 ( Table 1). The only appreciable difference in fuel loading between 1993 and 1994 and 2015 occurred in FWD, with reduced 1 and 10 h fuels ( Table 1). By 2015, FWD load was very low in all treatments, ranging from 0.2 kg m −2 in the FB to 0.3 kg m −2 in the NB (Figure 7A). Coarse woody debris did not differ among treatments ( Figure 7B) or from time step 1 to 22 (Table 1). Duff and litter loads were lower in the SB and FB compared to the CO at time step 12, while the NB did not differ from the other treatments ([ Figure 7C In the retention shelterwood, woody surface fuels followed a similar pattern to the thinning (Figure 7). FWD fuels did not differ among treatments in any time steps (Figure 7D). When the individual FWD components were examined, there were no treatment effects in 1, 10, and 100 h loads at time step 12. By time step 22, treatment differences had developed in the 1 h component [F (3,136) = 6.08, p = 0.006], where the NB and WB had lower loads than the CO, while the DB was not different than of the other three treatments. However, 1-hr loads were very low (0.002 to 0.023 kg m −2 ) across all the treatments. FWD loads were generally lower in 2015 than post-treatment in 1993 (Table 1), but statistically, the only change over time was in 10hr fuels, which were 57% lower in the WB (p = 0.0044) and 59% lower in the DB (p = 0.0138). There were no treatment effects in CWD in any time step (Figure 7E), but CWD tended to increase from time steps 1 to 22 (Table 1; WB: p = 0.0697; DB: p = 0.0563). Litter and duff loads were higher in the CO compared to the WB and DB 12 and 22 years post-treatment, with the NB mean load statistically the same at the CO, WB, and DB loads ( Figure 7F). Based on the nearest weather station, the median average daily wind speed during the months of July, August, and September is 10.7 km h −1 , with daily gusts averaging 20.9 km h −1 .

Crown Fire Hazard
In the thinning, although treatment means were similar, Torching Index was extremely variable (Figure 8A). For example, the mean Torching Index in the CO was 82 km h −1 , yet median values for two of the three units were approximately 50 km h −1 , while one of the units was 113 km h −1 and variability was highest in the FB, ranging from median unit-level values of 30 to 160 km h −1 . Crowning Index was generally lower (i.e., lower wind speeds needed to sustain a crown fire) in the CO than in the NB and cut+burn treatments, but again, there was high variability (Figure 8B). Average wind speeds and gusts in the area are typically strong enough to support torching on 19-22% of the plots in all treatments (Figure 8A), but never high enough to support active crown fire ( Figure 8B).
In the retention shelterwood, average Torching Index ranged from a low of 16 km h −1 in the NB to 35-38 km h −1 in the other treatments, but means were not statistically different due to high variability ( Figure 8C). In contrast to the thinning, average wind speeds and gusts are often high enough to support passive crown fires in the retention shelterwood in all treatments, with values below 10 km h −1 for 61-63% of the plots in the CO, WB, and DB and 75% in the NB (Figure 8C). Crown fire hazard covered the entire range of possible values from 0 to 160 km h −1 for units in the CO and both cut+burn treatments. Crowning Index means, at approximately 58 km h −1 in the CO and NB and 93 km h −1 in the WB and DB, were not statistically different and above the threshold needed to cause crown fire in all but a few plots in the NB (Figure 8D).

DISCUSSION
What elements contribute to the persistence or longevity of a fuel treatment's effectiveness? Jain et al. (2012) outlined three determinants of fuel treatment longevity: (1) fuel decay, the degradation of dead surface fuels; (2) fuel growth, the response of residual vegetation (especially advance regeneration) to the treated environment; and (3) fuel recruitment, the establishment of new vegetation (typically tree seedlings) which ultimately becomes live canopy and eventually dead surface fuel. Because forests are dynamic and responsive, treatment longevity depends on how heavily and where cuttings occur, as well as site quality and time (Jain et al., 2012). It is also important to quantify surface, ladder, and canopy fuel loads to help inform management decisions (Stephens et al., 2009). Our research questions focused on determining long-term differences in fuel loads and hazard among cutting and cutting in combination with prescribed burning treatments. The most notable changes in fuels at Lick Creek occurred in the posttreatment responses of the advance regeneration and recruited vegetation resulting from the silvicultural treatments in each study (Figure 3).
By 2015, saplings in the thinning had increased in all treatments except the FB, which was the only treatment with fewer saplings per hectare than the CO and NB (Figure 3). In the shelterwood, saplings steadily increased over time, with the largest increases being in the treated areas compared to the CO, such that no treatment effect remained by time step 22. Saplings consisted almost entirely of Douglas-fir in the shelterwood, whereas most saplings were ponderosa pine in the thinning (Figure 5). Retention shelterwood cuttings are intended to facilitate regeneration while retaining large trees of the desired species. That treatment strategy is often employed as a form of "transformation silviculture" (O'Hara, 2001) to initiate the conversion of even-aged stands toward a multi-aged structure, typically for the promotion of forest restoration objectives (O'Hara, 2014). While the treatments intended to recruit ponderosa pine, residual mature Douglas-fir in the overstory served as a seed source for natural regeneration. Additionally, seedlings and saplings (i.e., advanced regeneration) survived the treatments, especially in cut, unburned units. Arno (1999b) previously reported at Lick Creek that 5 years post-treatment, Douglas-fir advance regeneration averaged over 3,200 stems ha −1 in the NB units compared to near 0 stems ha −1 in cut-andburn (WB and DB) treatments. Even post-treatment Douglas-fir regeneration (seedlings <5 years old) exceeded 760 stems ha −1 in the cut-only units, more than 3-5 times higher than cut-burn units. Prescribed burning following cutting had a longer-term impact on sapling density than cutting alone, likely from killing smaller saplings that were not cut during harvest. In the case of the shelterwood, this prescription of cutting without additional surface fuel treatments, successfully promoted regeneration, but benefited the shade-tolerant species most and undermined the long-term treatment goal of increasing ponderosa pine as the dominant second cohort (Arno, 1999a). The emergent ladder fuels from seedlings and saplings in many units of the retention shelterwood study translated directly to increased canopy fuel loads and hazard from torching and/or crowning.
Harvesting operations inevitably generate some residual slash that affect surface fuel loads and fire behavior (Schwilk et al., 2009;Stephens et al., 2009;Prichard et al., 2010). When wholetree harvesting is utilized, thinning operations have negligible effect on surface fuels. However, when limbs and tree tops are left on site, whether masticated or scattered, the residual slash can increase fireline intensity (Stephens et al., 2009). At Lick Creek, whole trees were harvested but tree tops above 15.24 cm diameter were left on site after felling. Tree limbs from the merchantable bole were removed and pile-burned near roadsides, reducing the potential amount of activity fuels. The residual tree-top slash was either consumed in the prescribed burns or largely decomposed to rotten material by 2015.
Accordingly, there were few long-term responses of surface fuel loads between immediately post-treatment and 22-23 years later. Unfortunately, surface fuels were only measured in the prescribed burn treatments before 2005, so we were not able to examine how fuels changed in the CO and NB over the full time period compared to the cut+burn treatments. In both the thinning and retention shelterwood, FWD surface fuels were lower or unchanged. CWD tended to increase over time but accumulations were variable and statistically non-significant. The general decrease to no change in surface woody fuels in the cut+burn treatments may be because stands opened by harvesting exhibit less self-pruning and therefore drop fewer branches to the forest floor. Additionally, increased sunlight exposure and precipitation through-fall to the forest floor in open stands increases surface fuel decomposition rates (Keane, 2008). Our results are in agreement with Stephens et al. (2012a), who reported fine and coarse woody fuels and duff remained unchanged between 1 and 7 years after cutting and prescribed burning in mixed conifer forests in California. Crotteau et al. (2018) also found FWD and CWD was unchanged in a thinned, burned treatment from 1 to 14 years post-treatment in a ponderosa pine-Douglas-fir forest in Montana.
Thinning and burning in combination have been found to more effectively reduce surface and canopy fuels than either treatment alone (Fulé et al., 2012;Martinson and Omi, 2013;Kalies and Yocom Kent, 2016). Twenty-two years post-treatment, we did not find many differences between cutting vs. cutting followed by prescribed burning. Taken as a whole, the results from the retention shelterwood study weakly suggests that prescribed burning can increase treatment effectiveness longevity by controlling regeneration compared to cutting alone to result in more heterogeneity in crown fuel hazard. In the thinning, crown fire hazard did not differ between treatments or the control 22-23 years following cutting and burning. This is unsurprising, as few differences existed between the treatments in surface and canopy fuels. The large variation within and between treatment units suggests fuels and forest structure are conducive in some areas to support individual to small group torching with wind speeds lower than 50 km h −1 , a common occurrence during the fire season in this area. This variability is important to consider, as wildland fire behavior is influenced by local fuel and stand conditions, rather than average values. The area is likely not at high risk to a large crown fire, as sustained winds of >70 km h −1 would be required.
The shelterwood-which generally had higher surface fuel loads, lower canopy base height, and higher canopy bulk density-is more susceptible to passive crown fire compared to the thinning. Surface fuels were generally higher in the shelterwood but mostly in CWD, not FWD. A larger effect on crown fire hazard was likely from increased overstory tree density and ladder fuels, which increased canopy fuel load in the CO compared to the cut+burn treatments.
The silvicultural prescriptions affected understory and ladder fuels in contrasting ways over time. In the thinning, prescribed burning did not further reduce canopy fuel load or bulk density over cutting alone. In contrast, cutting combined with burning in the shelterwood generally lowered canopy fuel load and bulk density compared to cutting alone. The combination of open conditions, residual Douglas-fir in the overstory, and Douglas-fir advance regeneration in the NB treatment in the shelterwood undermined the longevity of the restoration treatments because the prescription allowed a large number of saplings (2,135 trees ha −1 from 0 to 12 cm DBH) to develop in the 23 years since treatment (Figures 5, 9). This set the stage for an explosion of shadetolerant seedling recruitment and advance regeneration release such that 88% of saplings and 71% of seedlings in the NB are Douglas-fir. These ladder fuels are eroding the treatment effectiveness, and the retention shelterwood is overdue for additional treatments.
Fuel reduction treatments are increasingly being tested by wildfires. In the northern Cascades, several 5-20 yearold fuel reduction treatments were tested by the Tripod Complex of 2006 and in most cases, previous prescribed fires (with and without prior mechanical treatments) reduced wildfire severity (Prichard and Kennedy, 2014). In northern Arizona, the Wallow Fire of 2002 burned through 10 yearold fuel treatments with much lower severity than in untreated stands, as measured by overstory mortality and herbaceous community response (Strom and Fulé, 2007;Waltz et al., 2014). These cases highlight the advantages of fuel reduction treatments for breaking canopy fuel continuity and mitigating fire severity, especially in the wildland-urban interface.
Treatment longevity varies by site productivity but in similar existing studies, fuels have recovered to near pre-treatment levels within 10 years (Martinson and Omi, 2013). However, the lack of available long-term data makes it difficult to determine treatment longevity beyond a decade (Fulé et al., 2012). In the Black Hills, prescribed fire within 10 to 15 years of mechanical treatments is necessary to ensure fuel treatment longevity by reducing ladder fuel development from pine regeneration following mechanical activities (Battaglia et al., 2008). Long-term monitoring plots in ponderosa pine forests of the southern Sierras of California showed 63-84% surface fuel recovery within 10 years of prescribed fire treatments, with fuel loads exceeding pre-treatment levels after 30 years (Keifer et al., 2006).
Our results corroborate the growing emphasis on the need to understand fuel treatment lifespan as well as effectiveness, and to plan for regimes of fuel treatments as vegetation grows and fuels accumulate. In forests that historically had a frequent, low-severity fire regime, fuel treatments and restoration activities align to promote open forest structure, dominated by shade-intolerant species. We assert that fire history studies combined with knowledge of site productivity can be informative for determining fuel treatment lifespans and when re-entry treatments are needed. At Lick Creek, the thinning treatments were conducted on a dry southfacing slope with low understory vegetation loads and modest ponderosa pine and Douglas-fir regeneration. In contrast, the retention shelterwood study occurred near the bottom of the Lick Creek drainage, with a cutting treatment intended to facilitate regeneration. The cutting treatments reduced Douglas-fir in the overstory but did not eliminate the species from the site. As a result, post-treatment and advance regeneration grew dense enough over 23 years for the shelterwood site to be vulnerable to wildfire. The historical 3 to 30 year fire return interval (Gruell et al., 1982) is in agreement with our findings that fuel treatments on the drier sites can remain effective for longer than moister sites.

CONCLUSIONS
In order to select from alternative fuel treatment options, managers need to know the expected longevity of their effectiveness and plan for subsequent entries, which requires an understanding of silvics and historical fire regimes. For our study, treatments were considered effective if they (1) maintained lower canopy and surface fuel loads than control units and (2) reduced crown fire hazard. After 22 to 23 years post-treatment, few differences in fuel load persist and all treatments have increased ladder fuels in the form of live saplings and seedlings. Surface fuels across the site were low for a northern Rockies ponderosa pine/Douglas-fir stand (Ottmar et al., 2007;Keane, 2008). Consequently, the potential for crown fire hazard is similar among treatments.
While the cutting and cutting + burning treatments in this study maintained 30-60% lower canopy fuel loads, the silvicultural strategy (thinning or retention shelterwood) had a significant effect on stand development over two decades. Where a thinning was used, treated stands maintained lower canopy fuel levels. However, when a retention shelterwood was used and shade tolerant species were left on site, regeneration was dominated by Douglas-fir thickets that resulted in much higher densities of saplings and seedlings, shortening the period of fuel treatment effectiveness. As has been demonstrated elsewhere, a fuel treatment cannot be successfully conducted as a single-entry event, but must be part of a regime of treatments employed over time in order to maintain effectiveness. Our findings support other studies showing fuel reduction treatments are most successful with a combination of cutting and burning strategies (Reinhardt et al., 2008;Fulé et al., 2012;Stephens et al., 2012a;Kalies and Yocom Kent, 2016), but also show that fuel treatments in low-elevation dry forests will likely not remain effective for much longer than historical mean fire return intervals.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in the online repository: Forest Service Research Data Archive; https://doi.org/ 10.2737/RDS-2020-0008.

AUTHOR CONTRIBUTIONS
SH and CK designed the research, secured funding, and led analysis. SH, CK, and KB co-wrote the initial manuscript and SH led writing of final version. DL led development of the database and fuel calculations. CS provided advice on data collection and fuel hazard. All authors contributed to further writing and editing.

FUNDING
We acknowledge funding from the Joint Fire Science Program under project JFSP 15-1-07-30. This was a study of the Applied Forest Management Program at the University of Montana, a research and outreach unit of the Montana Forest and Conservation Experiment Station, and the USDA Forest Service, Rocky Mountain Research Station.