Modeling the interactive effects of spruce beetle infestation and climate on subalpine vegetation

In the subalpine zone of the Rocky Mountains, climate change is predicted to result in an increase in the frequency and severity of spruce beetle outbreaks. Climate change itself may affect vegetation, potentially leading to changes in species composition. The direct and indirect effects of climate and disturbances on forest composition, biomass, and dynamics open the possibility for non-linear ecosystem responses. Modeling studies allow for the study of the interaction of these effects and their impact on the forest system. University of Virginia Forest Model Enhanced (UVAFME), an individual-based gap model that simulates forest dynamics and characteristics, is updated with a spruce beetle subroutine that calculates the probability for beetle infestation and potential mortality of each tree on a plot. The updated model is then run with multiple scenarios that combine beetle infestation with current or altered climate at sites across the southern Rocky Mountains. Results show that spruce beetle infestations acted to facilitate competition with invading lower-elevation species, resulting in an increase in the biomass of historically lowerelevation species and a further decline in Engelmann spruce biomass than occurred with solely bark beetle disturbance or solely climate change. We also found an initial enhancing effect between spruce beetle infestation and climate change; however, by the end of 100 yr of climate change and potential beetle infestation, climate had a dampening effect on spruce beetle infestation, through loss of host trees. These results are an important step in understanding the possible futures for vegetation of the Rocky Mountains as well as for spruce forests across the western United States and Canada.


INTRODUCTION
Disturbances such as fire, windthrow, and insect outbreaks are principal drivers of vegetation dynamics within high-elevation and high-latitude spruce forests and interact to affect forest composition and dynamics as well as ecosystem processes and biogeochemical cycling (Veblen et al. 1991a, Edburg et al. 2012, Goetz et al. 2012, Hansen 2013, Frank et al. 2014, O'Halloran et al. 2014. Outbreaks of the spruce beetle (Dendroctonus rufipennis Kirby), which infests spruce (Picea spp.) species and whose native range extends from the southern Rocky Mountains in Arizona to boreal Alaska and Canada (Wood 1982), have increased in recent years (O'Connor et al. 2015, USFS 2016, leading to widespread mortality and carbon losses throughout the western United States, Alaska, and Canada (Berg et al. 2006, Bentz et al. 2010). Many factors have been attributed to these recent outbreaks, including the availability of vast, contiguous areas of large-diameter spruce (DeRose et al. 2013, higher incidents of drought (Hebertson and Jenkins 2008, DeRose and Long 2012, Hart et al. 2014b, and increases in ambient temperatures (Bentz et al. 2010, Sherriff et al. 2011, DeRose et al. 2013.
Spruce beetle populations typically exist at relatively low, endemic levels, punctuated by occurrences of high, epidemic levels due to climate-, disturbance-, or forest structure-related triggers (DeRose et al. 2013). During endemic intervals, spruce beetles colonize downed spruce logs and may attack older or larger trees, though the success of these attacks is mediated by a combination of the health of the tree, its ability to defend itself, and the number of attacking beetles Frye 1977, Raffa et al. 2008). Trees with high vigor may fend off infestations by exuding resin, which trap and kill their attackers and their brood. A mass attack of many beetles, however, overwhelms trees' efforts, leading to successful infestations and subsequent tree mortality (Raffa et al. 2008). Weather and environmental/forest stand factors that increase susceptibility to spruce beetle attack (i.e., a high availability of weakened trees, high density and proportion of spruce, or droughts) can allow for more successful mass attacks that can build into widespread outbreaks (Schmid and Frye 1976, Berg et al. 2006, DeRose and Long 2012. In addition to forest and environmental characteristics, spruce beetle growth response to ambient temperatures also influences beetle population growth and thus outbreak success (Hansen et al. 2011), and this effect has been especially important in recent years (Sherriff et al. 2011, Hart et al. 2014a. Spruce beetles have a variable life cycle depending on temperature. Under lower to normal temperature conditions, spruce beetles complete their life cycle in two years, whereas anomalously warm temperatures allow completion of the life cycle in a single year (Werner andHolsten 1985, Hansen et al. 2001). This flexible voltinism results in faster growing populations under univoltine (one-year) life cycles compared to semivoltine (two-year) life cycles. With increasing ambient temperatures from climatic warming, more spruce beetles should switch from a semivoltine to a univoltine life cycle (Bentz et al. 2010). This switch in life cycles will likely lead to more beetles emerging and reproducing every year rather than every other year, potentially allowing for enhanced population growth relative to semivoltine populations (Hansen et al. 2001). This faster population growth may allow for more frequent mass attacks on spruce trees, and potentially more widespread outbreaks.
The combination of climate change's potential indirect effects on the bark beetle disturbance regime and its direct impacts on the performance of individual trees introduces the possibility for complex, non-linear ecosystem responses . Vegetation productivity and insect population growth both respond to short-term changes in weather, such as seasonal droughts, as well as long-term changes in climate, such as El Niño events or more directional climate change (Veblen 2000, Sherriff et al. 2011, O'Connor et al. 2015. Forest productivity and bark beetle population growth are accelerated by increasing summer and winter temperatures (Werner and Holsten 1985, Veblen et al. 1991a, Hansen et al. 2001, whereas droughts or periods of elevated hydrological demand can cause moisture stress in trees and compromise their defenses through loss of carbohydrate reserves (Fettig et al. 2013). Due to this confluence of increases in temperatures, drought, and evaporative demand, and the high availability of susceptible spruce trees, it is predicted that the frequency and severity of spruce beetle outbreaks will likely increase further with climate change (Bentz et al. 2010), and this appears to already be occurring (O'Connor et al. 2015). However, as bark beetles require adequate host material to reproduce, climate change may result in a decrease in insect outbreak simply by the reduction of suitable host trees as the range of spruce changes in response to altered climate (DeRose et al. 2013). The frequency and severity of wildfires within spruce systems are also predicted to increase with climate change (Westerling et al. 2006, Kashischke et al. 2010, potentially leading to increased fire-related mortality and a shift toward more fire-adapted species. The future of these ecosystems is thus becoming progressively unclear as climate change and disturbances act in concert to alter their structure, composition, and internal dynamics (Fettig et al. 2013). Mountains have experienced large-scale spruce  beetle outbreaks over the past decade, and infestations continue to expand across the region  (USFS 2016). Within these forests, fire, windthrow, and bark beetle outbreaks have the capacity to modify and interact with the structural and biological diversity of forest stands (Veblen et al. 1991a, Kulakowski and Veblen 2002, Sibold et al. 2007). This differential effect on tree size and species allows for unique disturbance interactions, through alteration of the resilience and/or resistance of post-disturbance forest stands to subsequent disturbances (Buma 2015). For example, high-intensity wildfire can increase stand resilience to windthrow by increasing the prevalence of smaller, more wind-firm stems which colonize the gaps created by stand-initiating fires (Kulakowski and Veblen 2002). Windthrow and lowintensity fires can decrease resilience to bark beetle outbreaks through tree damage and increases in downed trees that foster spruce beetle population growth (Schmid and Frye 1977, Geiszler et al. 1984, Christiansen et al. 1987, Rasmussen et al. 1996, Hood and Bentz 2007, Fettig et al. 2008, Mezei et al. 2014. In contrast, high-intensity fires decrease the likelihood of bark beetle outbreak by decreasing the availability of suitable host material (Veblen et al. 1994, Bebi et al. 2003, Kulakowski and Veblen 2006. Bark beetles may also interact with wildfire, though the exact nature of this relationship is still actively debated. Studies on the effect of mountain pine beetle (Dendroctonus ponderosae Hopkins) outbreaks in pine forests found increasing probability for active crown fires in the early stages of outbreak through increases in dry, flammable fuels , Jenkins et al. 2012, followed by decreasing likelihood of crown fire for as many as several decades following the loss of fine needles (Jenkins et al. 2014). A recent study by Andrus et al. (2016) found that spruce beetle outbreak had no effect on subsequent wildfire severity under extreme fire weather conditions but had a positive effect on subsequent fire severity under moderate fire weather conditions. These linked and compound interactions among climate, vegetation, and disturbances at multiple spatial and temporal scales may produce unexpected system responses and trajectories under future climate scenarios (Gustafson 2013, Buma 2015.

Many subalpine forests of the southern Rocky
Detailed vegetation modeling can be used to simulate these interactions and to assess the response of spruce beetles, and subsequently vegetation, to various climate change and disturbance scenarios (Seidl et al. 2007, Temperli et al. 2013, Buma 2015, Thom et al. 2017. Through high-resolution modeling, the direct and indirect impacts of climate change can be separated and combined , and potential disturbance-vegetation-climate interactions and their cascading effects can be analyzed (Buma 2015). In the case of interacting climate and spruce beetle infestations, plausible outcomes include an enhancing effect of climate on spruce beetle infestations through increased population growth and higher tree stress (Raffa et al. 2008), leading to greater spruce mortality than would be expected from simply the summation of climate and beetle-related mortalities, as well as a dampening effect of climate on infestations due to declining spruce hosts . Some combination of these interactions may also occur, and it is likely to change over time and with stand characteristics.
Individual tree-based models, or gap models, explicitly consider interactions of individual trees and individual stands with climate, vegetation, and various disturbances. The hierarchical nature of gap models allows for integration of tree-and stand-level responses to study emergent properties of forested landscapes Woodward 2011, Gustafson 2013). As such, they are valuable tools for answering such questions about the future of forested ecosystems and are capable of capturing the interactive dynamics between various vegetation drivers. In addition, gap models have been tested on their ability to hind-cast forest responses for past climate changes, increasing the confidence in predictions made for novel future climates. For example, Solomon et al. (1980), Solomon and Webb (1985), and Overpeck et al. (1990) used a gap model to reconstruct the pollen/climate data for eastern North American late Quaternary landscapes (Bartlein et al. 1986). They utilized palynological data and climate series built from these data to test gap model predictions of vegetation state based on this inferred climate. Bonan and Hayden (1990) used a paleo-climate derived from ocean sediment data as an independent proxy for climate to simulate the vegetation of Virginia ❖ www.esajournals.org during the last glaciation. These studies which tested the ability of gap models to predict vegetation response to past climate changes are a logical preamble to using these models to predict future climate responses (Shugart et al. 2018).
The University of Virginia Forest Model Enhanced (UVAFME) is an individual tree-based gap model, which simulates individual tree growth, mortality, and regeneration in response to climate, competition, and disturbances. The University of Virginia Forest Model Enhanced is an updated version of the gap model FAREAST (Yan and Shugart 2005) and has been successfully updated and tested for use in the southern Rocky Mountains (Foster et al. 2015. Previous studies with UVAFME in the Rocky Mountains predicted upward shifts in species zonation with increasing temperatures ); however, increasing disturbances were not considered. Climate change effects on wildfire as well as spruce beetle infestations are likely to further modify the species composition, structure, and biomass of the Rocky Mountains subalpine zone, potentially leading to novel species assemblages. Other forest models have been utilized within the Rocky Mountains region, ranging in scale from landscape to tree cohorts (Bugmann 2001, Schumacher et al. 2006, Keane et al. 2011, Temperli et al. 2015. This study builds on these applications through the simulation of individual tree growth and mortality, as well as vegetation response to and interaction with nutrient availability.
In the present study, a spruce beetle submodel is developed and implemented within UVAFME. Model simulations are conducted at sites within the southern Rocky Mountains subalpine zone with different combinations of spruce beetle presence and climate change to evaluate the relative and combined effects of beetle disturbance and changing climate on vegetation. These results are intended to advance the understanding of the possible futures for the study region, as well as for the greater range of the spruce beetle.

Model description and updates
UVAFME is an extension of the individualbased gap model FAREAST (Yan and Shugart 2005). The University of Virginia Forest Model Enhanced simulates the regeneration, growth, and death of individual trees on independent patches of a forested landscape, with each patch/ plot experiencing the same site and climate conditions. As a traditional gap model, UVAFME is only spatially distributed in the vertical dimension. Plots are assumed to have no direct horizontal spatial interactions with one another. A sample system for an actual forest area conforming to this assumptive framework would consist of a set of several hundred independent survey plots (in the present case, each 500 m 2 ) with species and size tallies of each tree on each plot measured annually. Through a Monte Carlo-style aggregation, these plots are averaged to characterize the forested landscape under specific site and climate conditions. Except for its more frequent update interval, output from the model resembles data from a system of replicate independent forest inventory plots Woodward 2011, Shugart et al. 2018). Optimal annual diameter increment growth for a tree of a given size and species is calculated using growth equations from the FAREAST model. Actual diameter increment growth for each tree is then calculated by adjusting this optimal annual diameter growth according to climate and site conditions on the site in conjunction with species-and tree size-specific tolerances (e.g., drought, temperature, and nutrient tolerance). This diameter increment growth is then used to update each individual tree's diameter at breast height (dbh) as well as other allometric tree characteristics (e.g., height, biomass, leaf area index) for each simulated year. Trees may die from stress if they have low diameter growth for several years, or from age-related factors. Trees may also die due to fire or windthrow disturbance. New trees can regenerate each year based on climate and environmental conditions, disturbance occurrence, species-specific tolerances, and species' seedbanks (which are also modified each year based on site and climate conditions). A complete description of the equations and subroutines within UVAFME can be found in Shuman et al. (2015) and in the UVAFME User's Manual in Foster et al. (2017). This model has been updated and parameterized for the southern Rocky Mountains (Foster et al. 2015, and it was found that model output on size structure and speciesspecific biomass agreed with inventory data within the subalpine zone. The model was also found to accurately predict successional dynamics ❖ www.esajournals.org following bare-ground initiation as well as the change in species composition with elevation present within the southern Rockies (Daubenmire 1943, Veblen 1986, Aplet et al. 1988. Following this demonstration of fidelity in the southern Rockies (Foster et al. 2015, UVAFME was chosen for this study to simulate forest response to shifting climate and disturbance regimes within the study region. The University of Virginia Forest Model Enhanced's windthrow and fire disturbance submodels were updated for this study to better represent the complex interactions between climate, disturbances, and forest stand structure. Windthrow in the model is stochastic and is based on a site-specific return interval. Previously, when windthrow occurred on a simulated plot, it would immediately kill all trees on the plot, regardless of size. Studies have shown that windthrow differentially affects trees of varying sizes, with larger trees having a higher probability of windthrow mortality than smaller, more wind resistant trees (Veblen et al. 1991b, 2001, Kulakowski and Veblen 2002, Rich et al. 2007). The model was updated to reflect these dynamics using equations from Rich et al. (2007): where p wind is the probability of tree mortality by windthrow, and dbh is the diameter at breast height (cm) of the simulated tree. This updated windthrow submodel allows for a more realistic simulation of the effect of windthrow on forest structure within the study region and also allows for better interaction between fire, windthrow, and bark beetles. Previous work with UVAFME in the study region involved updating and testing of the fire submodels to better reflect vegetation-fire interactions within the region . Fire in this model is stochastic, with a site-specific intensity and return interval (FRI). Each year of the simulation, each plot has a site-specific probability (i.e., 1/FRI) of a fire occurring, independent of other plots. Once a fire occurs, that fire's intensity value is randomly selected from a distribution based on the site-specific average fire intensity. Fire mortality is dependent on this intensity, as well as tree size and species-specific bark thickness parameters (Appendix S1: Table S1), with smaller and thinner-barked trees having a higher probability of mortality for a given fire intensity . Regeneration after fire is also modified by species-specific seedbank and seedling responses to fire, with some species experiencing a reduction and others an increase of their seedbank in response to fire occurrence. The addition of this fire submodel was found to produce forest characteristics and dynamics comparable to forest inventory data and descriptions of successional processes within the southern Rockies ; however, the fire module used in Foster et al. (2017) was not responsive to changes in temperature or evaporative demand. The fire submodel was updated for this study such that the fire probability of a site increases with an increasing site aridity index, defined as the ratio of precipitation to potential evapotranspiration (AI = P/PET) as in Feng and Fu (2013). For the first 100 yr of simulation, the fire probability is held at historical (i.e., input) values. A base aridity index (AI base ) for each site is calculated during this time (based on input distributions of historical monthly temperature and precipitation), and this base index is compared against each subsequent year's index. If the aridity index in any subsequent year is lower (i.e., drier) than the site's base index, the fire probability for that site that year is modified using the percentage difference between the base index and that year's index: Using this modification, the probability of fire occurrence can increase along with increasing evaporative demand, either due to lower precipitation or higher temperatures (and thus higher PET). This interaction between fire and climate has been widely predicted for various regions, including western North America (Jolly et al. 2015). With the addition of the effect of increasing aridity on fire probability in the model, the combined effect of changing climate, increasing fire, and potentially increasing insect infestations can be evaluated.

Spruce beetle submodel
A spruce beetle submodel was developed to calculate the probability of spruce beetle ❖ www.esajournals.org infestation of each host tree on each plot. This probability is based on three factor types, each operating at different scales: climate factors, which affect infestation probability at the site and plot level; plot factors, which affect infestation probability at the plot level; and tree factors, which affect each tree individually.
Climate factors include the probability of any one plot having spruce beetles with a univoltine (one-year) life cycle, based on spruce beetle phenology studies by Hansen et al. (2001Hansen et al. ( , 2011. Plot-level factors are based on individual plot characteristics, such as the average dbh of live spruce above 25.4 cm dbh, basal area, and percent spruce in the canopy, based on spruce beetle stand ratings from Schmid and Frye (1976). Other plot-level factors include those relating to recent windthrow events and the amount of downed, still green trees on the stand, based on studies implicating windthrow and downed logs as important precursors to bark beetle outbreaks (Wichmann and Ravn 2001). Finally, tree-level factors were developed based on tree size, stress level, and damage from recent fires. Once a tree becomes infested, it ceases growth (Frank et al. 2014), loses its needles after two years (Schmid and Frye 1977), and is marked as dead and added to the decomposition pool after five years. During the time when a tree is infested and still on the plot, it increases the infestation probability of directly and diagonally adjacent spruce trees on the same plot, based on studies finding proximity to other infested spruce trees as an important factor in determining infestation probability (Hart et al. 2014b). A full description of the equations and variables used in this spruce beetle submodel can be found in Appendix S1.
Previous work with UVAFME at these four sites  showed good agreement between model output and inventory data and descriptions of expected species composition and successional dynamics (Daubenmire 1978, Veblen 1986, Aplet et al. 1988). As in Foster et al. (2015Foster et al. ( , 2017, the model was parameterized for these sites using regional species information from Burns and Honkala (1990; see Appendix S1: Table S1 for species parameter values), and site, soil, and climate information from the US Forest Service and the National Climatic Data Center (Menne et al. 2012a, b). To incorporate the potential effects of upslope migration under future climate scenarios, range maps were used to determine which of the eleven major southern Rocky Mountain species (i.e., A. lasiocarpa, Juniperus scopulorum, P. engelmannii, Picea pungens, P. contorta, Pinus edulis, P. flexilis, Pinus ponderosa, Populus angustifolia, P. tremuloides, and P. menziesii) were eligible for colonization and growth at and in the vicinity around each site (Little 1971). Thus, even if a species was not present in the description for a site, it was eligible for colonization if it was shown as being present at that location in range maps by Little (1971). However, input species characteristics (Appendix S1: Table S1) and internal mechanics of the model determine whether those species are successful under the given climate, disturbance, and site conditions. Previous work with UVAFME in the southern Rocky Mountains  found that the model successfully produced the expected change in species composition with elevation, even when all major species were allowed to grow at all elevations.
To determine the response of subalpine vegetation to climate change, spruce beetles, and their interaction, several model simulations were conducted involving four different beetle/climate change scenarios at each site: (1) a control run with current climate and no spruce beetle disturbance; (2) current climate with potential recurring spruce beetle disturbance; (3) climate change without beetle disturbance; and (4) climate change with concurrent beetle disturbance. For all simulations, UVAFME was run with 200 independent, 500-m 2 plots. For the current climate simulations (scenarios 1 and 2), the model was run from bare ground for 800 yr. For climate  2 and 4), potential for beetle infestation was initiated at year 400, and run until year 800. All scenarios included fire disturbance with site-specific fire frequency and intensity, with climate-fire interactions (Eq. 2). Climate in UVAFME is driven by input state variables (average monthly temperature minima and maxima and precipitation, and their standard deviations). Each year, new monthly values are generated from the input distributions and converted into daily weather values for use in the calculation of hourly temperature for Appendix S1: Eq. S1 in the new spruce beetle submodel and in the daily soil moisture and nutrient subroutines (see UVAFME Users' Manual in Foster et al. 2017 for further description of these methods). Climate change input, in the form of changing monthly minimum and maximum temperatures and monthly precipitation, was derived from output from NCAR's Community Earth System Model (CESM) for the A1B and A2 IPCC climate change scenarios. These prescribed temperature changes resulted in a~3°C increase in temperatures for the A1B scenario, a~5.5°C increase in temperatures for the A2 scenario, and relatively little changes in precipitation for either IPCC scenario (Table 2).
NCAR's Community Earth System Model is considered suitable for the western United States due to its simulation of seasonal variation in temperature and precipitation as well as multi-year variability in Pacific sea-surface temperature (Westerling et al. 2011). Climate variability (in the form of input standard deviations of monthly temperature minima and maxima and precipitation) was held at historical values taken from original input climate data (Menne et al. 2012a, b). The 200-plot averages of size structure and speciesspecific biomass output were then compared across all model simulations at all four sites.

RESULTS
The response of subalpine vegetation to spruce beetles under current climate conditions was fairly similar across the four sites with the exception of FEF, where the spruce biomass loss was relatively lower (Fig. 2). Across all four sites, the addition of spruce beetle infestation under current climate scenarios resulted in about a 60-80% loss of spruce biomass at year 800 (at the end of the simulations) relative to year 800 biomass without beetles or climate change. There was also an increase in the biomass of non-host species (i.e., Abies lasiocarpa, Pinus contorta) between the control and beetle simulations, with site-specific differences. (see GLEES output in Fig. 3 and see Appendix S1: Figs. S3-S5 for other sites). At GLEES and Niwot Ridge, both subalpine fir and lodgepole pine increased in biomass, with fairly high proportional increases in lodgepole pine biomass at Niwot Ridge (~+98%; Appendix S1: Fig. S5). At Wolf Creek, quaking aspen and Douglas-fir also increased in biomass following spruce beetle outbreak (Appendix S1: Fig. S3). Additionally, there was a difference in the size structure of Engelmann spruce between the control simulation and solely beetle simulation (Appendix S1: Fig. S6). At GLEES, for example, the solely beetle disturbance simulation resulted in no large spruce stems (above 40 cm dbh) at year 800 (Appendix S1: Fig. S6b); it also had a much higher proportion of small spruce stems (below 10 cm dbh) than did the control simulation (Appendix S1: Fig. S6a). The simulated spruce beetle infestations thus resulted in a shift toward smaller trees and an increase in subdominant spruce trees (Appendix S1: Fig. S6) in addition to an increase in non-host species ( Fig. 3; Appendix S1: Figs. S3-S5). The addition of spruce beetles at GLEES also resulted in a fluctuation over time in spruce biomass (Fig. 4b), which corresponded to a fluctuation in spruce beetle-killed biomass over time aggregated across all plots (Fig. 5). A 10-yr running average shows that landscape-level beetle-killed biomass had a periodicity of about 100-150 yr. Variability within this average indicates smaller outbreaks occurring on individual plots within the larger landscape cycle (Appendix S1: Fig. S7). At the scale of an individual plot, infestations occur with higher frequency and at different times within each plot as spruce trees grow and become vulnerable to beetle attack. When beetlekilled biomass is aggregated across all 200 plots to the landscape scale, a longer-scale periodicity arises as an emergent property of the landscape.
Across all four sites, increasing temperatures without the addition of spruce beetle infestation resulted in a decline in Engelmann spruce biomass, an increase in more drought-tolerant subalpine species (i.e., P. contorta), and the introduction of lower-elevation species (i.e., Pseudotsuga menziesii, Pinus ponderosa; Figs. 2, 3; Appendix S1: Figs. S3-S5). In some cases, climate change was so detrimental as to completely or very nearly eradicate all subalpine species (i.e., Picea engelmannii, A. lasiocarpa), even without the addition of beetle infestations ( Fig. 2; Appendix S1: Fig. S5). Except for at Niwot Ridge, there was little difference in the biomass and species composition changes for the solely climate change simulations between Fig. 2. Biomass (tonnes C/ha) at year 800 of Engelmann spruce (Picea engelmanni) at all four subalpine sites for the control simulation (i.e., no beetle disturbance, current climate; "No BB, No CC"), the solely beetle disturbance simulation (i.e., beetle disturbance, current climate; "BB, No CC"), the solely climate change simulation (i.e., no beetle disturbance, climate change; "No BB, CC"), and the combination of beetle disturbance and climate change ("BB, CC"), along with 95% confidence intervals. Climate change simulations for graph (a) use the A1B climate change scenario; simulations for graph (b) use the A2 climate change scenario. those using A1B and those using the A2 IPCC scenarios (Figs. 2, 3; Appendix S1: Figs. S3-S5), even though the A2 scenario had larger overall temperature increases ( Table 2).
The effect of the combination of spruce beetles and climate change on Engelmann spruce biomass differed between the four subalpine sites. For GLEES, Wolf Creek Pass, and Niwot Ridge, the combination of spruce beetles and climate change resulted in lower Engelmann spruce biomass than did either factor alone (Fig. 2). At GLEES, beetle disturbance and climate change resulted in a further decrease in Engelmann spruce as well as the introduction of new lowerelevation species, ponderosa pine (P. ponderosa; Figs. 3, 4d). This shift in species dominance was also manifested in the changes in stand structure between simulations with and without beetle infestation. Without spruce beetle infestation, there were still many large-diameter spruce at year 800, even after climate change effects (Appendix S1: Fig. S8a). With climate change and spruce beetles, however, there were virtually no large-diameter spruce and a higher number of moderately sized trees of other species (Appendix S1: Fig. S8b). These results are similar for Wolf Creek Pass and Niwot Ridge (Fig 2; Appendix S1: Figs. S3, S5), however, at FEF climate change produced such a large loss of Engelmann spruce that the addition of spruce beetle infestation had little to no effect on spruce biomass (Fig. 2).
The relative effects of spruce beetles and climate change varied over time and among the different sites. Graphs of spruce biomass difference between the control simulation and all other simulations show the loss of spruce due to each factor individually and in concert (Fig. 6). The additive effect, which is simply the loss from beetle disturbance alone (scenario 2) plus the loss from climate change alone (scenario 3), represents what spruce biomass loss would have been if there had been no interaction at all between climate and beetle infestation. By comparing the Fig. 3. Biomass (tonnes C/ha) at year 800 of nonhost species (i.e., all but Picea engelmannii) at USDA Forest Service's Glacier Lakes Ecosystem Experimental Site for the control simulation ("No BB, No CC"), the solely beetle disturbance simulation ("BB, No CC"), the solely climate change simulation ("No BB, CC"), and the combination of beetle disturbance and climate change ("BB, CC"). Climate change simulations for graph (a) use the A1B climate change scenario; simulations for graph (b) use the A2 climate change scenario. (Fig. 3. Continued) ❖ www.esajournals.org additive and combined (i.e., loss from scenario 4) effects, we can evaluate whether there is an enhancing or dampening effect between beetle disturbance and climate change. When the additive effect is less negative than the combination effect, the spruce loss under the combination scenario is higher than would be expected and there is an enhancing effect between beetle disturbance and climate change. When the additive effect is more negative than the combination effect, spruce loss in the combination scenario is lower than would be expected, and there is a dampening effect between the two factors. At GLEES, beetle disturbance and climate seemed to enhance each other following the initialization of climate change (Fig. 6), resulting in greater biomass losses. However, toward the end of the simulation, around years 650-700, the additive loss was much greater than the combination loss, indicating that, as the effects of climate change progressed through time, there may have been an eventual dampening effect on spruce beetle outbreak. Beetle-killed biomass also declined over time under the combination scenario at GLEES and was less periodic than it was without climate change effects (Fig. 5). These results are similar at Wolf Creek (Appendix S1: Fig. S9); however, at Niwot Ridge, climate change had a mostly dampening effect on beetle disturbance (Appendix S1: Fig. S11), and at FEF climate change was so detrimental to Engelmann spruce biomass that mortality from spruce beetles became almost irrelevant by the end of the simulations ( Fig. 2; Appendix S1: Fig. S10).

DISCUSSION
Through the use of high-resolution modeling of individual tree interactions, we were able to evaluate the response of the Rocky Mountains subalpine system through the comparison of a combination of different scenarios. As direct and indirect interactions among disturbances such as wildfire and insect outbreaks with vegetation and climate have the potential to lead to unexpected and non-linear ecosystem responses (Buma 2015), it was important to investigate the drivers alone, as well as in combination. By using long time frames and a factorial approach, we were able to study these interactions as emergent phenomena of cross-scale forest and disturbance dynamics (Shugart et al. 2018). The biomass dynamics simulated in the solely beetle disturbance simulations are comparable to what has been found in various field studies (Figs. 2-4). A resampling study following spruce beetle infestation by Derderian et al. (2016) found that a recent spruce beetle infestation in a forest stand in northern Colorado resulted in a decline in Engelmann spruce biomass of about 70%, a slight increase in subalpine fir biomass (~2%), and a strong increase in the stem production of understory spruce (~50%) and fir (~80%). Another study investigated the effects of the 1940s spruce beetle outbreak across sites in central western and northwestern Colorado and also found that the growth rates of both subcanopy spruce and fir increased for several decades following infestation (Veblen et al. 1991a). As in these field studies, we found a strong decline in spruce biomass (Fig. 2), an increase in subalpine fir (Figs. 3, 4b; Appendix S1: Figs. S3-S5), and an increase in subcanopy stem density (Appendix S1: Fig. S6) as a result of spruce beetle infestation. This similarity between model output and independent inventory data shows the model's ability to simulate the impacts of spruce beetle outbreak on forest stand structure, an important effect that is feasible through the model's individual-based framework.
The periodicity seen in the biomass dynamics as a result of spruce beetle disturbance (Figs. 4b,5) is comparable to periodicities in spruce beetle outbreaks found by recent field studies, with outbreaks occurring every 50 yr in Alaska (Berg et al. 2006), and about every 100 yr in British Columbia and northwestern CO (Veblen et al. 1994, Zhang et al. 1999, and about every 40-60 yr in the Colorado Front Range (Hart et al. 2014a). In these simulations and in the field, once spruce beetles have killed most of the large-diameter spruce on a plot, the only remaining spruce trees are too small to be available for infestation, thus beetle-killed biomass decreases. These subcanopy trees grow, and eventually become large enough to be susceptible to beetle infestation, thus beetle-killed biomass increases (O'Connor et al. 2015). In this way, cycles of increasing and decreasing infestations arise over time, and are visible in landscape-level (i.e., 200plot averages) model output (Figs. 4b, 5), and emerge in field data as cycles of beetle outbreaks across the forest landscape. Within individual Fig. 6. Spruce biomass difference (tonnes C/ha) over time at USDA Forest Service's Glacier Lakes Ecosystem Experimental Site between the control simulation and the solely beetle disturbance simulation (red), the solely climate change simulation (green), and the combination of climate change and beetle disturbance (blue) for the (a) A1B and (b) A2 climate change scenarios. The additive (purple) line is the loss from beetle disturbance alone plus the loss from climate change alone. Shading represents 95% confidence intervals. plots/stands, smaller outbreaks occur with higher frequency (Appendix S1: Fig. S7); however, when averaged across all simulated plots (Fig. 5), the emergent cycle is longer. The timing of these infestation cycles may be important to consider when predicting the effects of droughts, El Niño events, or longer-scale climate change effects on spruce mortality from beetle outbreaks. Moisture stress has been shown to lead to increased tree susceptibility to beetle infestation, often leading to large-scale outbreaks when in conjunction with high availability of large-diameter spruce (Sherriff et al. 2011). Thus, if a severe drought or other outbreak-inducing event were to occur during a period of high host biomass, the availability of large-diameter spruce may allow for even more rapid growth of beetle populations following drought, potentially triggering largescale outbreaks across whole landscapes. In contrast, if a drought were to occur during a period of low host biomass, there may not be enough large-diameter spruce to sustain high populations of spruce beetles, even with the addition of drought stress and lowered tree defenses.
It is clear from the solely climate change simulations that local scale factors such as site characteristics, species composition, and climate play an important role in the response of subalpine vegetation to climate change. Even though the effect of spruce beetles on Engelmann spruce biomass was fairly consistent across all four sites, the effect of climate change was quite variable (Fig. 2). Increasing temperatures had only a moderate effect on spruce and fir biomass at GLEES and Wolf Creek (Figs. 2-3; Appendix S1: S3) and was not as detrimental to Engelmann spruce by itself as was beetle infestation (Fig. 2). In contrast, the A1B climate change scenario acted to strongly decrease spruce and fir biomass at FEF, and the A2 scenario resulted in strong losses in subalpine vegetation at both FEF and Niwot Ridge (Fig. 2; Appendix S1: Figs. S4, S5). Fraser Experimental Forest and Niwot Ridge are slightly drier sites (~70 and~50 cm annual precipitation, respectively) compared to GLEES and Wolf Creek (~100 and~110 cm, respectively), and this difference in overall climate may be driving the increased effect of elevated temperatures at Niwot Ridge and FEF. These results indicate that fine-scale patterns in climate, weather, and disturbance regimes should be considered when predicting the future state of vegetation within the Rocky Mountains.
In the combination scenarios, it seems that beetle infestation not only resulted in loss of spruce biomass, but that it may have also facilitated competition between Engelmann spruce and lower-elevation species (Appendix S1: Fig. S8). The reduction in large-diameter Engelmann spruce by spruce beetles and increase in larger diameter non-host trees modified the canopy dynamics in the combination scenario, therefore allowing successful colonization and growth of lower-elevation, non-host species such as ponderosa pine and Douglas-fir (Fig 3; Appendix S1: Figs. S3-S5). Thus, spruce beetle infestation in conjunction with climate change may help to open up the canopy to new, historically lowerelevation species. These low elevation species, which tend to have faster growth rates than the cold-adapted, shade-tolerant subalpine species (Burns and Honkala 1990), may then be able to outcompete and crowd out younger spruce trees from ever dominating the stand again.
This increase in lower-elevation species at the expense of subalpine species has been documented in other modeling studies (Rehfeldt et al. 2006, Crookston et al. 2010, Notaro et al. 2012, Jiang et al. 2013, Bell et al. 2014, Temperli et al. 2015. Notaro et al. (2012) utilized a dynamic global vegetation model (DGVM) to predict declines in suitable habitat for Engelmann spruce and subalpine fir with changing climate, and Rehfeldt et al. (2006) utilized a climate envelope approach to predict a decline in subalpine and alpine forests and an increase in lower-elevation forests and grasslands under climate change scenarios. A study utilizing a DGVM coupled with a global climate model also found a decrease in needle-leaved evergreen trees and an increase in shrubs and woodland with increasing temperatures (Jiang et al. 2013). The simulation results agree with these past studies and expand on their predictions by providing additional details on potential changes in size structure.
Other studies have also shown competition and species interactions to be a key factor in predicting species compositional change (Araujo andLuoto 2007, Zhang et al. 2015). A recent longitudinal study of over 27,000 trees found that competition accounted for the most variability in growth and mortality (Clark et al. 2011). While inhabitants of ❖ www.esajournals.org the lower elevations may be susceptible to higher moisture stress from increasing temperatures Breshears 1998, Breshears et al. 2005), those of the characteristically mesic subalpine zone may be more vulnerable to the negative effects of competition arising from those very same low elevation species escaping drought. These studies, along with the simulations presented here, suggest a grim future for Engelmann spruce and other hosts of the spruce beetle. It may be possible for spruce ranges to expand upward or northward (Hanberry andHansen 2015, Bretfeld et al. 2016); however, this possibility depends on adequate soil, available space for migration, and future spruce beetle outbreak dynamics (Raffa et al. 2008, Bell et al. 2014. The decline in beetle-killed biomass over time in the combination scenario (Fig. 5) and the eventual dampening effect between beetle disturbance and climate change (Fig. 6) found in these simulations are similar to findings by DeRose et al. (2013). They utilized forest inventory data and spruce beetle population metrics to predict future spruce beetle presence across the central and southern Rocky Mountains. Their results showed that while climate data were important, they were far outstripped by the importance of habitat variables like stand basal area and percent Engelmann spruce. As spruce beetles require adequate host material to survive and reproduce, it follows that even with the ability for accelerated population growth under increasing temperatures, flexible voltinism becomes a moot point in the absence of spruce biomass. A survey of disturbance/climate change literature by Seidl et al. (2017) also found a high number of insect-related studies which cited dampening effects between outbreaks and climate. Thus, from these results it seems that climate change and spruce beetles may initially enhance one another, through facilitation of competition and through increases in beetle population growth. However, eventually, due to the high infestationrelated mortality of large-diameter spruce and the effective suppression of small-diameter spruce by species like ponderosa pine and Douglas-fir, spruce biomass may decline so much that spruce beetles will become less and less important in their current ranges.
A recent study by Temperli et al. (2015) utilized the landscape model LandClim to investigate the response of subalpine vegetation in northern CO to increasing spruce beetle infestations and climate change. They similarly predicted a reduction in Engelmann spruce and an increase in Douglas-fir and ponderosa pine; however, their climate change scenarios were somewhat more extreme than the ones utilized here (+4.4°C and À9% precipitation; +5.2°C and +12% precipitation; +7.0°C and À29% precipitation for the Land-Clim study; see Table 2 for temperature changes employed in this study). They also predicted a decline in beetle-related mortality with climate change, and an eventual dampening effect from interactions between climate and spruce beetles at high elevations and found a more intense dampening effect between beetle disturbance and climate at low elevations.
The spruce beetle infestation model presented here and that of the Temperli et al. (2015) use similar spruce beetle susceptibility metrics, with notable differences in the modeling structure and unit of simulation for each model. Because Land-Clim uses a tree cohort-based forest modeling approach, trees that originated in the same decade are necessarily assumed to be the same size (Temperli et al. 2015). LandClim simulates tree mortality through reducing the stem number count of each age cohort via a mortality rate, whereas UVAFME directly simulates the growth and potential mortality of each individual simulated tree. By simulating each individual tree on several hundred independent plots, UVAFME allows for more fine-scale interactions between trees, their environment, and disturbances to be borne out on each individual plot, without having to average across tree ages or assume similar vegetation structure across plots. This is especially important for beetle and fire disturbance, as mortality due to these is strongly size dependent (Keane et al. 2011, Jenkins et al. 2014. As a result, individual tree interactions arising from differences in tree sizes, the facilitation of competition between Engelmann spruce and lowerelevation species, and the effects and interactions among non-stand-replacing disturbances (i.e., fire, windthrow, and beetles) are captured at the level of the individual tree within UVAFME. The differences in the length of outbreak cycles on an individual plot (Appendix S1: Fig. S7) and that averaged across 200 plots (Fig. 5) highlight UVAFME's ability to simulate fine-scale forest interactions with insect infestations at the standlevel, as well as those arising as emergent properties of larger landscapes. Though the spruce beetle model described in the present study and that in LandClim were developed and parameterized independently, the results are complementary (Temperli et al. 2015), increasing the confidence in the predictions of both models. The results of both studies suggest that interacting disturbances and their impacts on vegetation dynamics are important for understanding the trajectory of change in this region (Buma 2015. With an individual-based model, we see that these disturbances are variable across the different sites and scenarios, as a result of stand structure, composition, successional dynamics, and tree-level interactions. When considering the potential increase in non-host and lower-elevation species, it is important to note that lodgepole pine (Pinus contorta) and ponderosa pine (P. ponderosa) are principal hosts of the mountain pine beetle (Dendroctonus ponderosae), a close relative of the spruce beetle that has recently caused extensive tree mortality throughout the western United States (Logan andPowell 2001, Powell andBentz 2009). Douglas-fir and subalpine fir also have associated bark beetles, the Douglas-fir beetle (Dendroctonus pseudotsugae), and western balsam bark beetle (Dryocoetes confusus), respectively (McMilin et al. 2003, Negron andPopp 2009). Infestation by mountain pine beetles, Douglas-fir beetles, and western balsam bark beetles was not included in these simulations, and thus, projections of increased ponderosa pine, lodgepole pine, Douglas-fir, and subalpine fir biomass must be taken with these mortality agents in mind.
In these simulations, spruce beetle infestations do not influence tree susceptibility to fire probability and severity. Studies in pine forests have shown that forest stands recently attacked by mountain pine beetle (i.e., within one to four years of infestation) have a higher probability for fire ignition and spread due to the decreased moisture level in the leaves and increased levels of flammable defense chemicals of attacked trees (Jenkins et al. 2012). Later on in the stages of an outbreak, however, the potential for active crown fires has been shown to decline due to decreased surface-to-canopy fuel continuity . These interactions also seem to vary with location and climate. A recent study by Andrus et al. (2016) found that under moderate levels of fire weather spruce beetle outbreak had an enhancing effect on subsequent fire severity. A study by Hansen et al. (2016) in the Kenai Peninsula of Alaska found that spruce beetle outbreaks increased fire probability in northern Kenai, an area dominated by both black (highly flammable, beetle resistant) and white (low flammability, beetle mortality prone) spruce, but decreased fire probability in white-spruce-dominated southern Kenai. Future work with UVAFME which explicitly models forest fuels and wildland fire dynamics will seek to simulate these additional bark beetle-wildfire interactions. These interactions may be key in predicting forest species composition under a warming, drying, and potentially more flammable landscape (Jolly et al. 2015).

CONCLUSIONS
The forest dynamics simulated in this study and the climate-vegetation-disturbance interactions that shaped them can only be attained with a high-resolution model such as UVAFME. Other ecological models have studied the effects of climate and disturbances on vegetation but have often been limited by the lack of explicit consideration of individual species and individual trees. A comparison of 40 terrestrial biosphere models (which simulate ecosystem processes and the distribution of vegetation, generally at the level of plant functional types) found high uncertainty and variability in both the magnitude and sign of annual carbon flux over the Alaskan Arctic, a region strongly driven by climate and disturbances (Fisher et al. 2014). The results presented here demonstrate the importance of tree-level interactions between vegetation, climate, and disturbances. These simulations have shown that although the amount of infested spruce decreases over time under a warming climate due to declining host availability, there may be an initial enhancing effect between climate change and spruce beetle outbreak, as infestations facilitate competition between Engelmann spruce and invading lower-elevation species, further driving composition and structural change in subalpine systems. Spruce forests are important for their contribution to the global carbon budget, for their influence on slope stability, for ❖ www.esajournals.org commercial timber and water resources, and for fostering tourism in the region (Burns andHonkala 1990, Fettig et al. 2013). They are a source of great natural beauty and are home to many important wildlife species (Bentz et al. 2010). The loss of these valued ecosystems would thus carry with it environmental, economic, and social implications.
The disturbance submodels developed in this work, along with additional harvest and management submodels, can be used to inform and improve management techniques aimed at alleviating some of the predicted mortality and enhancing the recovery from such infestations (Jenkins et al. 2014). The infestation submodel described here can also be applied to the study of spruce beetles within boreal Alaska and Canada as well as other bark beetle species, such as the mountain pine beetle, the Douglas-fir beetle, the western balsam bark beetle, or Ips spp. bark beetles. Outbreaks of these additional bark beetles have similar effects on forest stand structure and interact with other disturbances and climate (Bentz et al. 2010). Thus, individual-based modeling will be important with these systems for effective prediction of how species and forest size structure might respond to changing conditions. As shifting climate and disturbance regimes continue to alter forest dynamics and interactions between vegetation and vegetation drivers, individual-based modeling will become a valuable tool to investigate the possible futures of spruce forests.

ACKNOWLEDGMENTS
This work was funded by the University of Virginia, the NASA Postdoctoral Program, and grants from the VA Space Grant Consortium Graduate Fellowship (VSGC FY 15-16 to A.C.F., project title "Understanding spruce beetle outbreak dynamics and their response to climate change through remote sensing and ecological modeling") and the National Fish and Wildlife Foundation (grant number: 0106.12.032847) to A.C.F. We thank Katherine Holcomb of Advanced Research Computing Services at the University of Virginia for providing computational resources and technical support that have contributed to the results reported within this paper. We also thank two anonymous reviewers for their valuable critiques and suggestions which contributed to the final version of this manuscript.