Landscape-scale analysis of aboveground tree carbon stocks affected by mountain pine beetles in Idaho

Bark beetle outbreaks kill billions of trees in western North America, and the resulting tree mortality can significantly impact local and regional carbon cycling. However, substantial variability in mortality occurs within outbreak areas. Our objective was to quantify landscape-scale effects of beetle infestations on aboveground carbon (AGC) stocks using field observations and remotely sensed data across a 5054 ha study area that had experienced a mountain pine beetle outbreak. Tree mortality was classified using multispectral imagery that separated green, red, and gray trees, and models relating field observations of AGC to LiDAR data were used to map AGC. We combined mortality and AGC maps to quantify AGC in beetle-killed trees. Thirty-nine per cent of the forested area was killed by beetles, with large spatial variability in mortality severity. For the entire study area, 40–50% of AGC was contained in beetle-killed trees. When considered on a per-hectare basis, 75–89% of the study area had >25% AGC in killed trees and 3–6% of the study area had >75% of the AGC in killed trees. Our results show that despite high variability in tree mortality within an outbreak area, bark beetle epidemics can have a large impact on AGC stocks at the landscape scale.


Introduction
North American forests are important components of the global carbon cycle, commonly acting as carbon sinks in which net primary production (NPP) exceeds heterotrophic respiration (R h ) (CCSP 2007). Forest disturbances have substantial effects on carbon cycling through the reduction of NPP and increase in R h following tree mortality, and effects can last decades (Odum 1969). In North America, biotic disturbances (insect and disease outbreaks) can cause substantial tree mortality (e.g. Meddens et al 2012). As a result, biotic disturbance agents can have significant impacts on local-to regional-scale carbon cycling . Mountain pine beetle (Dendroctonus ponderosae Hopkins) outbreaks are major biotic disturbances that affect carbon cycling (Kurz et al 2008). Trees killed by mountain pine beetles progress from green to red in the summer following attack, then lose their needles three to five years after attack, becoming gray (Wulder et al 2006). NPP of beetle-killed forests is reduced for several years (Pfeifer et al 2011, Brown et al 2010, and decomposition is enhanced with the increase in dead organic matter, causing a decrease in net ecosystem productivity and often resulting in a conversion from a carbon sink to a carbon source (Kurz et al 2008, Edburg et al 2011. Analyses of bark beetle outbreak impacts on carbon stocks and fluxes have been done on plot and stand scales and have found substantial ranges in mortality severity and resulting impacts Romme et al 1986, Pfeifer et al 2011, Edburg et al 2011, Brown et al 2010. For example, Pfeifer et al (2011 found that reductions in aboveground carbon (AGC) stocks and fluxes varied from 30% to 80% among 12 plots sampling one outbreak. Such variability suggests considerable uncertainty in the aggregated impact within one outbreak area, yet such information is needed for quantifying larger-scale (landscape to regional) impacts of bark beetle outbreaks. A landscape-scale study using remote sensing to document spatial variability of bark beetle-affected AGC stocks has not been done to our knowledge.
Remotely sensed imagery is an excellent tool for evaluating processes with high spatial variability. In a previous study, we demonstrated the utility of combining multispectral (to map tree mortality) and LiDAR remote sensing to quantify AGC stocks of trees affected by bark beetle outbreaks (Bright et al 2012). In this study, we applied this methodology to a 5054 ha study area in central Idaho, USA. Our objective was to perform a landscape-scale analysis of AGC stocks after a bark beetle infestation by (1) assessing spatial variability in tree mortality and AGC in trees within the study area and (2) estimating spatially aggregated impacts to AGC stocks across the study area.

Methods
Methods to classify tree mortality, estimate AGC stocks, and quantify AGC in trees killed by mountain pine beetles are described in detail in Bright et al (2012). Here we provide a brief summary of these methods.

Study location
The study area was a 5054 ha, predominately lodgepole pine (Pinus contorta var. latifolia Dougl. ex Loud.) forest in central Idaho, USA (44.3 • N, 115.1 • W) that has experienced a bark beetle outbreak since 2002. Field observations indicated that the outbreak was subsiding in the summer of 2010. Other tree species, which include Douglas-fir (Pseudotsuga menziesii Franco), subalpine fir (Abies lasiocarpa Nutt.) and Engelmann spruce (Picea engelmannii Parry ex Engelm.), also occurred in limited numbers. We selected a relatively flat study area to minimize errors that can occur as a result of sloping topography.

Field observations
The study area was stratified into low, medium, and high normalized difference vegetation index (NDVI) classes using a pre-outbreak Landsat TM 5 NDVI image. Nine plot locations were randomly located within each stratum. We established 0.04 ha (400 m 2 ) fixed-radius plots in the summer of 2010 and measured the location, diameter at breast height (DBH), height, crown diameter, health (live or dead), and canopy dominance of each tree >7 cm DBH and ≤ 11.3 m from plot center. We computed canopy cover within plots by summing crown areas of dominant trees. Following Standish et al (1985), we used allometry to predict aboveground biomass from lodgepole pine DBH and height measurements.
We then multiplied aboveground biomass estimates by 0.5 to estimate AGC content of each tree (Lamlom and Savidge 2003).
We excluded areas within 91 m (300 ft) of roads because firewood cutting was allowed within these areas, and we wished to minimize confounding effects from known human impacts. Areas of the study area >2150 m in elevation were also excluded from analysis because forest composition is less dominated by lodgepole pine above 2150 m. These areas were excluded from both plot sampling and mapping.

Mapping mortality and aboveground carbon with remote sensing
Watershed Sciences, Inc. acquired four-band (blue, green, red, near-infrared), 0.2 m digital aerial imagery across the study area. We aggregated the imagery to 2.4 m spatial resolution to maximize classification accuracy . Locations with vegetation less than 0.5 m high, as determined from LiDAR data, were excluded from classification, as well as areas of water and shadow. Field observations guided the on-screen selection of green, red and gray tree pixels on true-color imagery. Selected pixels were used to train and evaluate a maximum likelihood classification. Overall classification accuracy was 87.2% with a kappa of 0.79 (Bright et al 2012).
We then applied the classifier to the study area to map locations with green, red, and gray trees. We aggregated the 2.4 m mortality map to a 20 m per cent mortality map to overlay on the 20 m plot-level AGC map, and also aggregated the map to 1 ha grid cells to assess variability in mortality.
High-density (8.7 returns m −2 ) LiDAR data was also acquired across the study area. We normalized the intensity values for range and automatic gain control (BCAL LiDAR Tools 2011, Korpela et al 2010), and classified returns as ground and non-ground (Evans and Hudak 2007). Groundclassified returns were interpolated to produce a bare-earth digital terrain model (DTM), which was subtracted from the LiDAR point cloud to calculate the height of each non-ground return (McGaughey 2008). LiDAR metrics of returns within each tree and plot extent were produced with FUSION software (McGaughey 2008), and these metrics were related to field-observed AGC estimations using models developed in R (R Development Core Team 2011). We developed both treeand plot-level models relating AGC to LiDAR metrics (Bright et al 2012). Correlation coefficients of observed versus predicted AGC for green, red and gray tree-level generalized additive models (GAMs) using cross-validation were 0.53, 0.50 and 0.54, respectively, with root mean square error (RMSE) values of 28.5, 26.3 and 42.1 kg AGC (48%, 32% and 39%), respectively (Bright et al 2012). The correlation coefficient of observed versus predicted plot AGC for the cross-validated plot-level regression model was 0.84, with an RMSE of 9.2 Mg AGC ha −1 (12%) (Bright et al 2012).
We then applied the three tree-level models to 2.4 m LiDAR metric grids to predict AGC across the study area using the AsciiGridPredict() function of the yaImpute package (Crookston and Finley 2007) in R. Tree-level ln(C) maps were back transformed to AGC and corrected for back-transformation bias (Snowdon 1991). Similarly, we applied plot-level models to 20 m LiDAR metric grids to predict AGC at the plot resolution. Confidence intervals for maps were derived by multiplying landscape AGC estimates by %RMSE values of cross-validated models. We combined mortality maps with AGC maps to produce maps of AGC in killed trees at 2.4 m and 20 m spatial resolutions. We also produced maps of %AGC in killed trees at 1 ha spatial resolutions.

Field observations
Overall, 1222 out of 1379 trees measured were lodgepole pine (89%), of which 408 (33%) were killed by bark beetles (30% of all trees). Dead canopy cover averaged 38% across plots and ranged from 7% to 83%. AGC density of plots averaged 77 Mg AGC ha −1 and ranged from 27 to 121 Mg AGC ha −1 . Per cent AGC in beetle-killed trees averaged 41% across plots and ranged from 13% to 85%. Lodgepole pine-dominated plots with greater AGC experienced greater mortality in terms of %AGC in killed trees, although the relationship was weak (R 2 = 0.14). Average AGC content of green, red and gray lodgepole pine trees measured in the field was 46 ± 34, 76 ± 40 and 86 ± 65 kg AGC, respectively. Green lodgepole pines were smaller and hence contained significantly less AGC per tree than either red (p < 0.001) or gray trees (p < 0.001, Wilcoxon rank-sum test). Red trees did not contain significantly less AGC per tree than gray trees (p = 0.846, Wilcoxon rank-sum test).

Tree mortality
Tree mortality was positively spatially autocorrelated across the study area (p < 0.0001, Moran's I); i.e., areas of high mortality tended to be located adjacent to other areas of high mortality, suggesting clumped mortality ( figure 1(a)). Generally, the study area contained four areas of severe mortality in the northwest corner, center, southwest corner, and southeast panhandle of the study area. Red trees, indicating more recent beetle activity, were generally located in the northern part of the study area. The western edge of the study area, where elevation is greatest, was less affected by bark beetles because forest composition is mixed at higher elevations, thereby including nonhost tree species.

Aboveground carbon
Like tree mortality, AGC was positively spatially autocorrelated across the study area (p < 0.0001, Moran's I). Generally, the study area contained three areas of high AGC content: the north corner, center and southwest edge. AGC was lower in the southern and eastern part of the study area than the central and northern part.
The 2.4 m AGC map predicted a total of 0.295 ± 0.127 Tg AGC (58 ± 25 Mg AGC ha −1 ) for the study area. The 20 m AGC map estimated a total of 0.288 ± 0.03 Tg AGC (57 ± 6.8 Mg AGC ha −1 ), about 2% less AGC than the 2.4 m AGC map (table 1). The confidence interval of the 20 m AGC map was smaller (±12%) than the confidence interval of the 2.4 m AGC map (±43%) due to the higher %RMSE values of the tree-level models.

Aboveground carbon in killed trees
Many areas of high mortality coincided with areas of high AGC content. The north corner, center and southwest corner Table 1. Beetle-caused tree mortality area and aboveground carbon (AGC) for the entire study area as predicted by the 2.4 and 20 m AGC maps. Confidence intervals were derived by multiplying totals by per cent root mean square errors from cross-validated models. Because the green, red, and gray components of the 20 m AGC map were not derived in the same manner, confidence intervals for these components could not be calculated (see text for more details).  of the study area had both high AGC content and high mortality. The southeast part of the study area, where AGC density was lower, was less affected by bark beetles. The 2.4 m map predicted a total of 0.146 Tg AGC (50%, 29 Mg AGC ha −1 ) in beetle-killed trees, whereas the 20 m map predicted a total of 0.116 Tg AGC (40%, 23 Mg AGC ha −1 ) in beetle-killed trees (table 1). For comparison, total (landscape) AGC in green trees was 0.149-0.172 Tg AGC and 29.5-34.0 Mg AGC ha −1 (2.4 and 20 m maps). Per cent AGC in killed trees of 1 ha cells averaged 47% and 38% for aggregated 2.4 m and 20 m AGC maps, respectively (figure 2). In the 2.4 m AGC map, 89% of hectares contained >25% AGC in killed trees, 45% of hectares contained >50% AGC in killed trees and 6% of hectares contained >75% AGC in killed trees. In the 20 m AGC map, 75% of hectares contained >25% AGC in killed trees, 22% of hectares contained >50% AGC in killed trees and 3% of hectares contained >75% AGC in killed trees (figure 2).

Discussion
Although a previous plot-scale study documented large variability in beetle-affected AGC stocks (Pfeifer et al 2011), that study did not show how this variability aggregated across a landscape (i.e., the net effect). Our spatially explicit study showed that landscape-scale impacts to AGC stocks in trees were large. This large impact occurred despite the fact that only a small percentage of the area was severely impacted; instead, the large spatial extent of moderately affected stands resulted in a large landscape modification. Accounting for variability in mortality across a region is crucial for understanding and modeling carbon cycling of these beetle-infested landscapes, as outbreak severity is an important factor determining the impact of bark beetle-caused mortality on carbon fluxes of beetle-infested forests (Pfeifer et al 2011, Edburg et al 2011. Mortality severity (%AGC in killed trees) was greater in areas of higher AGC content because these areas likely contain (1) larger trees, which bark beetles preferentially attack (Hopping and Beall 1948), (2) older trees, which can be more susceptible to successful attack (Safranyik et al 1974), and (3) densely growing trees, where low vigor can lead to vulnerability and successful attack (Shore and Safranyik 1992). Plot-and tree-level AGC maps estimated similar amounts of total AGC across the study area. Total AGC estimated by the plot-level AGC map was more certain than total AGC estimated by the tree-level AGC map because of greater uncertainties in tree-level models used to produce tree-level AGC maps. The 2.4 m AGC map estimate showed a greater percentage of AGC contained in beetle-killed trees than the 20 m AGC map. Although the 20 m map predicted total AGC more accurately and with smaller confidence intervals compared with plot observations, the 2.4 m AGC map predicted %AGC in killed trees with less bias because the 2.4 m AGC map accounted for differing AGC densities of green and beetle-killed pixels (beetle-killed pixels tended to contain more AGC because beetles preferentially attack larger trees). The 20 m AGC map did not account for differing AGC densities of green and beetle-killed trees. Thus, AGC contained in beetle-killed trees was underestimated by the 20 m AGC map and more accurately estimated by the 2.4 m AGC map (figure 2).
Our estimates of mortality and AGC in trees are similar to estimates of other studies. At the plot scale, our range of 7%-83% mortality is similar to mortality reported by plot-scale studies of mountain pine beetle in lodgepole pine (Romme et al 1986, Pfeifer et al 2011, Meddens et al 2009. Our estimates of AGC densities for our study area (plots averaged 77 Mg AGC ha −1 and ranged from 27 to 121 Mg AGC ha −1 ) are similar to AGC densities of other studies in similar forest types (Pfeifer et al 2011, Chatterjee et al 2009. At the landscape scale, our estimate of 39% mortality across the study area was similar to but less than 59% mortality within a 94 km 2 forested area in Colorado (Meddens et al 2009). Our landscape AGC density estimates (57-58 Mg AGC ha −1 ) are similar to those reported by Hicke et al (2007), who found densities of 50-80 Mg AGC ha −1 for interior forests of the Rocky Mountains, and Turner et al (1995), who found densities of 50-60 Mg AGC ha −1 for the Rocky Mountains and east side of the Cascades.
The transition of 40%-50% of AGC from live to dead trees will affect the net carbon flux of the study area. The primary fate of the carbon in killed trees is to return to the atmosphere via decomposition over the course of the coming decades. This decomposition rate is highly dependent on a number of factors that include snagfall rates and nutrient dynamics (Edburg et al 2011). Remaining and newly established vegetation also affect net carbon exchange with the atmosphere through stand recovery, and may experience increased growth rates because of the increased availability of light, water and nutrients (Romme et al 1986, Brown et al 2010. Rates of decomposition and growth govern whether this area is a net carbon source, which is expected early after such a disturbance, or a sink, which will likely occur later in the stand trajectory (Edburg et al 2011). Thus, despite a large pulse of dead organic matter, bark beetle outbreaks only temporarily reduce forest carbon sequestration (similarly to other disturbance types).
Our study assessed AGC in trees, which has some limitations for a more complete understanding of the carbon cycle. A more complete analysis of the impacts to shrubs, herbaceous vegetation, and soils is needed to quantify impacts to the total ecosystem carbon pool. Furthermore, additional research at landscape scales is needed to quantify how this tree mortality affects carbon fluxes as well as temporal dynamics of carbon stocks and fluxes following the disturbance and would provide valuable evaluation data for ecosystem models.

Conclusions
Our spatially explicit, landscape-scale study illustrates the large spatial variability of beetle-caused tree mortality and its effects on AGC stocks. This variability needs to be accounted for to accurately quantify the impacts of bark beetle outbreaks on carbon cycling when aggregating to broader spatial scales. Field sampling designs should sample the entire range of mortality severity to minimize model extrapolations. Use of high-resolution remote sensing to map bark beetle impacts provides spatial context to effects on AGC stocks that are directly observable at the landscape scale, which helps to upscale AGC effects to larger spatial scales using coarser resolution imagery, where bark beetle effects are harder to see. Most of our 5054 ha study area was affected by some level of mortality, but only a small area was severely impacted. Despite this variability, we found large impacts to AGC stocks when aggregated across the study area. We estimated that 39% of forest cover was killed by bark beetles and between 40% and 50% of AGC was in beetle-killed trees. Such a large reduction of live AGC stocks will substantially alter carbon sequestration in these forests in the short term. Because of the large impacts to carbon stocks, beetle effects should be included in carbon models for more accurate carbon accounting. Given the large extent of beetle epidemics in western North America (Meddens et al 2012), we suggest these forest disturbances have a significant impact on regional carbon cycling. Furthermore, climate change is expected to continue to facilitate bark beetle outbreaks in North America (Bentz et al 2010), suggesting continued effects on the capability of these forests to sequester carbon.