A changing climate is altering ecosystem carbon dynamics with consequences for natural systems and human economies, but there are few tools available for land managers to meaningfully incorporate carbon trajectories into planning efforts. To address uncertainties wrought by rapidly changing conditions, many practitioners adopt resistance and resilience as ecosystem management goals, but these concepts have proven difficult to monitor across landscapes. Here, we address the growing need to understand and plan for ecosystem carbon with concepts of resistance and resilience.  Using time series of carbon fixation (n=103), we evaluate forest management treatments and their relative impacts on resistance and resilience in the context of an expansive and severe natural disturbance. Using subalpine spruce-fir forest with a known management history as a study system, we match metrics of ecosystem productivity (net primary production, g·C·m²·y^-1) with site-level forest structural measurements to evaluate (1) whether past management efforts impacted forest resistance and  resilience during a spruce beetle (Dendroctonus rufipennis) outbreak, and (2) how forest structure and physiography contribute to anomalies in carbon trajectories. Our analyses have several important implications. First, we show that the framework we applied was robust for detecting forest treatment impacts on carbon trajectories, closely tracked changes in site-level biomass, and was supported by multiple evaluation methods converging on similar management effects on resistance and resilience. Second, we found that stand species composition, site productivity, and elevation predicted resistance, but resilience was only related to elevation and aspect. Our analyses demonstrate application of a practical approach for comparing forest treatments and isolating specific site and physiographic factors associated with resistance and resilience to biotic disturbance in a forest system, which can be used by managers to monitor and plan for both outcomes. More broadly, the approach we take here can be applied to many scenarios, which can facilitate integrated management and monitoring efforts.

Remote sensing procedures

We used multitemporal NPP trajectories to develop our resilience and resistance metrics. The Numerical Terradynamic Simulation Group at the University of Montana has developed a vegetation productivity product that estimates yearly vegetation productivity at 30m resolution for the contiguous US from 1986 to 2018 (Robinson et al. 2018). We downloaded the annual Landsat NPP time series data from 1986 to 2018 across the landscape of interest that experienced a spruce beetle outbreak from 2008-2015. The product includes annual NPP and 16-daily gross primary productivity (GPP) estimates. GPP (gC/m²/y^-1) is derived using a light use efficiency model based on meteorological, Landsat data, and biome parameter lookup tables, following this equation:

GPP = LUE_max × f_Tmin × f_vpd × 0.45 × SW_rad × FPAR                         (Eq. 1)

 LUE is the biome specific maximum light use efficiency (in gC/MJ), f_Tmin and f_vpd are biome-specific temperature and vapor pressure deficit down regulators. SW_rad is the incoming shortwave radiation (45% of that is available for photosynthesis) and FPAR is the fraction of photosynthetically active radiation, which is derived from the Landsat imagery. This GPP product is calculated for every 16-day period throughout the year. Annual NPP is calculated as the total sum of the 16-day GPP estimates within a given year minus the autotrophic respiration (Robinson et al. 2018). As the Landsat data record is inherently noisy due to atmospheric effects, sensor differences, clouds, and retrieval errors, we adopted the methods of Robinson et al. (2018) to perform a gap-filling method, employing climatology and smoothing approaches to harmonize the satellite data record (for greater detail see Robinson et al. 2017).

To assess the resistance and resilience functions of the plot locations following spruce beetle disturbance, we developed NPP anomalies from an undisturbed reference mean. First, we randomly selected a set of 60 locations in the vicinity of our plot locations (i.e., southcentral Colorado) that did not show any mortality. These undisturbed locations were selected using aerial survey data (i.e., selection of locations that were outside areas identified as affected by the spruce beetle) and placed outside known management areas, but within Engelmann spruce dominated forest as determined by the Individual Tree Species Parameter Maps (U. S. Forest Service 2020). Next, we plotted these reference sites over time and removed plots (n=35) that showed a deviation from a temporal mean at the end of the time series; such locations likely experienced other unaccounted disturbances (such as harvest) or experienced spruce beetle mortality but might have been omitted by the aerial surveys.  

To compare sites and calculate our resilience and resistance metrics (Figure 1), we calculated an NPP anomaly at the pixel-level, i.e., the difference of NPP from a temporal mean. The NPP anomaly (dNPP) for each plot was then calculated as follows:                                                         

dNPP_{plot, year} = NPP_{plot, year} - Mean(NPP_undisturbed_year)                               (Eq.2)

 Where NPP_{plot, year} was the NPP (gC/m²/y^-1) for a given year (1986 - 2018) for a given plot and Mean(NPP_undisturbed_{plot, year}) was the mean NPP for the undisturbed plots for a given year. We corrected for differences in productivity between disturbed and undisturbed plots by adjusting the anomalies with the mean difference of NPP over the five years before the outbreak event. We compared the maximum NPP reduction (dNPP_max) with the reduction of the percent biomass from field observations, and this comparison verified that the NPP data product was sensitive to the effects of the spruce beetle outbreak and reflected changes in carbon fixation due to loss of photosynthetic biomass (R²=0.804, RMSE= 1195 gC/m²/y^-1; Appendix S1: Figure S1).

Following the establishment of the anomalies, we calculated the maximum decrease in dNPP in a single year after the initial outbreak; we used this metric as the resistance metric, i.e., the maximum reduction of productivity (dNPPmax). For most stands this occurred between 2011-2012. Therefore, dNPPmax reflects the degree to which the ecosystem departs from equilibrium due to a disturbance (in this case, a spruce beetle outbreak). Subsequently, we determined the year at which the maximum reduction of NPP occurred and fitted a linear model of dNPP across the years from the maximum reduction of NPP to the final year in the time series. The slope of that fitted line was analyzed as the resilience metric (Figure 1). In other words, we treat resilience as the rate (gC/m2/y2) at which ecosystem processes recover from disturbance towards the prior steady state. These metrics where then used to assess ecosystem resistance and resilience for the plot locations and we explored how different management, environmental, and topographical characteristics affected the resilience and resistance metrics.