Simulated impacts of mountain pine beetle and wildfire disturbances on forest vegetation composition and carbon stocks in the Southern Rocky Mountains

Forests play an important role in sequestering carbon and offsetting anthropogenic greenhouse gas emissions, but changing disturbance regimes may compromise the capability of forests to store carbon. In the Southern Rocky Mountains, a recent outbreak of mountain pine beetle ( D ndroctonus ponderosae ; MPB) has caused remarkable levels of tree mortality. To evaluate the long-term impacts of both this insect outbreak and another characteristic disturbance in these forests, high-severity wildfire, we simulated potential changes in species composition and carbon stocks using the Forest Vegetation Simulator (FVS). Simulations were completed for 3 scenarios (no disturbance, actual MPB infestation, and modeled wildfire) using field data collected in 2010 at 97 plots in the lodgepole-pine-dominated forests of eastern Grand County, Colorado, which were heavily impacted by MPB after 2002. Results of the simulations showed that (1) lodgepole pine remained dominant over time in all scenarios, with basal area recovering to pre-disturbance levels 70–80 yr after disturbance; (2) wildfire caused a greater magnitude of change than did MPB in both patterns of succession and distribution of carbon among biomass pools; (3) levels of standing-live carbon returned to pre-disturbance conditions after 40 vs. 50 yr following MPB vs. wildfire disturbance, respectively, but took 120 vs. 150 yr to converge with conditions in the undisturbed scenario. Lodgepole pine forests appear to be relatively resilient to both of the disturbances we modeled, although changes in climate, future disturbance regimes, and other factors may significantly affect future rates of regeneration and ecosystem response.


Introduction
The global atmospheric carbon pool (estimated for 2003) holds about 705 petagrams of carbon (Pg C; 1 Pg = 10 15 g), of which 535 Pg C are from non-anthropogenic sources and 170 Pg C are from anthropogenic sources (King et al., 2007).Anthropogenic emissions from fossil-fuel combustion and land-use/land-cover change have been increasing; for the United States, recent emissions estimates suggest a 0.5 % annual rate of increase between 1990 and 2010 (US Environmental Protection Agency, 2012).Such rates of change incite concern because increases in atmospheric carbon are thought to be a primary driver of global climate change (Intergovernmental Panel on Climate Change, 2007).Some anthropogenic emissions are offset by ecosystem carbon sequestration, especially in forests, which cover about a third of Earth's land mass and sequester and store large amounts of carbon in soils, live biomass, and dead biomass.In North America, forests contain approximately half of the carbon stored in ecosystems and carbon sequestration by forests offsets 0.21 Pg C of the 1.68 Pg C emitted by natural and anthropogenic sources each year (King et al., 2007).Disturbances, such as wildfire and insect outbreaks, impact forests periodically, altering carbon stocks and sequestration rates as dead vegetation decomposes and newly established vegetation Published by Copernicus Publications on behalf of the European Geosciences Union.
grows.Understanding the role played by forested ecosystems and their associated disturbance regimes in constraining or contributing to carbon sequestration and offsetting anthropogenic carbon emissions into the future represents a fundamental challenge in global change research (Liu et al., 2011;Running, 2008;Goetz et al., 2012;Hicke et al., 2012a).
In recent years, an outbreak of the mountain pine beetle (Dendroctonus ponderosae; MPB) has resulted in extensive tree mortality in North American lodgepole pine (Pinus contorta) forests.Historically, MPB has persisted at endemic population levels that periodically erupted into large-scale outbreaks (Amman, 1977;Baker and Veblen, 1990;Raffa et al., 2008).However, the extent and severity of the current outbreak is remarkable, affecting nearly 3 700 000 ha in the conterminous United States by 2009 (Mann, 2012).There is debate and concern regarding the extent to which this outbreak will alter both short-and long-term wildfire hazard and risk (Jenkins et al., 2008;Bentz et al., 2010;Hicke et al., 2012b), influence carbon storage capacity (Kurz et al., 2008;Hicke et al., 2012a), and induce shifts in vegetative composition (Collins et al., 2011;Diskin et al., 2011).Given the unusual magnitude of this epidemic compared to previous outbreaks, few precedents or sources of information exist to guide resource managers, researchers, or the general public in understanding its short-and long-term impacts and initiating appropriate management responses.
The potential effects of MPB and other bark beetle outbreaks on forest structure and carbon cycling differ from those of other disturbances such as wildfires.MPB preferentially selects and attacks large-diameter individuals of its host tree species (lodgepole and other pines) and largely ignores smaller diameter individuals (Amman, 1977;Klutsch et al., 2009).Trees killed by MPB gradually lose their needles and branches over a period of 6 or more years (Klutsch et al., 2009).After attack, the proportion of biomass in standinglive, standing-dead, and downed-dead wood pools changes as needles, twigs, branches, and eventually the trees themselves fall, but for the most part total biomass (and carbon stored) in lodgepole pine forests may change little, as decomposition rates in these ecosystems are generally low (Son et al., 2010;Brown et al., 2004;Fahey, 1983).However, substantial reductions in carbon sequestration rates can occur following an MPB outbreak because of the loss of photosynthetic capacity following tree mortality (Kurz et al., 2008).The rate at which forest structure and carbon stocks and sequestration rates will return to pre-outbreak levels is not known, making it difficult to quantify the long-term impact of MPB on carbon cycling at regional and national scales.
The occurrence of wildfires, as well as insect outbreaks, has been increasing in recent years in the western US (Westerling et al., 2006).Wildfire is also a selective disturbance agent: tree mortality following fire can vary among size classes and species, depending on forest type, fire intensity, and numerous other factors.The effects of fire are generally more variable than those of an MPB epidemic, both within the perimeter of the disturbance and among disturbances.For instance, wildfires in lodgepole pine and subalpine coniferous forests typically occur during conditions favorable to high-severity crown fire initiation and spread, and few trees survive (Lotan and Perry, 1983).In contrast, wildfires in ponderosa pine (Pinus ponderosa) forests range from low to high severity, with variations in severity both within and among fires (Baker, 2009).When fire intensity is low to moderate, the thick bark of older trees effectively insulates them from fire damage, while the younger trees with thinner bark suffer mortality (Ryan and Reinhardt, 1988;Michaletz and Johnson, 2007).When fire intensity is high, a large amount of live and dead biomass may be consumed.Thus, the amount of standing-live and -dead biomass remaining after fires depends on forest composition and structure prior to the fire, as well as factors influencing fire behavior such as weather and topography.Similar to insect outbreaks, the proportion of biomass in the standing-live, standing-dead, and downeddead pools changes as needles, twigs, branches, and trees fall.However, unlike insect outbreaks, wildfires typically consume a proportion of live vegetation, dead surface fuels, and organic soil layers.The immediate result of fires is a release of greenhouse gases to the atmosphere as live and dead fuels are consumed (Seiler and Crutzen, 1980).Shortly after fires, biomass may be transferred among standing-dead and downed-dead biomass pools as dead trees fall.Fires also have long-term effects on forest processes as well as structure, altering rates of nutrient cycling, decomposition, photosynthesis, and regeneration such that years to decades may pass before carbon stocks return to pre-fire conditions (Dale et al., 2001;Swift, 2001;Turner et al., 1998;Kashian et al., 2013;Williams et al., 2012).
Disturbances drive successional trajectories by altering species composition, seedling establishment, and numerous abiotic and biotic factors in the post-disturbance environment (Turner et al., 1998;Collins et al., 2011), and different disturbance types clearly have varying effects on certain components of the ecosystem.Serotinous lodgepole pine regenerates abundantly in mineral soil that is exposed after stand-replacing fires (Lotan and Perry, 1983;Fahey and Knight, 1986).However, following MPB outbreaks, mineral soils may remain covered by a thick organic layer and litter that may hinder lodgepole pine seedling recruitment and regeneration (Collins et al., 2011).Therefore, forest vegetation recovery may be more dependent on existing seedlings and saplings (advance regeneration) and trees (Klutsch et al., 2009;Collins et al., 2011).Lodgepole pine is also relatively shade intolerant, and canopy cover remaining after an MPB outbreak may inhibit growth of existing lodgepole pine, while favoring more shade-tolerant species, such as Engelmann spruce (Picea engelmannii) and subalpine fir (Abies lasiocarpa) (Claveau et al., 2002;Collins et al., 2011).In addition to individual species' responses to the conditions created by either fire or MPB disturbance, direct or indirect competitive or facilitative interactions among species will influence the trajectories of succession (Callaway et al., 2002;Connell and Slatyer, 1977).Because carbon storage in forested ecosystems varies with forest age, composition, and management history (Bradford et al., 2008;Pan et al., 2011), these successional pathways and changes in community structure and composition should be considered in efforts to quantify the long-term impacts of different disturbances on carbon stocks and fluxes.
Several field studies in lodgepole-pine-dominated ecosystems have evaluated the effects of MPB (Collins et al., 2011;Diskin et al., 2011;Klutsch et al., 2009;Pfeifer et al., 2011;Kayes and Tinker, 2012;Pelz and Smith, 2012), wildfire (Anderson and Romme, 1991;Buma, 2011;Turner et al., 1997), or multiple disturbances (Buma and Wessman, 2011;Sibold et al., 2007) on forest structure and species composition.Some studies of MPB disturbance have also simulated changes in succession (Collins et al., 2011;Diskin, 2010) and carbon stocks (Pfeifer et al., 2011) over time.Numerous studies at various scales have evaluated the effects of fire on carbon cycles (Kashian et al., 2006(Kashian et al., , 2013;;Smithwick et al., 2009).In this paper, we characterized how species composition and carbon stocks in lodgepole pine forests might change in response to each disturbance: MPB and wildfire.Our approach built on the strengths of previous studies by comparing changes over time in multiple variables across two major disturbance types simultaneously.We used vegetation simulation modeling, initialized with field data collected after an MPB epidemic in Grand County, Colorado, to compare trajectories of forest growth and carbon stocks over 200 yr for three scenarios: (1) no disturbance, (2) MPB epidemic, and (3) wildfire.We assessed differences in species composition and aboveground carbon stocks among scenarios over 200 yr and asked (1) when do species composition and carbon stocks recover to pre-disturbance levels?and (2) when do species composition and carbon stocks converge among scenarios?

Study area
Our 1000 km 2 study area is located in eastern Grand County, Colorado (105 • 43 32 W to 106 • 0 47 W and 39 • 54 58 N to 40 • 18 2 N; Fig. 1).Forests in the study area are between approximately 2300 and 3400 m in elevation and are generally composed of even-aged stands of lodgepole pine, sometimes with a secondary stand structure composed of subalpine fir and Engelmann spruce seedlings and saplings.Quaking aspen (Populus tremuloides) is also common.A large percentage of the study area is public land with a multidecade history of wildfire suppression (Sibold et al., 2006).Until recently, the area's disturbance history consisted primarily of high-severity fires, occasional mixed-severity fires, and episodic bark beetle infestations at endemic population levels (Sibold et al., 2006).Beginning in 1996, an extensive and severe MPB outbreak started in the Southern Rocky Mountains and reached Grand County, Colorado, by 2002 (Tishmack et al., 2005).The peak of tree mortality in Grand County occurred between 2005 and 2008 (Klutsch et al., 2009), making our study area an ideal location for studying the impacts of MPB outbreaks on vegetation and carbon storage.

Field methods
In 2010, we collected field data to characterize forest composition and structure after the MPB outbreak.To ensure that the range of variability in biophysical gradients and overstory mortality in our study area was represented, we selected plot locations using a stratified random sampling approach.Strata included the following layers: years since peak overstory mortality (derived from the Forest Health and Monitoring Aerial Surveys; 1 yr, 2-3 yr, 4-5 yr, and 5+ yr); elevation (elevation quartiles between 2300 and 3400 m); and aspect (north, south, east, west, and flat), for a total of 80 different combinations of strata.We restricted plot locations to public lands within the study area that were greater than 90 m from roads and classified as Rocky Mountain Lodgepole Pine Forest in the LANDFIRE existing vegetation type layer (Rollins, 2009).We attempted to randomly position 2-3 plots in each combination of strata, but manually revised plot locations based on accessibility; plots located in potentially dangerous and inaccessible locations were moved within strata.We ultimately collected field data at a total of 119 plot locations.
At each plot, we measured trees, seedlings, saplings, and surface and canopy fuels using the Fire Effects Monitoring and Inventory Protocol: FIREMON (Lutes et al., 2006).Within fixed (8 m) radius plots, we measured tree diameter at breast height (DBH), total height (measured using a Haglöf Vertex Laser), status (live or dead), number of years dead, proximate cause of mortality (e.g., MPB; disease; unknown), and bearing and distance from plot center.We used www.biogeosciences.net/10/8203/2013/Biogeosciences, 10, 8203-8222, 2013 the criteria described by Klutsch et al. (2009) to estimate year of death for all trees killed by MPB based on the color and amount of needles and branches remaining.Trees were defined as having woody stems with a diameter at breast height (DBH) ≥ 12 cm and total height ≥ 1.4 m.Seedlings and saplings were counted within a 3.6 m-radius subplot and classified according to species and diameter (saplings) or height (seedlings) classes.Saplings were defined as any woody stems with a DBH < 12 cm and height ≥ 1.4 m, while seedlings had a DBH of < 12 cm and a height of < 1.4 m.We measured the depth of litter and duff, recorded downed woody debris by size class, visually estimated herbaceous ground cover, and measured overstory canopy cover using a spherical densitometer at standard locations along three Brown's transects between 5 and 25 m from plot center (Brown, 1974).We entered all data on trees, seedlings, saplings, live and dead fuel loads, and plot attributes such as slope and aspect into the Forest Vegetation Simulator (FVS) model, Central Rockies Variant .

Modeling methods
We conducted analyses with data from 97 plots that had an overstory dominated by lodgepole pine and an understory dominated by lodgepole pine, aspen, or subalpine fir.We excluded 22 plots that had an overstory dominated by subalpine fir or Engelmann spruce.We initialized FVS with data from the 97 lodgepole-dominated plots, and ran simulations of forest growth over a 200 yr time period to model three scenarios described in detail below.
The FVS, created and applied widely by the United States Department of Agriculture (USDA) Forest Service, is an aspatial growth and yield model that estimates forest growth based on tree recruitment, growth rates, and mortality, informed by algorithms for tree-(e.g., mean annual increment) and stand-level (e.g., basal area) production (Keyser and Dixon, 2008).FVS uses separate growth models for small and large trees, where the breakpoint DBH between small trees and large trees varies among species.The smalltree model predicts height growth first, and diameter is predicted from height growth, as height is a driver in tree regeneration when trees compete for space and resources.FVS predicts small-tree height growth based on equations using stand density, crown ratio and site characteristics (Keyser and Dixon, 2008).The large-tree model predicts diameter growth and changes in height using a suite of functions with predictors including species-level coefficients, DBH, tree age, stand-level basal area, and site index (Keyser and Dixon, 2008).Background mortality levels are predicted for individual trees based on DBH.Density-related mortality is predicted for the stand based on stand density index and the maximum stand density index.After stand-level densityrelated mortality is predicted, it is allocated to individual trees as a function of the contribution of each tree to stand-level basal area, species weights, and crown ratio (Keyser and Dixon, 2008).
Carbon stock estimates were made for 4 pools, representing total aboveground carbon, standing-live, standing-dead, and downed-dead biomass using the Fires and Fuels Extension (FFE) to FVS (Rebain, 2010;Reinhardt and Crookston, 2003).The standing-live carbon pool measured carbon in live trees, including stems, branches, and foliage, but not roots.The standing-dead pool included stems and branches and foliage, but not roots of dead trees.The downed-dead wood pool included all downed-dead wood, regardless of size.The total aboveground carbon pool included all of the above categories as well as carbon contained in the biomass of herbs, shrubs, roots, litter, and duff.FVS estimates carbon in these pools (except litter and duff) by multiplying biomass by 0.5 (Penman et al., 2003).Litter and duff in the downed-dead biomass pool is converted to carbon by multiplying estimates by 0.37 (Smith and Heath, 2002).Biomass of tree boles and crowns is estimated in FVS using national biomass equations (Jenkins et al., 2003) or the default FVS-FFE methods (Rebain, 2010).We selected the default FVS-FFE methods, which are considered to be more accurate because they use region-specific allometric equations from the National Volume Estimator Library (see Keyser and Dixon, 2008) to predict volume.Volume is converted to biomass with speciesspecific density values for boles (Brown et al., 1977) or using the equations in Brown and Johnston (1976) for crowns.During each cycle of an FVS simulation, some crown and bole material is transferred to fuel and litter pools.Transfers to fuel, litter and duff biomass pools in the model are based on ecological processes such as tree growth and mortality and snag dynamics (Rebain, 2010).FVS simulates decomposition of surface fuels and litter over time using a constant proportional loss model from four size-specific decay rates.Duff decay is estimated from a single decay rate where two percent of decayed matter from each fuel and litter pools are added to duff each cycle (Rebain, 2010).Snags are also modeled to decompose, decaying from hard dead wood to soft dead wood over time, or reach 64 % of the snags' original density (Rebain, 2010), which is a logarithmic relationship to decay rate for a particular size class.

Modeling scenarios
Three simulations were initialized in FVS: (1) a control scenario free of disturbance; (2) the actual MPB outbreak that affected the study area; and (3) a simulated high-severity wildfire scenario.For the control scenario, we reclassified the status of all MPB-killed lodgepole pine trees to live prior to running the simulation.This scenario was intended to represent baseline undisturbed conditions of vegetation composition and carbon over time, to provide a benchmark for comparison with possible changes after disturbances.For the MPB scenario, the plot data were unchanged and represented actual conditions in 2010, by which time 68 % of lodgepole pine basal area had been killed by MPB.
For the wildfire scenario, we parameterized the model using the same input data as was used in the control scenario, started the simulation in 2005, and scheduled a severe wildfire in 2006.Parameters set in the FFE to FVS (Reinhardt and Crookston, 2003;Rebain, 2010) for the wildfire included 29 • C temperature, very low fuel moistures, and 8 m s −1 wind speed; these represent extreme fire weather conditions.FVS fire parameters also include the proportion of the stand area that is burned.We set this value at 0.93, which was the maximum proportion of area classified with low, moderate, and high burn severity within fire perimeters in Colorado in the Monitoring Trends in Burn Severity (MTBS) database between 1984 and 2010 (Eidenshink et al., 2007).The results of the wildfire simulation were severe and mortality of understory and overstory trees was nearly complete.
In all three scenarios, we adjusted a number of FVS parameters to better reflect conditions in our study area.We disabled aspen resprouting, using the NoSprout FVS keyword, as it tends to result in unrealistic aspen densities in the Southern Rocky Mountains (personal communication Don Vandendriesche, 7 March 2013).We classified our plots into 3 types based on the dominant understory species -lodgepole pine, aspen, or subalpine fir (the overstory was dominated by lodgepole pine in all cases) -and then set values for maximum basal area (BAMax), stand density index (SDIMax), and species-specific mortality rates (TreeMax) for each plot type based on estimates derived from Forest Inventory and Analysis (FIA) data for the Southern Rocky Mountains (Tables 1 and 2; Don Vandendriesche, unpublished data; Vandendriesche, 2010b).
To generate species-specific regeneration rates, we used the REPUTE program in FVS (Vandendriesche, 2010a) with the Southern Rocky Mountains FIA data to calculate (impute) the numbers of new seedlings and saplings to add, based on overstory conditions, for each species in each time step in the simulations.In the MPB scenario, we assumed that the regeneration estimates generated by REPUTE were realistic for the simulation period because MPB-killed trees lose their needles and fall at relatively slow rates over time in this area and recent field studies have not identified a strong regeneration pulse in MPB-affected stands (Collins et al., 2011;Diskin et al., 2011;Kayes and Tinker, 2012;Pelz and Smith, 2012).However, Rocky Mountain lodgepole pine is known to have high and often variable regeneration rates after fire, depending on factors such as cone serotiny and burn severity (Turner et al., 1997;Kashian et al., 2004).Predicting post-fire seedling densities can be difficult, especially when a limited amount of data exist to train models (Robinson, 2008;Vandendriesche, 2010a).We assumed that post-MPB outbreak seedling and sapling densities measured in our plots in 2010 were representative of minimum post-fire seedling densities that might be possible given the biophysical setting of each plot.We also assumed that the species composition in the post-fire cohort of seedlings would more reflect pre-fire overstory species composition, plus a post-fire influx of seedlings from the serotinous cones of lodgepole pine.Therefore, we took a two-step approach to estimating the size and composition of the "post-fire" seedling cohort.First, we allocated the total counts of seedlings and saplings recorded on each plot in 2010 among species based on the proportion of stems of each tree species on the plot before the fire (i.e., representation of that species in the overstory).Second, we multiplied the numbers of lodgepole pine seedlings in this simulated post-fire cohort by a factor reflecting the probable proportion of lodgepole pine trees with serotinous cones in the study area, which would be expected to release additional seeds after fire.Past studies have indicated that 64 % (Aoki, 2010) and 65 % (Turner et al., 1997) of lodgepole pines are serotinous in Rocky Mountain National Park, Colorado (an area partially overlapping our study site), and Yellowstone National Park, Wyoming, respectively.Based on these results, we increased lodgepole pine seedling densities such that additional serotinous seedlings represented 64 % of the total lodgepole pine seedlings.This resulted in an average post-fire lodgepole pine seedling density of 8060 seedlings ha −1 in our first simulated time step (2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019).This density is lower than some previously published field measurements of lodgepole pine seedlings in the Yellowstone area 1 yr after a fire (14 000-91 000 seedlings/ha; Anderson and Romme, 1991) or 2 yr after a fire (50 000-140 000 seedlings/ha; Turner et al., 1997), and greater than other published estimates of conifer seedling densities 8 yr after a fire in subalpine forest (including lodgepole pine) in Colorado (317-1400 seedlings/ha; Buma, 2011;Buma and Wessman, 2011).We also allowed every aspen present in our plot to resprout after the simulated wildfire.Initial sensitivity analyses suggested that the FVS simulation results were relatively insensitive to post-fire seedling densities as long as enough post-fire seedlings were introduced to establish a new cohort of trees.

Statistical analyses
FVS simulation results for all three scenarios included attributes of individual trees, as well as stand-level characteristics such as trees per hectare (TPH) and basal area (BA; m 2 ha −1 ) per species, and carbon stocks in the various biomass pools.Simulation results were generated in 10 yr increments and are referred to hereafter by the first year of the decade (e.g., 2010 represents 2010-2019).Simulation results were not normally distributed; therefore, we used quantiles (e.g., median, 95th percentile, etc.) to describe the data and nonparametric statistical tests to determine if differences among scenarios were significant.To address our first question, we compared live tree density, basal area, and carbon stocks among scenarios for each decade (2010-2210) using Mann-Whitney U tests (Mann and Whitney, 1974)  pre-outbreak conditions, (i.e., values from the 2010 decade of the control scenario).Because 240 tests were performed (3 scenarios × 20 decades × 4 species), alpha levels were adjusted from 0.05 using a Bonferroni correction, such that differences between scenarios were assumed to be significant when p values were less than 0.0002 (0.05/240).We used the same approach to address our second question and determine if and when values for live tree density, basal area, and carbon stocks in the MPB and wildfire scenarios would converge with (i.e., have statistically insignificant differences from) those values in the undisturbed scenario.

Initial species composition and basal area
In the three scenarios that we simulated, four primary species were present in stands in 2010: lodgepole pine, Engelmann spruce, subalpine fir and quaking aspen (Table 3).A few other species were present in limited quantities, such as white fir (Abies concolor) and Douglas fir (Pseudotsuga menziesii; results not presented).The control scenario had the greatest initial density and basal area of live trees and the fewest standing-dead trees.Tree species composition was largely dominated by lodgepole pine (median density 800 TPH, median basal area 27.1 m 2 ha −1 ; Table 3).Aspen, Engelmann spruce, and subalpine fir were present in a very limited num-ber of plots and had low density and basal area relative to lodgepole pine (Table 3; Fig. 1; Fig. 2).Sapling and seedling composition was also dominated by lodgepole pine, with median densities of 500 TPH and 250 TPH, respectively.Aspen, subalpine fir, and Engelmann spruce saplings and seedlings were present on less than 1/2 the plots (median value of 0 TPH).On a few plots where these species were present (e.g., in the 95th and 100th percentiles), subalpine fir and aspen sapling and seedling densities were substantial (10-67 % of lodgepole pine sapling and seedling densities).Initial conditions for the MPB scenario were based entirely on the measured field data.Lodgepole pine was dominant and had a median tree density and basal area of 200 TPH and 4.5 m 2 ha −1 , respectively (Table 3).Lodgepole pine tree density and basal area had a wide range of variability, reaching as high as 1400 TPH and 34.1 m 2 ha −1 in spite of the extensive MPB mortality.The median tree density and basal area of dead lodgepole was 550 TPH and 21.2 m 2 ha −1 , respectively (Table 4).Lodgepole pine trees killed by MPB had significantly larger diameters (mean DBH = 22.1 cm) than live lodgepole pine trees (mean DBH = 16.9 cm; 2-sided t test p value < 0.05).Tree density and basal area of other tree species and advance regeneration (seedlings and saplings) in the MPB scenario were the same as in the control scenario because they were based on the same, unmodified field data.
Initial conditions for the wildfire scenario were created using the control scenario inputs and simulating a high-severity fire in 2006.Not surprisingly, the fire killed the majority of seedlings, saplings, and trees (Tables 3 and 4); fire-caused mortality (94 % of total BA) was greater than MPB-caused mortality (76 % of total BA).A very small proportion of trees survived the fire, mostly lodgepole pine, but also some subalpine fir and Engelmann spruce.The cohort of seedlings established after the fire was largely dominated by lodgepole pine, which was present in all plots and had highly variable densities in 2010-2019, with a median of 3220 TPH, minimum of 0, and maximum of 46 160 TPH (Table 3).Lesser amounts (up to around 590 to 3410 TPH) of subalpine fir and Engelmann spruce were present in 14 and 13 out of 97 plots, respectively (Table 3).Aspen was also present in only 13 out of 97 plots in 2010-2019, but had relatively high sapling densities (around 2700-18 500 TPH; Table 3) in the plots where it was present.

Initial carbon stocks
Total aboveground carbon was similar among the three scenarios in the first decade of the simulations (2010-2019); and median values were 8.4, 7.6, and 6.5 kg C m −2 for the control, MPB, and wildfire scenarios, respectively (Table 5; Fig. 4).Variability was high, and the maximum values for total aboveground carbon were greatest in the control scenario (21.3 kg C m −2 ), followed by the MPB (20.6 kg C m −2 ), and wildfire (16.8 kg C m −2 ) scenarios.Differences in total aboveground carbon in the 2010 decade were only significant between the control and wildfire scenario, based on Mann-Whitney U tests.
The differences among scenarios in the amount of carbon in the 3 separate pools we compared (standing-live, standingdead, and downed-dead) highlighted the different impacts of the MPB and wildfire disturbances (Table 5, and Fig. 4).The amount of carbon in the standing-live pool was greatest in the control scenario (median of 4.7 kg C m −2 ).Relative to the control, median standing-live carbon was 80 % less in the MPB scenario and 97 % less in the wildfire scenario; roughly half of the plots in the wildfire scenario had no standing-live carbon in the 2010 decade.The amount of standing-live carbon was significantly different for all possible comparisons of the various scenarios (e.g., control vs. MPB, control vs. wildfire, and MPB vs. wildfire).
Standing-dead carbon in the first decade of the simulations was similar in the MPB and wildfire scenarios (median of 2.9 and 2.5 kg C m −2 , respectively), and both scenarios had significantly more standing-dead carbon than the control scenario (Table 5 and Fig. 4).Downed-dead carbon was similar between the control and MPB scenarios (Table 5 and Fig. 4); largely because both scenarios started with the same initial values.However, by 2019, some MPB-killed trees had fallen, resulting in slightly larger values of downed-dead carbon in the MPB than control scenario (median of 1.0 vs. 0.8 kg C m −2 ).The wildfire scenario had the largest amount of carbon in the downed-dead pool in the initial decade of the simulations (median of 1.9 kg C m −2 ) and it was significantly larger than both the control and MPB scenarios (Table 5 and Fig. 4).

Simulated changes in species composition
Lodgepole pine remained the dominant tree species in all three scenarios throughout the simulations when measured in terms of both live tree density and basal area; however, the rates and magnitudes of change varied among scenarios and when comparing pre-disturbance to post-disturbance conditions (Figs. 2 and 3).In the control scenario, the median tree density of lodgepole pine gradually decreased and basal area steadily increased over most of the 200 yr simulation period (Figs.2a and 3a).In the MPB scenario, lodgepole pine tree density increased until 2110, after which it decreased; basal area showed a more continual increase over time at a more rapid rate.Lodgepole pine tree density was significantly lower in the MPB scenario than in the control scenario until 2060, but similar thereafter; basal area remained significantly lower until 2110 (Figs.2b and 3b).Relative to the pre-outbreak conditions (i.e., 2010 in the control simulation), lodgepole pine tree density and basal area returned to pre-outbreak levels in 2090.
In the wildfire scenario, there was an initial large increase in lodgepole pine density (seedlings and saplings), but it took 50 yr of growth before individuals reached the size (DBH ≥ 12 cm) at which they were counted as trees.Thus, the density of lodgepole pine trees remained extremely low (median of 63 TPH) through the 2060s, increased rapidly as the initial sapling cohort matured, reaching a median value of 1500 TPH by 2090, and then gradually declined to a median of 560 TPH by 2210 (Fig. 2c).Post-fire lodgepole pine tree densities were significantly lower than both the control and MPB scenarios until 2070 and significantly greater thereafter.Relative to pre-fire conditions, lodgepole pine tree densities in the wildfire scenario were significantly lower before 2070 and then significantly greater between 2090-2119 and 2180-2219 (Fig. 2c).Changes in lodgepole pine basal area in the wildfire scenario were similar to changes in tree density during the first 100 years of the simulation; basal area remained low during the first five decades and then increased rapidly from 2060 to 2100 (Fig. 3C).Basal area remained high in the last 100 yr of the simulation, unlike tree density, which ultimately declined.Decadal differences in basal area were significant between the control and wildfire scenario through the 2070s and between the fire and MPB scenario through the 2060s.Basal area returned to pre-fire levels by 2080 (Fig. 3c).
Live tree density and basal area of quaking aspen, Engelmann spruce, and subalpine fir did not differ as strongly among the three different scenarios or relative to preoutbreak levels as did lodgepole pine (Figs. 2 and 3).These three species were present in very low densities, if at all, in most plots, and present in large densities in a small number www.biogeosciences.net/10/8203/2013/Biogeosciences, 10, 8203-8222, 2013  of plots; median values for density and basal area were often at or near zero, but the 95th percentile and maximum values were often high (Figs. 2 and 3).The responses of these species to the 3 modeling scenarios were variable and visible in the results even if the differences were not statistically significant (Figs. 2 and 3).The density and basal area of quaking aspen trees were low initially, and then increased near the middle of all simulation periods as saplings matured, but with subtle differences among the scenarios in the timing and magnitude of the increase.Quaking aspen tree density and basal area peaked earliest in the wildfire and MPB scenarios (around 2060) and latest in the control scenario (around 2100).The maximum and final values for tree density and basal area of quaking aspen were greatest in the wildfire scenario (maximum 2900 TPH; final 930 TPH), followed by the MPB and then control scenarios (Figs. 2 and 3d-f).Engelmann spruce densities were relatively low throughout all scenarios (maximum value 2000 TPH; Fig. 2g-i).Tree density increased up until 2060 to 2100 depending on the scenario and then slowly declined (Fig. 2g-i).The maximum and 95th percentile of Engelmann spruce basal area increased over time under all the scenarios, but with the greatest magnitude of increase in the wildfire scenario, followed by the MPB and then control scenarios (Fig. 3g-i).Significant differences among the scenarios in Engelmann spruce density and basal area were found in the initial and final decades of the simulations, and between the fire and control scenarios between in the 2120s, 2130s, and 2140s.Neither density nor basal area of spruce was significantly different from predisturbance levels at any time.
Subalpine fir tree densities also increased in the initial half of the simulation period, and maximum density values reached relatively high levels in both the control and MPB scenarios (3020 and 3150 TPH, respectively; Fig. 2jk).In the MPB scenario, subalpine fir tree density was significantly greater than pre-outbreak levels in the 2070 and 2080 decades.After reaching maximum values, tree density steadily decreased over time in all scenarios.In contrast, basal area increased over time in all scenarios, reaching maximum values similar to those of lodgepole pine in both the control and MPB scenarios (Fig. 3j-l).However, in the wildfire scenario, subalpine fir tree density and basal area had markedly different magnitudes and rates of change compared to the other two scenarios.After the wildfire, subalpine fir tree density reached a maximum of only 1750 TPH in the 2070 decade, then showed a very slow decline, with significantly lower values than in the control scenario from 2090 to 2130.Maximum values for basal area increased at a steady rate to 53 m 2 ha −1 by 2210, although subalpine fir basal area was never significantly different from pre-fire levels (Fig. 3l).

Simulated changes in carbon stocks
Total aboveground carbon increased for all scenarios over time, but the range of variability over time was greatest under the wildfire scenario (Fig. 4a-c).Decadal differences between the control and MPB scenarios were significant starting in 2020 and ending by 2090, largely because total aboveground carbon was lower in the MPB scenario.From 2090 on, decadal differences between the control and MPB scenario were not significant.Total aboveground carbon in the MPB scenario was not significantly different than the preoutbreak levels until the 2060s, after which pre-outbreak levels were exceeded by the MPB scenario.Total carbon in the wildfire scenario, compared to pre-outbreak levels, was significantly lower in the 2010-2030s, was not significantly different in the 2040-2060s, and was significantly greater in the 2070s and thereafter (Fig. 4c).It took longer for levels of total aboveground carbon in the wildfire scenario to converge with the control scenario (2120) than it took for the MPB scenario to converge with the control scenario (2090; Fig. 4b  and c).
There was considerable variability among the 3 scenarios in standing-live carbon in the first decade of the simulations, but standing-live carbon increased and ultimately leveled off to similar values in all three scenarios (Fig. 4d-f).After MPB disturbance, standing-live carbon stocks were significantly smaller than pre-outbreak levels through the 2040s and significantly larger after 2080 (Fig. 4e).Decadal differences between the MPB and control scenarios were significant through the 2120s.In the wildfire scenario, standing-live carbon was significantly lower than pre-outbreak levels through the 2050s, increased at the most rapid rate of all scenarios and was significantly greater than pre-disturbance levels from 2090 and thereafter (Fig. 4f).Decadal differences between the wildfire and control scenario were not significant after 2150 (Fig. 4d and f).Differences between the wildfire and MPB scenario were only significant in the first 8 decades of the simulations (Fig. 4e and f).In both the wildfire and MPB scenarios, there was a decrease in the 5th percentile of standing-live carbon starting around 2120 and 2130 (Fig. 4e  and f), corresponding to decreases in quaking aspen densities and basal area (Fig. 2e and f; Fig. 3e and f).
Median levels of standing-dead carbon were low in all 3 scenarios.In the control scenario, standing-dead carbon increased initially in some plots, likely because FVS applied "background" mortality rates to a small percentage of the trees that had been killed by MPB but we recoded as live (Fig. 4g).Overall, standing-dead carbon in the control scenario increased until the 2150s before slightly decreasing over the remaining 60 yr.As expected, standing-dead carbon was initially high in the MPB scenario.It declined to a minimum in the 2040s as standing-dead carbon was transferred to the downed-dead carbon pool, increased from 2050 to 2200, and then declined slightly in the 2210 decade (Fig. 4h).Standing-dead carbon in the MPB scenario was significantly greater than pre-outbreak levels for all decades of the simulation.Standing-dead carbon in the wildfire scenario was also initially high, declined rapidly to a minimum in the 2030s, then increased slowly through the 2140s before decreasing again (Fig. 4i).As in the MPB scenario, standing-dead carbon in the wildfire scenario was significantly greater than pre-outbreak levels for all decades of the simulation.Differences between the wildfire and control scenario were not significant in the 2020s and after 2070 (Fig. 4g and i).Standingdead carbon in the wildfire and MPB scenario differed significantly during certain decades; values were greater in the MPB scenario from 2020 to 2039 but greater in the wildfire scenario from 2110 to 2169 (Fig. 4h and i).
Median levels of downed-dead carbon were intermediate between those of standing-live and standing-dead carbon in all simulations, and showed a general increase over time.In the control simulation, downed-dead carbon (Fig. 4j) was significantly lower than downed-dead carbon in the MPB scenario (Fig. 4k) and wildfire scenario (Fig. 4l) until the 2080s, after which differences among the scenarios were largely insignificant.Downed-dead carbon in the MPB scenario had an early increase from 2020 to 2039, remained relatively stable for the next 8 decades and then increased slightly over the remainder of the simulation period (Fig. 4k).Downed-dead carbon in the MPB scenario was significantly greater than pre-outbreak levels after the 2010 decade (Fig. 4k).In the wildfire scenario, downed-dead carbon was also high initially and increased in the 2020s, then decreased slightly, and then increased again in the 2070s and over the remainder of the simulation (Fig. 4l).Downed-dead carbon in the wildfire scenario was significantly greater than pre-outbreak levels throughout the entire simulation period (Fig. 4l).

Discussion
Our results identified both similarities and differences in the potential future response of lodgepole pine forest to the two major disturbances we examined, an actual MPB epidemic and a simulated wildfire.In answer to our first question (at what time might species composition and carbon stocks return to pre-disturbance levels?), our simulations predicted that lodgepole pine would remain dominant after both dis-turbances, with basal area returning to pre-disturbance levels after 70-80 yr.Other tree species (quaking aspen, Engelmann spruce, and subalpine fir) showed variable responses among scenarios but were not very abundant in our plots either before or after disturbance.Standing-live carbon also demonstrated a relatively rapid return to pre-disturbance levels after both MPB and wildfire (40 and 50 yr, respectively).In answer to our second question (at what point might species composition and carbon stocks in all scenarios converge?),we found that the basal area of lodgepole pine after MPB and wildfire was similar to that in the undisturbed scenario after 100 yr and 80 yr, respectively.The proportions of other tree species present rarely diverged from undisturbed conditions.Levels of standing-live carbon in the MPB and wildfire scenarios converged with those in the undisturbed control scenario after 120 and 140 yr, respectively.Standing-dead carbon varied among scenarios over time, but downed-dead carbon pools converged among all scenarios after 70 yr.In general, wildfire appeared to cause a greater initial magnitude of change in most of the variables we examined than did MPB disturbance, but the timeline of recovery to pre-disturbance and (or) undisturbed conditions was relatively similar after both disturbances.

Disturbance and vegetation composition
The MPB epidemic killed 65 % of lodgepole pine trees, equivalent to 76 % of basal area, at our study sites.Past stud-ies in the Southern Rocky Mountains have reported similar changes: Klutsch et al. (2009) reported a 71 % decrease in live lodgepole pine basal area in the Arapaho National Forest, CO; Collins et al. (2011) reported a 68 % decrease in live lodgepole pine basal area in Fraser Experimental Forest, CO; Diskin et al. (2011) reported a 64 % reduction in live basal area in Rocky Mountain National Park, CO; and Kayes and Tinker (2012) reported a 70 % decrease of live lodgepole pine basal area in the Medicine Bow National Forest, WY.The mortality we recorded may have been slightly higher than in these studies either because a longer period of time passed between the outbreak and data collection, or because the overstory in our plots contained greater proportions of the MPB's primary host, lodgepole pine.
Our simulations predicted that lodgepole pine will remain dominant after MPB disturbance over time, supporting the findings of several recent studies (Klutsch et al., 2009), especially in cases when severe MPB disturbances opened the canopy and when advanced regeneration was also dominated by lodgepole pine (Sibold et al., 2007;Amman, 1977).A recent study by Teste et al. (2011) suggested that lodgepole pine may be resilient to insect outbreaks as seeds held within cones may remain viable and can be released during a persistent seed rain lasting 9 or more years after MPB attack.Others have demonstrated that changes in overstory species composition are possible, but are largely determined by the dominant species in advanced regeneration in addition to non-host species present in the overstory (Collins et al., 2011;Diskin et al., 2011;Pelz and Smith, 2012).Sites with conditions that favor the growth of other species such as quaking aspen, subalpine fir, and Engelmann spruce are more likely to experience shifts in composition than sites that favor lodgepole pine (Collins et al., 2011;Diskin et al., 2011;Kayes and Tinker, 2012;Peet, 1981;Pelz and Smith, 2012;Romme and Knight, 1981).Our plots contained low levels of aspen (> 1 TPH present on only 50 % of plots in the absence of disturbance), but aspen did increase after MPB, occupying twice as many of the MPB-infested plots as control plots when at its peak.A similar increase in aspen after MPB was also reported or predicted by Collins et al. (2011), Diskin (2010), Diskin et al. (2011), Kayes and Tinker (2012), and Pelz and Smith (2012).In contrast, subalpine fir increased in the control scenario relative to the MPB scenario, probably because of its greater shade tolerance under the denser lodgepole pine canopy (Collins et al., 2011); however, these differences were only statistically significant for a few decades in the middle of the simulation period.Our results, and a growing body of literature, suggest that the impacts of MPB disturbances on forest structure may be relatively short-lived, that the majority of northern Colorado stands affected by MPB currently meet minimum stocking standards for density of advance regeneration (Collins et al., 2011;Kayes and Tinker, 2012), and that basal area will recover to pre-outbreak conditions within 80 yr or less (this study and Collins et al., 2011) in the absence of forest management or additional disturbances.
The simulated wildfire killed 94 % of all trees' basal area, representing a more severe disturbance to the forest overstory than the MPB epidemic.In addition, it essentially killed all seedlings and saplings (which were not impacted in the MPB scenario).However, the long-term results of the wildfire simulation reflected those of numerous field studies that have documented rapid recovery of Rocky Mountain lodgepole pine after fire (Lotan et al., 1985;Lotan and Perry, 1983;Lotan, 1976;Fahey and Knight, 1986;Romme and Knight, 1981;Peet, 1981;Turner et al., 1997).Our results predicted that lodgepole pine basal area would recover to pre-disturbance levels in 70 yr, a time frame similar to recovery in the post-MPB scenario (80 yr).The decadal outputs in FVS showed that the primary process driving this re-sponse was rapid growth of the large post-wildfire cohort of seedlings, dominated by lodgepole pine, that we established on the plots based on the findings of field studies in similar ecosystems (Anderson and Romme, 1991;Buma, 2011;Turner et al., 1997).Moreover, the growth of trees after the simulated fire was probably less affected by shading from residual canopy trees than in the MPB scenario, which may have contributed to the slightly more rapid recovery of basal area after the fire.Smaller cohorts of seedlings (also dominated by lodgepole pine) added by REPUTE at later time steps made a much more minor contribution to total basal area.Other species also represented only a small proportion of the overstory, although some had a stronger response to the fire scenario than to either the MPB scenario or the undisturbed conditions.In particular, aspen reached an earlier peak and maximum levels on post-fire plots, as found in other studies (see Lotan and Perry, 1983;Romme et al., 1995).Engelmann spruce also established at higher densities and reached greater basal areas during the fire simulation, a pattern consistent with the findings of some post-wildfire field surveys and dendrochronological analyses in subalpine forests (Johnson and Fryer, 1989;Knapp and Smith, 1982).In contrast, the shade-tolerant subalpine fir remained at much lower levels in the wildfire scenario compared to the other scenarios.However, these three species were present on less than 50 % of our plots before disturbance, and reached high densities on less than 5 % of plots after disturbance.

Disturbance and carbon stocks
In the absence of disturbance, carbon stored in standing-live, standing-dead and downed-dead biomass pools largely increased during the 200 yr time span of our control simulation.Compared to our disturbance-free scenario, the immediate impacts of disturbances on carbon storage were substantial, and resulted in 78 % and 97 % less carbon in the standing-live biomass pool after the MPB outbreak and simulated wildfire, respectively.In the wildfire scenario, some of that carbon was lost to the atmosphere through biomass combustion and the resulting emissions.In both the wildfire and MPB scenarios, the remaining biomass was transferred to the standingdead pool and subsequently to the downed dead pool in the first few decades of our simulations.The amount of carbon in the standing-live pool recovered to pre-disturbance conditions after 40 and 50 yr in the MPB and wildfire scenarios, respectively.However, standing-live carbon stocks did not converge with the control scenario until 120 and 140 yr following the MPB and wildfire disturbances, respectively.Both the MPB and wildfire disturbances resulted in an increase of biomass in the standing-dead and downed-dead pools; consequently, total carbon was not significantly different between the control and MPB scenarios initially (through the 2050s), and was only different for the first 3 decades in the wildfire scenario.This finding suggests that much of the live carbon lost by forests in these disturbances is retained within the system (e.g., in the MPB scenario) or quickly regained (e.g., in the wildfire scenario), although some time may pass before they reach the point where they would have been in the absence of disturbance.
In spite of the recognized importance of insect outbreaks on carbon cycling (Goward et al., 2008;Hicke et al., 2012a;Kasischke et al., 2013;Running, 2008;Liu et al., 2011), few of the previous studies examining the impacts of MPB outbreaks on forest vegetation structure and composition in Southern Rocky Mountain lodgepole pine forests have explicitly evaluated long-term changes in carbon stocks.Kurz et al. (2008) used the carbon budget model of the Canadian Forest Sector to simulate cumulative impacts of MPB outbreaks from 2000 to 2020 in British Columbia, Canada.Their methods account for carbon lost through decomposition, fire emissions, and harvest, and carbon gained through regrowth of the remaining, undisturbed portions of stands.Their results suggest that forests that were formerly atmospheric carbon sinks (17.2 g C m −2 yr −1 ) will switch to serving as sources (−42.4g C m −2 yr) and remain sources past 2020.Pfeifer et al. (2011) simulated changes in carbon following an MPB outbreak in a northern Idaho forest.They found that, on average, carbon returned to pre-outbreak levels within 25 yr, but their plots experienced lower mortality than in our region (33 % of lodgepole pine killed) and contained higher proportions of non-host tree species.They did identify variation in recovery rates among plots, which they attributed primarily to differences in size and growth rates of the surviving trees.Edburg et al. (2011) used the Community Land Model version 4 to examine sensitivities of carbon and nitrogen to outbreak severity, outbreak duration and tree fall rates.They predicted that 80 to 100 yr would pass before carbon stocks in vegetation recovered to pre-outbreak conditions.Our estimates of a 40 yr recovery time are within the range identified by these previous studies, and taken together they indicate that recovery of forest carbon storage to preoutbreak conditions is not only possible, but is likely to occur over a relatively short time span.Presumably, the growth of new regeneration as well as existing seedlings, saplings, and surviving overstory trees permits rapid exploitation of the increased resources available after an MPB outbreak (Hicke et al., 2012a;Pfeifer et al., 2011;Rhoades et al., 2013).
The wildfire we simulated caused a greater initial magnitude of change and a different pattern of recovery in carbon stocks compared to the MPB infestation, but the overall effects on standing live carbon lasted only slightly longer (recovery to pre-disturbance conditions after 50 vs.40 yr; convergence with undisturbed conditions after 140 vs. 120 yr).Throughout the wildfire simulation, carbon stored in the standing-live pool tended to be less than in the MPB scenario, but differences in the amount of carbon stored in the standing-dead and downed-dead carbon pools were mostly minimal between the two scenarios.Our results showed that after advance regeneration was removed by the simulated wildfire, a lag of several decades occurred until new seedlings became established, demonstrated rapid growth, and created a lodgepole overstory with density and basal area greater than that in the undisturbed forest by 80 yr postwildfire.At the end of the 200 yr simulation, the standinglive carbon pool still showed a steady increase, in contrast to the trends in the MPB and undisturbed scenarios, which exhibited stagnation and decline, respectively.
In-depth studies of patterns of carbon recovery after wildfires in lodgepole pine forests have been completed in Yellowstone National Park, WY.Kashian et al. (2006) examined patterns of ecosystem productivity following fires and made projections to estimate that the total carbon lost from the 1988 fires would be recovered within 230 yr.However, a more recent study suggested that 80 % and 90 % of total carbon would recover within 50 and 100 yr, respectively, but large variability was observed among stands (Kashian et al., 2013).They concluded that forests are resilient to disturbance unless there is vegetation type conversion following disturbances.Our results differ from Kashian et al. (2006Kashian et al. ( , 2013) ) in that we project a more rapid recovery of carbon in the standing-live biomass pool; however, our study did not explicitly incorporate carbon losses through decomposition of downed-dead woody biomass, and tree growth rates may be greater in our study area than in Yellowstone National Park.

Uncertainties and limitations
Several limitations and sources of uncertainty in our analyses could be improved upon in future studies; these include assumptions about, and the methods to simulate, tree regeneration and tree fall rates through time, components of the carbon cycle not simulated by FVS, and the potential for additional future disturbances to occur over time and affect vegetation composition and carbon cycling.
Although regeneration models have been developed within FVS for some regions of the US, this has not been done for the Central Rockies Variant, presenting a significant limitation to users who are interested in long-term projections of vegetation change in this area.Users have 3 options: determining and scheduling the appropriate number of seedlings to add to simulations at the appropriate time steps themselves; using Climate-FVS (Crookston et al., 2010); or using the REPUTE post-processor in FVS to impute future regeneration from their own data or additional sources of data (Vandendriesche, 2010a).Scheduling seedling establishment requires making assumptions about the species and density of seedlings likely to establish at given times in the future.Climate-FVS introduces functions to initiate forest regeneration based on stand conditions and climate predictors, but implies that users understand probable climate-related vegetation dynamics in their study system and are explicitly incorporating climate-change scenarios into their analyses.For the third option, REPUTE imputes future seedling and sapling counts based on forest type and stand structure categories.be acquired from the nearest regional FIA data set if the user's own data are not sufficient.However, regardless of the data source, use of REPUTE assumes that the forest structure conditions under which future seedlings will establish are represented in the input data set.If users wish to model significant future changes to forest conditions, this assumption may not be valid.
We used REPUTE for this study and found that it generated reasonable predictions of seedling establishment in the post-MPB scenario, given our and others' observations that forest structure changes (such as rates of needle fall, tree fall, and new seedling establishment) after MPB-caused mortality in our study area have been relatively gradual in recent years (Collins et al., 2011;Kayes and Tinker, 2012;Pelz and Smith, 2012).Moreover, because our field data used to parameterize REPUTE were drawn from locations selected using a randomly stratified sampling design, we are confident the full range of variability of conditions affecting forest regeneration after an MPB outbreak were included.However, because few wildfires have occurred and been surveyed recently in unmanaged lodgepole-dominated forests in our study area, we were less confident that either REPUTE or the regeneration data published after fires in the Yellowstone area represented realistic post-fire regeneration patterns in Colorado.Additionally, since REPUTE is not stochastic, there is no variability in the number of seedlings and saplings introduced for a given combination of forest type and stand structure conditions.Thus, REPUTE may not fully represent the range of variability observed in regeneration studies (Kashian et al., 2005;Sibold et al., 2007;Turner et al., 1998).FVS and REPUTE may be further limited for longterm simulation of forest vegetation studies in that species' occurrence in future time steps is restricted to predefined plot-types; dispersal dynamics are not included.Finally, in this and other aspatial models, it is not possible either within or among plots to account for spatial patterns of forest structure and how they influence seed dispersal and regeneration rates.Modeling potential future regeneration dynamics after significant, complex patterns of disturbance in complex ecosystems remains a considerable challenge.
Key pools and fluxes in the carbon cycle are not addressed by FVS, which primarily quantifies aboveground carbon pools.However, belowground carbon pools are important and can be substantial (Kashian et al., 2013).Nearly 20-30 % of the biomass of lodgepole pine is belowground in the root system, and allocation of primary productivity to root systems can be high (Comeau and Kimmins, 1989;Litton et al., 2004Litton et al., , 2007;;Jackson et al., 1996).Soils can also contain a significant amount of carbon; past studies estimated that organic and mineral carbon in soil represents 25 to 45 % and 20 to 48 % of the total carbon in lodgepole pine forests in SW Colorado (Kueppers and Harte, 2005) and NW Wyoming (Litton et al., 2004), respectively.Biogeochemical changes to soil carbon resulting from insect disturbances occur on a slower temporal scale than do changes in the car-bon in the aboveground pool, requiring different modeling and measurement techniques for accurate assessment (Edburg et al., 2012;Rhoades et al., 2008).Our research was aimed toward quantifying changes in the recovery of aboveground carbon pools, and could only account for a subset of all carbon pools and fluxes.However, carbon allocation patterns between above-and belowground pools vary little with stand age and we would expect carbon stocks in belowground pools to track those in aboveground pools (Litton et al., 2004).
Not only are carbon storage and sequestration affected by stand density, composition and age, but they also depend on rates of litter decomposition and tree fall (Kashian et al., 2004(Kashian et al., , 2006;;Edburg et al., 2011).Our analyses also did not address the impacts of disturbances on respiration by soil microbes and how altered levels of downed-dead wood or decomposition rates would influence their activity as this was beyond the scope of our study.However, past studies have demonstrated that respiration rates following disturbances are important and can shift a forest ecosystem from a carbon sink to a carbon source (Edburg et al., 2011;Goetz et al., 2012;Kashian et al., 2006;Kasischke et al., 2013;Kurz et al., 2008).
Our analysis simulated changes in forest vegetation composition, structure, and carbon storage after two primary disturbances, but did not incorporate the effects of future disturbances, such as additional insect outbreaks (e.g., spruce or pine beetles), parasites (e.g., mistletoe), blow down, management activities, or wildfires.Potential interactions between past and future disturbances (Sibold et al., 2007) can impact conifer recruitment in a species-specific manner (Buma and Wessman, 2012), and may be required before a shift in species composition occurs (Amman, 1977).Improvements to our modeling methodology could include belowground carbon pools, carbon fluxes, a more dynamic (and, ideally, field-validated) approach to modeling regeneration and tree fall rates, and a focus on the potential effects of additional future disturbances over time.Using models like Climate-FVS (Crookston et al., 2010), the effects of potential climate change on species regeneration, growth rates, and mortality could be incorporated, in addition to climate-driven changes in disturbance frequency.These improvements were beyond the scope of this analysis, but could be incorporated into future analyses and studies.
Strengths of our approach include the simultaneous consideration of changes in both carbon stocks and vegetation composition after two major disturbances in lodgepole-pinedominated forests, and a more complete effort to model regeneration after these disturbances than has previously been published (Smithwick et al., 2009;Diskin, 2010;Pfeifer et al., 2011;but see Collins et al., 2011).Given the lack of feasibility of conducting replicated, long-term "natural experiments" to compare forests' response to MPB infestation and wildfire, and the confounding factors and complexities encountered by retrospective studies that have examined multiple past occurrences of these two disturbances (Axelson et al., 2009;Dordel et al., 2008;Sibold et al., 2007), our approach represents a useful evaluation of the resilience of lodgepole-pine-dominated forest in the Southern Rockies to each of these major disturbance processes, and our results suggest that recovery is relatively rapid after both.

Summary
We simulated potential changes in forest vegetation and carbon pools for a scenario without disturbance and two scenarios with disturbances: MPB and wildfire.Our results showed an immediate impact of disturbance on forest vegetation structure and carbon storage, but minimal changes in species composition were found and lodgepole pine remained the dominant tree species under all three scenarios.Losses of carbon and changes among carbon pools in lodgepole pine forests affected by MPB and wildfire disturbances were relatively short-lived under scenarios that did not incorporate further disturbances.In our MPB scenario, standing-live carbon rebounded 40 yr after disturbance.However, carbon pools were impacted more severely and for a slightly longer time period after a simulated wildfire; standing-live carbon rebounded to pre-fire levels within approximately 50 yr after disturbance.Substantial differences in standing-dead and downed-dead carbon persisted throughout 200 yr of simulation.These results support those of other recent field-based and disturbance modeling studies and emphasize that lodgepole pine forests are largely resilient to disturbances.
Disturbances are projected to increase in frequency and severity with climate change (Dale et al., 2001;Bentz et al., 2010;Raffa et al., 2008;Westerling et al., 2011;Littell et al., 2010), and the resilience of carbon pools in lodgepole pine forests to disturbances may be threatened (Smithwick et al., 2009;Kashian et al., 2006).If this occurs, the relationships between carbon storage, species composition, and trajectories of succession over time will become increasingly important to resource management as efforts to mitigate anthropogenic greenhouse gas emissions face greater uncertainties.

Fig. 1 .
Fig. 1.State boundaries and county boundaries for Colorado.Grand County, our study area, is highlighted in black.

Fig. 2 .
Fig. 2. Boxplots of live tree density by decade for the 3 simulation scenarios and 4 primary species.Species codes are PICO = lodgepole pine (Pinus contorta); ABLA = subalpine fir (Abies lasiocarpa); PIEN = Engelmann spruce (Picea engelmannii); and POTR = quaking aspen (Populus tremuloides).Whiskers span 100th to 95th and 0th to 5th percentiles; boxes span 95th to 5th percentile; and horizontal lines show median (50th percentile).Letters above whiskers indicate significant differences in values for a decade using a Mann-Whitney U test between the (a) control and MPB scenarios, (b) control and fire scenarios, and (c) MPB and fire scenarios.Asterisks indicate significant differences between control 2010 values for the MPB and fire scenarios.

Fig. 3 .
Fig. 3. Boxplots of live tree basal area by decade for the 3 simulation scenarios and 4 primary species.Species codes are PICO = lodgepole pine (Pinus contorta), ABLA = subalpine fir (Abies lasiocarpa), PIEN = Engelmann spruce (Picea engelmannii), and POTR = quaking aspen (Populus tremuloides).Whiskers span 100th to 95th and 0th to 5th percentiles; boxes span 95th to 5th percentile; and horizontal lines show median (50th percentile).Letters above whiskers indicate significant differences in values for a decade using a Mann-Whitney U test between the (a) control and MPB scenarios, (b) control and fire scenarios, and (c) MPB and fire scenarios.Asterisks indicate significant differences between control 2010 values for the MPB and fire scenarios.

Fig. 4 .
Fig. 4. Boxplots of (a) total aboveground carbon, (b) standing-live carbon, (c) standing-dead carbon, and (d) downed-dead carbon from 2010 to 2210 for the three scenarios assessed in this paper.Whiskers span 100th to 95th and 0th to 5th percentiles; boxes span 95th to 5th percentile; and horizontal lines show median (50th percentile).Letters above whiskers indicate significant differences in values for a decade using a Mann-Whitney U test between the (a) control and MPB scenarios, (b) control and fire scenarios, and (c) MPB and fire scenarios.Asterisks indicate significant differences between control 2010 values for the MPB and fire scenarios.
It requires field data to parameterize, which can sometimes www

Table 1 .
Adjusted values used for the FVS maximum basal area (BAMax) and stand density index (SDIMax) keywords for different plot types.

Table 2 .
Adjusted values used for the FVS species-specific mortality rate (TreeMax) keywords for different plot types and species.

Caldwell et al.: Simulated impacts of mountain pine beetle and wildfire disturbances 8211Table 4 .
Summary statistics of dead tree density and basal area by species for each scenario in the first time step (2010-2019) of the simulations.Species codes are PICO = lodgepole pine (Pinus contorta); ABLA = subalpine fir (Abies lasiocarpa); PIEN = Engelmann spruce (Picea engelmannii); and POTR = quaking aspen (Populus tremuloides).

Table 5 .
Summary statistics for carbon (kg C m −2 ) in various carbon pools for each scenario in the first time step (2010-2019) of the simulations.