Historically, fire has been essential in Southwestern US forests. However, a century of fire-exclusion and changing climate created forests which are more susceptible to uncharacteristically severe wildfires. Forest managers use a combination of thinning and prescribed burning to reduce forest density to help mitigate the risk of high-severity fires. These treatments are laborious and expensive, therefore optimizing their impact is crucial. Landscape simulation models can be useful in identifying high risk areas and assessing treatment effects, but uncertainties in these models can limit their utility in decision making. In this study we examined underlying uncertainties in the initial vegetation layer by leveraging a previous study from the Santa Fe fireshed and using new inventory plots from 111 stands to interpolate the initial forest conditions. We found that more inventory plots resulted in a different geographic distribution and wider range of the modelled biomass. This changed the location of areas with high probability of high-severity fires, shifting the optimal location for management. The increased range of biomass variability from using a larger number of plots to interpolate the initial vegetation layer also influenced ecosystem carbon dynamics, resulting in simulated forest conditions that had higher rates of carbon uptake. We conclude that the initial forest layer significantly affects fire and carbon dynamics and is dependent on both number of plots, and sufficient representation of the range of forest types and biomass density.

Initial Communities Data

The initial communities layer is the base vegetation layer that sets the starting conditions for the exchange of carbon, water, energy, species interactions, disturbance effects, and other landscape processes. The initial treatment optimization study in this landscape (Krofcheck et al., 2019) used 68 Forest Inventory and Analysis (FIA) plots from within the Santa Fe National Forest that had been inventoried in 2010 or later and had not burned since 2005. Forest types represented by the FIA plots were piñon-juniper, ponderosa pine, Douglas-fir (Pseudotsuga menziesii), Engelman spruce (Picea engelmannii) and limber pine (Pinus flexilis). The latter three were grouped into a general mixed-conifer forest type. The authors then used elevation, transformed aspect using Topographic Radiation Aspect Index, TRASP (Roberts & Cooper, 1989), and a tasseled cap transformation of spectral data from Landsat 8 (available at https://www.usgs.gov/landsat-missions/landsat-8) as predictors for Random Forest models and used the rfUtilities library (Evans & Murphy, 2018) to select the most parsimonious model for each forest category separately. Existing vegetation classification from the Southwest Regional Gap Analysis (SWReGap, https://swregap.org/) using the ‘yaImpute’ library (Crookston & Finley, 2008) was used to stratify the measured plots for the imputation. We determined plot sampling intensity by calculating the relative area of land each plot represents within its forest type.

To evaluate the influence of additional plot data on the initial communities layer and its effects on model behavior, we used data collected as part of the planning process by the US Forest Service. These data were located entirely within the study area and included 1072 plots from 111 stands inventoried in 2011, where each stand included between 3 and 31 plots. Plot data were collected using a common stand exam protocol using variable radius plots. The specific Basal Area Factor (BAF) prism chosen for each stand was a function of stand density and they ranged from 10 to 30 BAF. We had coordinates for the centroid of each stand, but not for each individual plot, which effects the imputation process. We used the tree data from each plot to determine a specific forest type and generalized category (e.g. piñon-juniper, ponderosa pine, mixed-conifer), corresponding to the FIA classification, and added an aspen forest type, resulting in a total of four generalized forest types. We also defined non-forested areas and included two generalized species parameterizations to represent shrubs that resprout and shrubs that do not resprout following fire. We used all FIA and common stand exam (CSE) plot data (n=1140) to generate a new initial communities layer following the same method as Krofcheck et al. (2019) in R v4.1.2 (R Core Team, 2021).

Simulation Analysis

To estimate the uncertainty in the initial communities layer that is due to not having coordinates for all CSE plots, we re-ran interpolations randomly selecting one plot for each set of stand coordinates. This led to 31 initial communities layers, which we used to initialize the model with the five climate projections, for a total of 155 simulations. We compared the aboveground carbon following model initialization of these initial communities layers with the initial communities layer that we created using all plot data and that we used for our management simulations. We calculated the difference in aboveground carbon between each layer and the one we used in our simulations to determine how much the initial communities layer is influenced by this source of uncertainty.

To determine the influence of the number of plots used in the development of the initial communities layer, we produced five additional initial communities layers with differing numbers of underlying plot data. For four of the five layers, we halved the number of CSE plots used in the interpolation each time (e.g. 536, 268, 134, 67) and combined those with the FIA data. For the fifth initial communities layer, we only used the 68 FIA plots. For each of the layers, we randomly selected plots from each forest type proportional to the prevalence of each forest type on the landscape. We then initialized the model using each of these initial communities layers using the five climate projections and compared the aboveground carbon following model initialization to that of the initial communities layer that we used for our simulations.

We quantified differences between our primary initial communities layer and that of Krofcheck et al., (2019) by comparing the difference in quantity and distribution of aboveground carbon at the beginning and end of the simulations. We ran an independent t-test to assess the difference in carbon between the two studies at each site every 10 years for each of the climate models, and computed the percent of area with a significant difference (p < 0.01) in aboveground carbon. We compared treatment location as determined by the probability of high-severity fire between our initial communities layer and that of Krofcheck et al. (2019). We calculated Net Ecosystem Carbon Balance (NECB) by subtracting carbon lost from the system (treatment and wildfires) from carbon gained (photosynthesis) and then relativized the treatment scenario NECB values to the no-management scenario for both our simulations and those of Krofcheck et al. (2019). Data processing and analysis was conducted using R v4.1.2 (R Core Team, 2021).

Treatment scenarios

To develop the optimized treatment placement scenario, we first ran simulations that included no management to identify locations where landscape conditions were such that there was a high probability of high-severity wildfire. We ran the no-management simulations using the same five projected climate data sets and fire weather data described above. We ran 25 replicate simulations using each of five projected climate data sets, for a total of 6250 simulation years. We used fire severity raster data from these model outputs to quantify the probability of high-severity by dividing the number of years with high-severity fires by the total number of fire years per site. We then identified sites with a probability of high-severity fire greater than 0.3 and targeted those locations in the treatment scenario simulations, assigning treatment to those areas first. 

To determine the type of treatment we used the probability of high-severity fire in combination with slope and forest type. We limited our management simulations to the ponderosa pine and dry mixed-conifer forest where the combined ponderosa pine and Douglas-fir aboveground carbon was at least 65% of the total. We used the same thinning and prescribed burning treatments as Krofcheck et al. (2019), which were designed to approximate common treatments for the region. Thinning treatments simulated thinning from below by removing approximately 30% of the biomass, preferentially removing the youngest cohorts (Hurteau et al., 2011, 2016) and was only applied to ponderosa pine forest and confined  to slopes <30%, to account for a common limitation on mechanical thinning. We simulated prescribed burning based on historic mean fire return intervals, with all ponderosa pine burned using a 10-year return interval and forests co-dominated by ponderosa pine and Douglas-fir burned using a 15-year return interval. The forest type and probability of high-severity fire are highly dependent on the initial communities layer which defines the initial forest conditions. As a result, our treatment placement map differed substantially from the one in Krofcheck et al. (2019).

To examine the effects of the treatment on the landscape we produced a new probability of high-severity fires raster and calculated the difference in aboveground carbon between the management and no management scenarios of this study at the end the simulations.

We ran simulations over a 50-year period, using climate model projections for years 2000-2050. We ran 25 replicates for each of the five climate projections, totaling 125 simulations each for the no-management and management scenarios. Fire weather distributions tracked projected climate and were updated each decade to account for changes in temperature and precipitation.

R