Introduction

Restoration is a top management priority for dry forests in the western U.S. due to widespread threats from wildfire, insects, and drought that challenge forest resilience, and because of shifting social expectations that place higher value on the quantity, quality, and reliability of non-timber ecosystem services. Changes in land use over the past century, including logging, grazing, and fire suppression, have led to high tree density in many dry conifer forests (Allen et al. 2002). Many forests of the western US are experiencing uncharacteristic wildfires, leading to negative ecological and social consequences (Westerling et al. 2006; Fornwalt et al. 2016). Forest restoration and fuels reduction treatments that alter forest structure and composition (e.g., thinning or prescribed fire) are often implemented to increase forest resilience to wildfire, insects, and anticipated future climate (Agee and Skinner 2005; Fulé et al. 2012; Negrón-Juárez et al. 2014; Strahan et al. 2016).

Several federal initiatives support forest restoration on public and private lands through focused investments in at-risk landscapes such as the Collaborative Forest Landscape Restoration Program (CFLRP; Schultz et al. 2012), Joint Chiefs’ Landscape Restoration Partnership (JCLRP; Cyphers and Schultz 2019), and the Regional Conservation Partnership Program (RCPP; Villar and Seidl 2014). These programs all have a common goal to improve outcomes for ecological goals that transcend ownership boundaries, such as water and air quality, wildlife habitat, and ecosystem resilience to disturbance. They also aim to increase the pace and scale of forest restoration and fuels reduction to achieve meaningful outcomes at the landscape-scale. This raises the need for tools to evaluate and prioritize opportunities based on current conditions and to monitor and report on program accomplishments (Bestelmeyer et al. 2010; Schultz et al. 2012).

Notable advancements in landscape assessment and prioritization have been made in the field of wildfire risk assessment, which has adopted flexible, quantitative frameworks for estimating potential wildfire risks to critical resources and infrastructure, terrestrial and aquatic habitat, and timber (Ager et al. 2013; Thompson et al. 2013a). Other studies have extended the risk assessment framework to estimate the avoided impacts from fuel treatment in terms of post-fire erosion and sedimentation, smoke production, and avian richness, among others (Stevens et al. 2016; Jones et al. 2017; Gannon et al. 2019). These studies often emphasize effects on landscape-scale processes (e.g., fire, erosion, wildlife habitat) that depend on treatment extent, placement, spatial configuration, and type (Ager et al. 2013; Stevens et al. 2016; Jones et al. 2017; Gannon et al. 2019). A goal of many conservation programs is to quantify the impact of investments in conservation efforts on a broad suite of natural resources and consider the effects of spatial scale and landscape context (Bestelmeyer et al. 2010; Briske 2011; Cannon et al. 2019; Metz and Rewa 2019). Landscape-scale simulation studies can guide prioritization, design and placement of conservation treatments within large planning areas, and provide information on potential ecological effects and tradeoffs (Bestelmeyer et al. 2010; Briske 2011; Thompson et al. 2013a; Jones et al. 2017; Gannon et al. 2019; Cannon et al. 2019).

Along with objectives to reduce the consequences of wildfire, many restoration treatments in dry conifer forests focus on diversifying forest structure, composition, and age classes to increase landscape resilience to disturbance. High levels of structural heterogeneity, a hallmark of many fire-frequent pine forests including those dominated by Pinus ponderosa (Lawson & C. Lawson), P. palustris (Mill.), and P. resinosa (Aiton), drives many aspects of forest dynamics, diversity, and resource availability (Platt and Rathburn 1993; Peck et al. 2012; Matonis and Binkley 2016). Guidelines for restoration of ponderosa pine-dominated forests, for example, emphasize restoration of variability in forest structure including increasing horizontal spatial complexity at both fine and coarse scales (Reynolds et al. 2013; Addington et al. 2018). Maintaining high structural variability with mechanical treatments can be relatively challenging, and outcomes can range from heterogeneous to homogeneous depending on implementation methods (Maher et al. 2019), initial conditions (Ziegler et al. 2017), and logistical constraints (Cannon et al. 2018). However, studies of treatment effects on landscape-scale wildfire risk and other ecosystem characteristics typically simulate fuel treatments with post-treatment target conditions (e.g., Thompson et al. 2013b; Stevens et al. 2016), or proportional changes in fuel structure (Jones et al. 2017; Gannon et al. 2019) to mimic completed, planned, or hypothetical fuel mitigation treatments without explicitly capturing spatial goals. Incorporating spatially complex treatment patterns is important as it impacts fire intensity and effects (Cochrane et al. 2012), diversity of herbaceous communities (Matonis and Binkley 2016), resource availability and successional trajectories (Sanchez Meador et al. 2009), and wildlife habitat quality (Reynolds et al. 2013).

In addition to aims to reduce fire hazard and enhance landscape heterogeneity, conservation practices in dry conifer forests can emphasize decreasing the potential for post-fire erosion (Jones et al. 2017; Gannon et al. 2019). Post-fire erosion is of particular concern in watersheds that serve large population centers and have soils with high erosion susceptibility (Moody and Martin 2001; Graham 2003).

Design of forest conservation practices in ponderosa pine forests can vary dramatically, depending on whether the focus is reducing crown fire potential or restoring forest structure to within the historical range of variability. Hazardous fuels reduction treatments favor removal of small diameter ladder fuels and retention of large trees of fire resistant species to decrease crown fire potential and increase the odds of residual tree survivorship in future fires (Agee and Skinner 2005; Fulé et al. 2012; Ziegler et al. 2017).

Landscape-scale prioritization and effects assessments may be improved by addressing these additional goals of restoration. However, uncertainties remain regarding the extent to which purported benefits of various conservation approaches may differ when spatial complexities of restoration treatments and benefits are explicitly incorporated. Here we demonstrate a framework for dry forests to (a) simulate forest restoration, (b) assess how ecological effects of restoration vary spatially, and (c) compare the effects of two conservation practice scenarios. As an example, we present an analysis of ecological effects on multiple natural resources in the Upper South Platte Watershed of Colorado (Fig. 1) in response to two variations on forest thinning and fuels reduction: one that simulates proportional reduction in fuels, and a second that incorporates objectives to restore historical variation in forest structure. We also map ecological outcomes and a composite outcome metric to examine the effects of multiple conservation practice scenarios for prioritization.

Fig. 1
figure 1

a Vegetation zones of the Upper South Platte study area as defined by Kauffman (Kaufmann et al. 2006), and National Hydrography Dataset Plus catchment delineations. b Delineation of 300 forested catchments within the upper and lower montane zones. Forest type is reclassified from existing vegetation type (LANDFIRE; Rollins 2009). Coordinates are in UTM Zone 13N (m)

Full size image

Methods

Study area

To inform landscape-scale conservation decisions, we used catchments (small watersheds) from the medium resolution National Hydrography Dataset Plus (NHDPlus v2; USEPA and USGS 2012) as the primary analysis units. A total of 1010 catchments fall within the Upper South Platte Watershed (USP). We included catchments in the analysis if they were (1) ≥ 75% forest cover (as mapped by LANDFIRE, Rollins 2009), (2) within the upper and lower montane zones (defined below), (3) unimpacted by the 2002 Hayman fire where forest thinning is unnecessary due to widespread effects of high severity fire, and (4) ≥ 40.5 ha (100 acres) in size. The final set of 300 catchments range in size from 41 to 5000 ha (mean = 570 ha; Fig. 1a). Montane forests in the study area are primarily ponderosa pine and Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) with minor components of aspen (Populus tremuloides Michx.) and lodgepole pine (Pinus contorta Douglas ex Loudon) (Fig. 1b) (Rollins 2009). Because management of lodgepole pine forests differs considerably from dry conifer forests, we did not simulate treatments in lodgepole pine forests, although they were included when assessing landscape-scale metrics for heterogeneity, fire behavior, and erosion modeling.

Conservation practice scenarios

The framework for the analysis is to first simulate the primary effects of management on forest and fuel structure and to then estimate treatment effects on selected ecological metrics using process-specific models (Fig. 2). Management activity and research in ponderosa pine forests often centers around reducing hazardous fuels (Fulé et al. 2012), yet more comprehensive restoration treatments which also aim to simultaneously restore structural complexity are increasingly included in desired restoration conditions (Underhill et al. 2014). Thus, we compared three scenarios including an untreated landscape (UNT) and two forest management options: fuel hazard reduction (FHR) and spatially complex restoration (RES) (Fig. 2a).

Fig. 2
figure 2

Schematic diagram demonstrating framework for comparing conservation outcomes of two restoration approaches. Landscape-scale data on forest fuel structure are used as inputs to simulate restoration treatments by (a) adjusting fuel traits based on empirical data or desired conditions. Resulting simulated fuelscapes are then (b) used as inputs to an array of ecological models to (c) estimate conservation impacts on landscape heterogeneity, fire behavior, and soil erosion. Figure adapted from Cannon et al. (2019)

Full size image

Each scenario is represented by a set of raster forest structure and fuel attributes that become the inputs to an array of landscape heterogeneity, fire behavior, and erosion models to estimate the associated ecological effects (Fig. 2b and c). We used 30-m resolution forest structure, fuel, and topography data representing 2014 landscape conditions from LANDFIRE (2016) to characterize baseline conditions for assessing the effects of fuel treatments (Fig. 2a). Forest structure and fuels are described by canopy cover, canopy bulk density, canopy base height, canopy height, and categorical fire behavior fuel model (Scott and Burgan 2005). We adjusted the baseline data for lodgepole pine by lowering the canopy base height by 20% and changing the fire behavior fuel model from moderate load conifer litter to high load conifer litter (Scott and Burgan 2005) to better reflect recent fire behavior observations (Moriarty et al. 2019). Starting with the baseline fuelscape, we simulated each conservation practice design by adjusting fuel structure data according to the expected changes in fuel structure for each treatment. This approach is commonly used to explore landscape-scale impacts of management practices on fire behavior and risk (Scott et al. 2013; Thompson et al. 2013a; Stevens et al. 2016; Jones et al. 2017; Gannon et al. 2019). To better understand the hypothetical benefits of complete restoration at the catchment-level, we make the simplifying assumption that there are no feasibility constraints on treatment. However, real-world application of this model for prioritization will necessarily consider what areas are suitable for treatment based on landownership, accessibility, operability, budget, and regulatory constraints.

Simulating fuel hazard reduction (FHR) treatments

The primary goals of FHR treatments are to reduce fire hazard by lowering tree density and surface fuels coupled with increased canopy base height (Agee and Skinner 2005). To simulate landscape-scale FHR treatments, we multiplied baseline LANDFIRE canopy fuel variables by adjustment factors based on typical outcomes of fuel hazard reduction. We applied an adjustment factor of 0.7 for canopy cover across the study area (e.g., 30% reduction, Fulé et al. 2012). Mechanical thinning typically removes small trees, thus we applied an adjustment factor of 1.2 for canopy height and canopy base height (20% increases, Ziegler et al. 2017). We adjusted canopy bulk density using an adjustment factor based on the reduction of canopy cover [(CCfinal/CCinitial) × 0.8 + 0.18, based on figures in Keane et al. 2005)]. To simulate changes in surface fuel loads, we adjusted the fire behavior fuel model (FBFM) to account for the addition of activity fuels from thinning (Stephens et al. 2009; Fulé et al. 2012). Briefly, FBFMs are used to predict fire behavior across a range of weather conditions for various fuel complexes (Scott and Burgan 2005). For timber litter and slash–blowdown FBFMs, we assumed fuel models would transition to a similar model with a higher flame length (Table S1). If thinning in a timber litter model led to canopy cover reduction below 30% cover, a significant understory response may be expected (Jameson 1967); in these cases, we transitioned fuels to a low load timber-grass-shrub model. We assumed none of the treatments would alter fire behavior in grass–shrub models, except for areas with a very high load timber-shrub model, which is expected to be reduced with thinning treatments. In this case, FBFMs were transitioned to a low load timber-grass-shrub model. For other grass–shrub, grass, and non-burnable fuel types, we assumed that none of the treatments would alter FBFMs. Transitions between all FBFMs can be found in Table S1.

Simulating restoration (RES) treatments

Restoration objectives in ponderosa pine forests aim to increase resilience to disturbances such as fire, insects, and drought by reducing tree density and maintaining or increasing variability in forest structure (Addington et al. 2018). Debates continue about the use of historical range of variability (HRV) to define restoration targets, especially where future climate conditions are predicted to differ from historical ones. In many circumstances, HRV is still a reasonable and objective foundation for developing desired restoration conditions (Keane et al. 2009). Furthermore, restoration of historical conditions in fire suppressed systems are generally consistent with objectives to reduce susceptibility to catastrophic wildfire (Fulé 2008). Local forest reconstructions suggest that current forest densities are higher than they were historically, but historical densities were varied and generally increased with elevation (Battaglia et al. 2018). Regional forest restoration planning guides also emphasize maintaining a diversity of forest densities across topographic moisture gradients (Addington et al. 2018). Thus, we simulate RES treatments with an algorithm that reduces tree density toward historical levels, while ensuring a historically appropriate range of variability between and within topographic and elevation gradients. We (1) classified catchments within the study landscape based on physiographic parameters (elevation and moisture gradients, discussed below), (2) developed desired forest structure distributions by physiographic setting based on published forest reconstructions (ca. 1860 from Brown et al. 2015 and Battaglia et al. 2018) and regional restoration guides (Addington et al. 2018), and (3) simulated restoration treatments such that forest structure and variability in each catchment approached desired conditions for the physiographic setting.

We classified each catchment within the study area as either upper montane or lower montane based on elevation criteria (Kaufmann et al. 2007) (Fig. 1 left). We then developed a simple soil moisture index (SMI) that classifies pixels within catchments into wetter or drier regions by combining standard topographic metrics (contributing area and slope) and expected solar insolation (Western et al. 1999), which can predict patterns in soil moisture attributable to topography. We calculated topographic wetness index (Beven and Freer 2001) using the dynatopmodel package in R. We calculated a raster model of solar radiation (kW m−2) during the growing season (May–September) using ArcGIS 10.3 (Fu and Rich 2002) using a solar transmissivity of 0.625 (mean of Denver and Colorado Springs during the growing season; Knapp et al. 1980). We define an index of soil moisture as

$$SMI=R\left(\mathrm{log}\left(TWI\right)\right)*\left(1-R\left(SR\right)\right),$$
(1)

where SMI is soil moisture index, TWI is topographic wetness index, SR is solar radiation, and R is a function to linearly rescale TWI and SR between the ranges of 0 to 1. We classified areas above the study area median (SWI = 0.2) as wet and areas below the median as dry (Fig. 3a).

Fig. 3
figure 3

Demonstration of restoration simulation (RestSim) algorithm used to simulate landscape-scale thinning to achieve desired levels of canopy cover mean and variability. a Classification of wetter and drier portions of catchments using soil moisture index (SMI, Eq. 1). b Assignment of desired canopy cover mean and variability to wet and dry portions of catchments based on SMI. c Canopy reduction to meet target objectives and incrementally increase target canopy cover as needed to meet silvicultural constraints

Full size image

To develop landscape-scale prescriptions for the RES treatment, we collated published data on historical density and estimated canopy cover from available literature (Tables 1 and 2 in Brown et al. 2015) and developed a predictive equation linking estimated canopy cover and stand density (r2 = 0.911).

Table 1 Parameterization of desired canopy cover conditions for combinations of elevation and moisture gradients using a gamma distribution to mimic landscape scale desired conditions for canopy cover mean (± s.d.) based on Brown et al. (2015), Battaglia et al. (2018), and Addington et al. (2018)
Full size table
$$\mathrm{\%}CC=0.1438*d+6.9035,$$
(2)

where %CC is estimated canopy cover and d is tree density in trees ha−1. We set mean desired catchment scale canopy cover to 20% for the lower montane region and 30% for the upper montane region corresponding to historical (ca. 1860) densities of 97 and 163 trees ha−1 for the lower and upper montane region, respectively, reported in Battaglia et al. (2018). Restoration guidelines often recommend lower densities, but relatively higher overstory retention in wetter aspects, and lower retention in drier aspects, as is congruent with soil moisture availability (Addington et al. 2018; Cannon et al. 2018). Canopy cover variability is expected to be greater on wetter relative to drier slopes due to interactions between moisture, topography, and fire behavior (Addington et al. 2018). To capture these desired conditions, we specified desired canopy cover and variability by catchment using a gamma distribution parameterized by physiographic setting, and we selected parameters α and λ to obtain a distribution that results in patterns of canopy cover and variability to parallel desired conditions (Table 1 and Fig. S1). Desired conditions vary dramatically by region, forest type, and management objectives and constraints can be formulated and modeled using this framework.

Restoration simulation (RestSim) algorithm

To simulate restoration of desired forest structure in each catchment, we developed Restoration Simulation R code functions (RestSim algorithm available online at https://github.com/jbcannon/landscape_restoration_simulator). The RestSim algorithm aims to reduce forest canopy cover to the specified distribution by catchment and physiographic setting. The algorithm first assigns each catchment to an elevation zone using a majority pixel assignment (lower vs. upper montane; Fig. 1). Pixels within the catchment are then assigned into wet and dry categories based on soil moisture index (wetter or drier) (see Eq. 1; Fig. 3a). For each set of wet or dry pixels, RestSim assigns a target canopy cover distribution using gamma distribution parameters corresponding to the elevation zone and moisture class for wet and dry portions of each catchment individually (Table 1; Fig. 3b).

To adjust LANDFIRE canopy cover pixels to target levels, a target canopy cover distribution is generated using the appropriate gamma distribution parameters, and target canopy cover values from the gamma distribution are assigned to pixels in each wet or dry set of pixels by ranking the continuous SMI values present in each set, and assigning the lowest canopy cover pixels from the distribution to the lowest (driest) SMI values. Likewise, the highest canopy cover pixels are assigned to the highest (wetter) SWI pixels (Fig. 3b). We imposed several basic silvicultural constraints: mechanical restoration treatments can only remove trees (not add them), and canopy cover cannot be negative. In some cases, retention of canopy cover in wet areas was not high enough to offset reductions in drier areas, and prescription targets could not be met. In these cases, we iteratively increased the target mean canopy cover (i.e., increasing α by 2% iteratively) until targets were achievable (Fig. 3c). In some cases, existing mean canopy cover was lower than target mean canopy cover before simulation, thus some catchments remained untreated. Finally, we merged canopy cover layers for wet and dry regions for all catchments to generate simulated fuelscapes (Fig. 3c). A comparison of simulated FHR and RES treatments for one catchment is shown in Fig. 4. To simulate changes in forest and fuel metrics besides canopy cover, we used the same methods as for FHR treatments.

Fig. 4
figure 4

Comparison of example catchment in the Upper South Platte watershed showing (a) existing canopy cover, (b) canopy cover following homogeneous fuel hazard reduction (FHR) treatments, and (c) canopy cover following restoration (RES) simulation with the goal of restoring historic canopy cover while maintaining variability.

Full size image

Treatment effects on landscape-scale forest structure

Restoration of variable landscape structure is generally thought to enhance biodiversity. For example, Matonis and Binkley (2016) found higher herbaceous cover and richness in the interior of large gaps following restoration treatments in ponderosa pine ecosystems. Landscape complexity can also benefit avian species. For example, forest treatments that benefit Northern Goshawk (Accipiter gentilis) emphasize interspersion of closed canopy forests with open herbaceous areas, and the retention of older tree groups, to benefit the various organisms in the goshawk’s food web (Reynolds et al. 2013). Such benefits may apply more broadly to avian communities as Latif et al. (2020) found that avian richness increased with measurements of forest complexity like the perimeter-to-area ratio of open forests. To evaluate the effects of each conservation approach on landscape-scale forest complexity, we quantified changes in canopy cover, large gap cover, landscape complexity, and contagion index from the canopy cover layers resulting from each treatment simulation using the raster package in R (Hijmans and van Etten 2016). We quantified large gap cover by identifying all areas with < 10% canopy cover on patches with a radius of 30 m (0.3 ha) or more. To estimate landscape complexity, we used the Shannon diversity index (Shannon and Weaver 1963). This index is typically used to assess diversity in species assemblages but can be applied to other aspects of ecological communities to assess complexity (Turner et al. 2001). We calculated landscape complexity using proportions of each catchment represented by ten classes of canopy cover (e.g., 0–10%, 10–20%, etc.) using the following equation.

$$LC= {\sum }_{i=1}^{n}\frac{{p}_{i}\mathrm{ln}{p}_{i}}{\mathrm{ln}n}$$
(3)

where LC is landscape complexity, pi is the proportion of a catchment with canopy cover within the range of class i, and n represents the number of canopy cover classes (n = 10). This metric results in an estimate of landscape complexity ranging from 0 to 1, with 1 representing equal representation across each canopy cover class. In addition, based on desired conditions related to generating a landscape mosaic with variability of patch types, we assessed changes in spatial configuration of landscape patches using a contagion index (RC2 in Li and Reynolds 1993)

$$RC=1+ \sum_{i=1}^{n}\sum_{j=1}^{n}{P}_{ij}\mathrm{ln}{P}_{ij}/2\mathrm{ln}n$$
(4)

where RC is relative contagion index, Pij is the probability that two randomly chosen adjacent pixels belong to canopy class i and j, respectively, and n is the number of canopy cover classes (n = 10). The index ranges from 0 to 1 and quantifies adjacency of different combinations of patch types, and is sensitive to the number of patch types and the spatial configuration (Li and Reynolds 1993). Larger values of RC represent larger continuous patches of the same canopy cover classes. We calculated RC using custom functions (available at https://github.com/jbcannon/spatial_heterogeneity_metrics).

Treatment effects on fire hazard

To quantify the influence of management on fire metrics relevant in assessing ecological effects and future potential for prescribed fire, we simulated fire behavior for each management scenario under both prescribed fire and extreme wildfire weather conditions using FlamMap 6.0 (Finney et al. 2019). We simulated prescribed fire weather conditions using a scenario with moderate fuel moisture and fuels two-thirds cured (D3L2 scenario, Scott and Burgan 2005) and a 4.5 m s−1 wind speed based on input from prescribed fire managers in the region. For simulating extreme conditions, we processed weather data from 15 Remote Automated Weather Stations (RAWS) from within and nearby (< 20 km) the Upper South Platte Watershed using FireFamilyPlus 4.1 (Bradshaw and McCormick 2000) to characterize 97th percentile (extreme) fire weather conditions for the April 01 to October 31 fire season (Table 2). Given the uncertainty in wind direction, we used the wind blowing uphill option in FlamMap with a default wind direction of 225° to model a worst-case scenario. The 10 min average 1 pm wind speeds from RAWS were converted to 1 min average wind speeds to better reflect potential fire behavior (Crosby and Chandler 1966), and we generated estimates of fireline intensity and fire type (surface, passive crown fire, or active crown fire; Scott and Reinhardt 2001).

Table 2 Weather conditions used to run prescribed fire and extreme fire weather simulations
Full size table

Based on model outputs, we derived three ecologically relevant fire behavior metrics including an integrated hazard index, an index of large high severity patches, and the proportion of the landscape expected to burn under surface fire conditions as a proxy for the ease with which prescribed fire may be incorporated into future management. We multiplied spatially explicit estimates of fireline intensity under extreme conditions (Table 2) (in kW m−1) from FlamMap by modeled burn probability (Short et al. 2016) to produce an integrated measure of fire hazard (Scott et al. 2012). To quantify changes in the probability of large high severity burn patches, we calculated a largest high severity patch index as the size of the largest contiguous patch (in ha) modeled to burn as active crown fire under extreme conditions (Table 2) for each catchment and multiplied this by the average burn probability for the catchment. Lastly, as a proxy for the potential use of prescribed fire or managed wildfire, we calculated the percent cover of each catchment predicted to burn as surface fire under prescribed fire conditions (Table 2).

Treatment effects on erosion hazard

Post-fire erosion was modeled using a spatial implementation of the Revised Universal Soil Loss Equation (Renard et al. 1997) following methods described in Gannon et al. (2019). RUSLE predicts gross erosion (Mg ha−1 yr−1) as the product of factors for rainfall erosivity (R), soil erodibility (K), length and slope (LS), cover (C), and support practices (P). Rainfall erosivity is the product of storm kinetic energy per unit area and maximum rainfall intensity (Renard et al. 1997). Here, we use the median annual rainfall erosivity of 615 MJ mm ha-1 h-1 from weather stations representative of the regional climate (Gannon et al. 2019). We summarized baseline K for the top 15 cm of soil from the Soil Survey Geographic Database (SSURGO) where available and the State Soil Geographic Database (STATSGO) as needed to fill gaps (NRCS Soil Survey Staff 2016). LS was calculated with terrain analysis of a 30-m resolution digital elevation model (USEPA and USGS 2012) following the methods of Winchell et al. (2008) with modifications to restrict hillslope length to 300 m and LS to 72.15 to not exceed the range of values suggested by Renard et al. (1997). We assigned baseline C by mapping previously reported values from the literature (Gannon et al. 2019) to existing vegetation type from LANDFIRE (2016).

To estimate post-fire increase in erosion, C and K factors were modified based on empirical observations from fires in the region (Larsen and MacDonald 2007) to reflect the large increase in erosion with burn severity (Benavides-Solorio and MacDonald 2001). Then, baseline and post-fire (pf) erosion were differenced per Eq. 5 to predict the 1st-year post-fire increase in erosion.

$$A =R\times LS\times \left[({K}_{pf}\times {C}_{pf}) - (K\times C)\right]$$
(5)

Hillslope erosion is typically elevated for 2–5 years after fire in the Colorado Front Range (Benavides-Solorio and MacDonald 2001; Wagenbrenner et al. 2006; Robichaud et al. 2013a, b), so we multiplied the 1st year sediment yield by a factor of 2.1 to estimate the total cumulative post-fire yield (Gannon et al. 2019). RUSLE can predict unrealistically high erosion on very steep slopes, so we limit erosion predictions to the maximum of 200 Mg ha−1 yr−1 observed at hillslope scale in the western U.S. (Moody and Martin 2009). RUSLE predicts gross hillslope erosion, much of which is captured in uplands and not delivered to streams. We estimated the proportion of hillslope erosion contributed to streams using an empirical model of post-wildfire sediment delivery ratio (SDR) from the western U.S. (annual length ratio model in Wagenbrenner and Robichaud 2014). Streams were defined as all pixels with a contributing area greater than 10.8 ha (Henkle et al. 2011). Post-fire SDR was estimated with Eq. 6 based on the flow path length from each pixel to the nearest stream channel as the “catchment length” and the flow path length across the pixel as the “plot length.” We calculated flow path length from a 30-m resolution DEM (USEPA and USGS 2012) in ArcGIS 10.3. The resulting catchment net sediment yield predictions fell far short of post-fire observations from the study watershed (Robichaud et al. 2008, 2013b), so we doubled the predicted SDR as a rough calibration to the limited field observations, and we assigned channel pixels a SDR of 1.

$$log\left(SDR\right) = -0.56 - 0.0094\times \left(\frac{\text{Flow path length to nearest channel}}{\text{Flow path length across pixel}}\right)$$
(6)

The total mass of fire-related sediment (Mg) delivered from a catchment to the stream network (\(TS\)) is the sum product of the hillslope erosion rate (\(A\)), multi-year empirical correction factor, and SDR for burned pixels in the catchment (Eq. 7).

$$TS={\sum }_{i=1}^{N}{A}_{i}\times 2.1\times {SDR}_{i} \times 0.09\frac{ha}{pixel}$$
(7)

Statistical methods

We used linear mixed models to examine the effects of each conservation approach on variables related to forest structural heterogeneity, fire hazard, and erosion risk. For each model, we included the fixed effect of the treatments (untreated, fuel hazard reduction, and restoration) as well as a random effect for each catchment using the lmerTest package in R. To test for differences between individual treatment means, we used post-hoc Tukey’s tests with a Bonferroni–Holm correction to account for multiple comparisons using the glht package (Holm 1979). To test whether tradeoffs or complementary responses among ecological variables differed among conservation practices, we used a linear mixed model to examine whether the relationship between the change in any two ecological variables differed among conservation practices using the following formula (Δvar2–Δvar1 × treatment) including a random effect for catchment. For all analyses, an individual catchment was the unit of observation.

Results

Effects of conservation practice designs on forest structural heterogeneity, fire hazard, and erosion risk

We found that both simulated fuel hazard reduction (FHR) and restoration (RES) treatments significantly affected several ecological variables related to forest structure, fire behavior, and erosion risk. RES treatments were more effective at moving forest structural heterogeneity towards desired conditions; whereas FHR treatments and RES were similar in their effects on variables related to fire behavior and erosion risk (Fig. 5). Treatment simulations reduced mean canopy cover from 41.9% in untreated catchments to 30.3% (padj < 0.001) and 23.9% (padj < 0.001) in FHR and RES treatments respectively (Fig. 5a). FHR treatments had no discernable impact on gap cover relative to untreated catchments (p = 0.999), whereas RES treatments tripled gap cover from 2.3% in untreated catchments to 7.0% (p < 0.001) (Fig. 5b). Whereas FHR treatments decreased landscape complexity from 0.69 to 0.55 (p < 0.001), RES treatments increased landscape complexity to 0.84 (p < 0.001) (Fig. 5c). In addition, FHR treatments increased relative contagion index from 0.39 to 0.51 (p < 0.001), while RES decreased landscape contagion index to 0.21 (p < 0.001) (Fig. 5d).

Fig. 5
figure 5

Changes in modeled ecological variables after simulating conservation strategies on untreated (UNT) landscapes using a fuel hazard reduction strategy (FHR), and a restoration strategy (RES). Bars sharing a letter do not differ significantly at the α = 0.05 level after Bonferroni–Holm correction.

Full size image

FHR and RES treatments consistently differed in their impacts on forest structural heterogeneity, but treatment effects on fire behavior and post-fire erosion potential were similar. Integrated fire hazard index decreased from 102 in untreated catchments to 37.9 in catchments treated with either FHR and RES (p < 0.001), yet there was no difference in percent cover of high severity fire between treatments (p = 0.280) (Fig. 5e). In a similar pattern, treatments reduced the largest high severity patch index from 0.798 to 0.078 (p < 0.001), but there was no difference between maximum high severity patch size between FHR and RES treatments (p = 0.809) (Fig. 5f). Surface fire potential increased from 40.5 to 73.4% in catchments treated with either FHR or RES (p < 0.001), yet there was no difference between treatment types (p = 0.984). Expected post-fire sediment delivery to streams decreased from 25.5 to 13.3 Mg ha−1 in the FHR treatment (p < 0.001), and to 15.1 Mg ha−1 in the RES treatment (p < 0.001).

Tradeoffs among ecological responses

We found significant correlations among changes in many simulated ecological outcomes indicating both contrasts and similarities in effects of conservation practice designs (Fig. 6 and Fig. S2). With greater decreases in canopy cover, we found increases in gap cover and landscape complexity, and greater potential for prescribed fire use. Furthermore, greater decreases in canopy cover were correlated with decreases in relative contagion, integrated fire hazard, largest high severity patch index, and post-fire erosion. In general, changes among ecological outcomes were similar where desired outcomes increased together (e.g., reductions in canopy cover were correlated with reduced integrated fire hazard). However, the strength of the correlations differed and in some cases tradeoffs existed where one ecological benefit occured at the expense of another as indicated by several conservation practice design × variable interactions (Fig. 6). For example, RES treatments generally had greater reductions in canopy cover for a similar reduction in integrated fire hazard, large high severity patch size, and expected sedimentation (Fig. 6d, e and g) relative to FHR treatments. Other ecological variables showed a distinct tradeoff in FHR treatments but were complementary in RES treatments (Fig. 6o, q, and r). For example, with greater reductions in integrated fire index and post-fire sedimentation, landscape complexity decreased in FHR treatments, but increased with RES treatments (Fig. 6o and r). Similarly, increases in prescribed fire potential were correlated with increases in landscape complexity with RES treatments, yet landscape complexity was reduced with increasing surface fire in FHR treatments (Fig. 6q). These contrasts are important because they indicate that more spatially complex treatments meet landscape complexity objectives without strong tradeoffs in effectiveness of wildfire and erosion avoidance objectives.

Fig. 6
figure 6

Scatterplot matrix visualizing correlations among ecological response variables with each point representing the change in an ecological response in untreated catchments and treated catchments for FHR and RES. Plots with linear fits indicate significant correlations between changes in ecological variables. When interactions between response variables and conservation practice designs were significant, linear fits are shown separately for FHR treated catchments (in green) and RES treated catchments (in blue). Axis values are omitted for simplification but are included in Supplemental material (Fig. S2).

Full size image

Ecological effects differ spatially by treatment design, which has implications for prioritizing the type and location of treatments. As an example, changes in forest spatial complexity, integrated fire hazard, and expected sedimentation (post-fire) for catchments of the Upper South Platte watershed are shown in Fig. 7. FHR treatments generally reduce forest complexity across the watershed, while RES treatments increase complexity, with the greatest gains clustered in southern portions of lower montane forests (Fig. 7a and b). Although study-wide changes in integrated fire hazard were similar among treatments (Fig. 5e), RES treatments showed larger reductions (more negative) in lower montane forests and were less impactful (less negative) in upper montane forests compared to FHR treatments (Fig. 7c and d). This trend was also observed for reductions in sedimentation (Fig. 7e and f). Lastly, an example prioritization based on equal weighting among these three ecological variables shows that spatial prioritization differs among treatments. In FHR treatments, the largest effects can be gained in the northern portions of the watershed across both upper and lower montane regions (Fig. 7g). However, for RES treatments, the most impactful catchments for increasing potential benefits are located in the southern and western portions of the lower montane region (Fig. 7h).

Fig. 7
figure 7

Changes in modeled ecological effects following simulated FHR (a, c, e) treatments and RES treatments (b, d, f) for landscape complexity, integrated fire hazard, and expected post-fire sedimentation. Prioritization map (g, h) indicating percentile rank for where catchments with the greatest modeled benefit summed across three ecological indicators. Expected sedimentation (e, f) is shown in Mg ha−1, prioritization (g, h) is shown as percentile rank with a score of 1 indicating where greatest benefits can be achieved based on equal weighting of the three ecological indicators. Landscape complexity (a, b) and integrated fire hazard (c, d) are unitless.

Full size image

Discussion

Unsurprisingly, our modeling results suggest that simulated restoration treatments led to the greatest increases in measures of spatial heterogeneity thought to be important for animal and plant biodiversity, age class diversity, and resilience to disturbance. However, both the fuel hazard reduction and the restoration treatments achieved similar fire hazard and watershed risk outcomes. Changes in forest structural heterogeneity were markedly different between modeled conservation scenarios. RES treatments showed greater reductions in canopy cover, and greater effectiveness in enhancing spatial complexity with larger increases in gap cover, landscape complexity, and reductions in contagion (Fig. 5a–d). One striking finding was that RES treatments increased gap cover threefold, whereas gap cover was unchanged in FHR, treatments (Fig. 5b). We note that our definition of gaps allowed for areas with low canopy cover (<10%). However, FHR treatments reduced canopy cover proportionately, such that canopy removal was modest in areas with the potential to reduce canopy cover below the 10% threshold. By contrast, RES treatments reduced canopy cover by as much as 35% in dry areas with high canopy cover. Proportional reductions in FHR treatments tended to homogenize forest structure by reducing areas with high canopy cover by higher absolute amounts than areas with low canopy cover. This led to reductions in landscape complexity and increases in contagion, whereas RES treatments retained many areas of high cover while simultaneously increasing abundance of low cover, leading to higher complexity and lower contagion. These findings are somewhat unsurprising as the restoration algorithm was designed to create such effects at the landscape scale. However, it is valuable to understand what unique effects these treatments have compared to FHR. Several ecological attributes of interest are expected to increase with landscape heterogeneity of ponderosa pine forests, such as understory plant diversity and use by some wildlife species (Matonis and Binkley 2016; Stevens et al. 2016). Thus, it is important to emphasize spatial heterogeneity in project planning, as uniform FHR treatments will not generally achieve these objectives. The spatial metrics measured in this study are a first step towards linking conservation practice design to spatially driven biodiversity outcomes, but there is a recognized need for more research to make quantitative estimates of these outcomes.

Outcomes related to fire hazard and post-fire sedimentation were similar between the two conservation approaches examined here. We found that both approaches similarly attained goals to reduce fire hazard, diminish large patches of potential high severity fire, and increase prescribed fire suitability (Fig. 5e and f). By contrast, Stevens et al. (2016) found that RES treatments were less effective at achieving fire risk reductions relative to FHR treatments. This difference may be attributable to the relatively aggressive canopy removals specified for lower montane forests in our study based on historical conditions for the Colorado Front Range and the fact that all areas within catchments were considered for treatment based on historical conditions. Although reductions in post-fire sedimentation were not as sizable in RES treatments (−41%) relative to FHR treatments (−48%), expected post-fire sedimentation exhibited considerable reductions with both approaches. Factors enhancing the risk of post-fire sedimentation are high fuel loads, steep slopes, and close proximity to streams. The RES approach reduces mean canopy cover overall, but generally retains relatively higher canopy cover specifically in areas with higher soil moisture (e.g., north facing slopes) consistent with restoration guidelines (e.g., Addington et al. 2018), which may constrain treatment intensity in the highest erosion risk locations. Together, these results indicate that restoration treatments that explicitly incorporate spatial objectives are similarly effective at achieving objectives related to wildfire hazard reduction, without the landscape homogenizing effects found for treatments without explicit spatial objectives.

Modeling limitations

Though useful for providing information for landscape-scale decision making, there are several limitations to landscape-scale effects assessment. Restoration guidelines for dry conifer forests emphasize goals of generating open-canopied forest with high variability (e.g., Addington et al. 2018). However, the extent to which these goals are achieved in practice varies widely, in part due to constraints of land ownership, management restrictions, accessibility, operability, and economics (Underhill et al. 2014; North et al. 2015; Colavito 2017). These factors should be considered with our estimates of treatment benefits to identify high priority areas that are suitable for treatment. This framework offers a flexible method of specifying desired conditions as treatment input targets (e.g., Table 1) to allow evaluation of a wider spectrum of possible outcomes. Second, for simplification, our treatment simulation assumes that conservation practices are implemented over entire catchments averaging approximately 500 ha in size. Although this size is somewhat large for a single restoration treatment, it is within the size range of long-term planning units and thus offers insight to prioritization and planning at larger temporal and spatial scales. Third, our analysis relies on a series of coupled models (Fig. 2) and may be subject to multiple model and model linkage uncertainties. We used fire modeling software that is routinely used in hazard assessment and fuel treatment planning (Finney et al. 2019). The watershed model made reasonable predictions of post-fire sediment yields. None of our catchment mean hillslope erosion predictions for the first year after fire (range 0.9–68.1 Mg ha−1) exceeded the first year Buffalo Creek Fire average of 72.0 Mg ha−1 (Moody and Martin 2001; assumed 1.6 Mg m−3 bulk density) and the study wide mean of 23.4 Mg ha-1 was is close to the 22.0 Mg ha−1 observed in the 1st year after the Hayman Fire (Robichaud et al. 2013a). The predicted 1st-year post-fire catchment mean net sediment load to streams averaged 12.2 Mg ha−1 for the study catchments and only the upper quartile of study catchments (range 16.2–36.7 Mg ha−1) approached the higher annual catchment sediment yields of 22.2–38.6 Mg ha–1 from the Hayman Fire (Robichaud et al. 2008, 2013a). The modular nature of the assessment framework (Fig. 2) allows for integration of novel and process-based models as they develop to further improve understanding of conservation impacts. Fourth, our framework incorporates current burn probability (Short et al. 2016) into our estimate of integrated fire hazard and high severity patch index, but it does not capture changes in local and neighborhood burn probability from treatment. Studies that estimate changes in burn probability as a result of treatment are limited to a small set of scenarios due to the computational expense of burn probability models which account for fire spread among neighboring cells (Ager et al. 2007, 2016; Cochrane et al. 2012; Thompson et al. 2016). We chose not to account for burn probability effects in this study because it would be computationally impractical to evaluate these effects independently for each catchment, let alone combinations of treated catchments. Future expansions of the effects assessment framework would benefit from representation of fire spread processes.

Conclusions

Mapping modeled changes in ecological indicators reveal how spatial prioritization of treatments may vary when considering multiple outcomes and scenarios (Fig. 7). Overall impacts of fire behavior and erosion effects were similar among the treatment approaches simulated here (Fig. 5e and h). However, the treatment locations of greatest benefit varied depending on the ecological indicator and treatment approach under consideration. In some cases, resource planners can utilize landscape scale treatment simulations in developing alternatives for prioritization (Huayhuaca 2016; Addington et al. 2018; Gannon et al. 2019). Our results suggest that both the ecological effects of restoration as well as the most effective places for treatment vary based on the conservation practices used. We note that our example prioritization based on equal weighting of three variables does not necessarily reflect agency and public preferences for achieving these outcomes, which should be incorporated into a multi-resource prioritization.

A number of landscape assessments have been developed to evaluate the potential ecological impacts of various conservation scenarios that differ in the treatment types simulated, ecological impacts modeled, and geography considered. Our assessment framework was influenced heavily by commonly used wildfire risk assessment methods (Scott et al. 2013; Thompson et al. 2013a), which integrate probabilistic estimates of fire likelihood and intensity with relative response functions for a broad range of resources. Extensions of this framework have focused on improving rigor of the effects assessment by replacing relative response functions with quantitative models and differencing baseline and post-treatment risk to estimate management effects. For example, Gannon et al. (2019) linked potential wildfire behavior to quantitative erosion and sedimentation models to estimate the avoided water supply impacts from mechanical thinning and/or prescribed fire. Stevens et al. (2016) conducted simulations of different landscape configurations, including treatments designed to uniformly reduce fuel hazards and protect the wildland-urban interface, along with other complex restoration treatments that varied by topography, and showed how the different configurations impacted a number of processes (wildlife diversity, fire behavior, smoke production). Here we build on these previous efforts by (1) developing a restoration simulation algorithm that is flexible to parameterize a wide range of conditions specified by canopy cover distribution (Fig. 3, Table 1), and (2) quantifying non-fire-related metrics of landscape heterogeneity to explore potential tradeoffs with other ecological processes. We include several categories of ecological effects related to common objectives for dry forest restoration including wildlife habitat improvement, fire hazard reduction, and watershed protection. Assessment frameworks that offer flexibility to simulate a wide range of treatment types and scenarios on multiple ecological metrics can facilitate comparison of potential treatment strategies for planning and prioritization.