Modeled interactions of mountain pine beetle and wildland fire under future climate and management scenarios for three western US landscapes

Background: Mountain pine beetle (MPB) is a native disturbance agent across most pine forests in the western US. Climate changes will directly and indirectly impact frequencies and severities of MPB outbreaks, which can then alter fuel characteristics and wildland fire dynamics via changes in stand structure and composition. To investigate the importance of MPB to past and future landscape dynamics, we used the mechanistic, spatially explicit ecosystem process model FireBGCv2 to quantify interactions among climate, MPB, wildfire, fire suppression, and fuel management under historical and projected future climates for three western US landscapes. We compared simulated FireBGCv2 output from three MPB modules (none, simple empirical, and complex mechanistic) using three focus variables and six exploratory variables to evaluate the importance of MPB to landscape dynamics. Results: We found that inclusion of MPB (empirical or mechanistic) in the simulations significantly changed past and future landscape dynamics and that the mechanistic MPB module had more cross-scale interactions that increased variability, and perhaps realism, of simulation results. We also evaluated impacts of fire and fuel management on MPB dynamics and found that fire suppression influenced fuel loadings more than MPB disturbance, but at a landscape scale, most fuel treatment programs did little to change fuel loadings, MPB dynamics, and burned area, except under high fire suppression. Conclusions: Synergistic interactions of climate, MPB, and wildfire catalyzed landscape-scale changes in vegetation distributions, fuels, and fire regimes in FireBGCv2 simulations. Models that simulate climate change on pine-domi-nated landscapes may be improved by including mechanistic MPB simulations to account for potentially important ecological interactions.


Introduction
Mountain pine beetle (Dendroctonus ponderosae Hopkins) (MPB) and wildland fire are two of the most significant natural disturbances in pine-dominated forests of the western US (Jenkins et al. 2008). As interacting and co-occurring disturbances, MPB outbreaks and wildland fire can cause highly visible, rapidly occurring, and persistent changes in landscape composition and structure . For example, bark beetle-caused tree mortality can affect the amount and configuration of surface and canopy fuels, with the potential for increased wildfire intensity and severity and subsequent, fireinduced changes in the availability of host trees for MPB (Meigs et al. 2015;Page et al. 2013). Beetle-caused tree mortality may also affect the amount of light, water, and nutrients available for overstory and understory plants, including exotics (Runyon et al. 2020), thereby altering species composition and stand structure (i.e., fuel characteristics) for years after an outbreak (Audley et al. 2020;Hansen 2014).
Wildland fire and MPB disturbances are responsive to patterns of climate and vegetation (Kane et al. 2017;Scheller et al. 2018). Climate, as a top-down driver, influences frequency and timing of MPB and fire disturbance (Bentz et al. 2010;Heyerdahl et al. 2001;Swetnam and Betancourt 1998), while vegetation and fuel characteristics act as bottom-up drivers, influencing MPB and fire severity and extent (Archibald et al. 2013;Fettig and Hilszczański 2015;Whitlock et al. 2008).
Anthropogenic climate changes-in particular drought and increasing temperatures-have been recognized as key drivers of increased wildfire (Littell et al. 2009) and bark beetle activities (Berner et al. 2017;Fettig et al. 2019). Warming temperatures may facilitate MPB population growth in some forests (Bentz et al. , 2019, and projected moisture deficit may increase drought stress, increasing vulnerability of host trees to MPB attack (Chapman et al. 2012;Kolb et al. 2016). Warming temperatures may also increase wildland fire frequency (Abatzoglou and Williams 2016;Westerling et al. 2006), size (Running 2006;Shive et al. 2014), intensity (Flannigan et al. 2009), and sometimes, severity (Miller et al. 2009). Predictions of warmer temperatures and increased drought suggest that the total area affected by large fires may increase in the coming decades (Williams, 2013). Acting independently or synchronously in space and time, wildland fires and MPB outbreaks can significantly influence forest structure, composition, and function; abruptly reorganize landscapes; and alter biogeochemical processes such as carbon cycling, water supply, and nutrient cycles (Edburg et al. 2012;Fettig et al. 2019;Kurz et al. 2008). Influences of MPB on wildfire occur in the context of other drivers, including topography, forest type, and disturbance history (Hicke et al. 2012). The complexity of these interactions among drivers and responses increases substantially when climate change and land management actions affect these processes, especially in this modern era of fire exclusion when approximately 96% of fires are suppressed each year .
The potential for increased wildland fire following MPB outbreaks has been observed since the early twentieth century (Evenden and Gibson 1940;Weaver 1943). Few studies, however, have addressed long-term influences of wildland fires on MPB disturbance regimes, and fewer still have examined reciprocal interactions of MPB and fire through several disturbance cycles (but see Parker et al. (2006)). Flammability and crown fire potential is highest in the first 3 years after MPB attack, when needles of killed pine trees have lower foliar moisture than green trees (Stalling et al. 2017) and are still attached to branches, maintaining aerial fuel continuity (Jenkins et al. 2014). Once standing dead trees have lost their needles and small branches (4 to 10 years after attack), active crown fire potential is much lower than in non-attacked stands, but increased fine surface fuel loads may result in higher surface fire rates of spread, flame lengths, and torching potential (Jenkins et al. 2014;Schoennagel et al. 2012). The highest fire hazard is assumed to occur decades after attack as a result of accumulation of heavy, large-diameter fuels (snags and large branches), dense regeneration, and increased wind speeds and drying of fuels resulting from the loss of pine canopy (Hansen et al. 2015;Jenkins et al. 2008).
Reciprocally, fire effects on MPB population activity are also time-dependent. Although fire-weakened trees are most susceptible to attack in the months after severe fire injury (Davis et al. 2012;Lerch et al. 2016;Parker et al. 2006;Powell et al. 2012), wildland fires can affect beetle activity indirectly over longer time periods by altering forest composition and structure (e.g., removing fire intolerant species) and providing increased water, light, and nutrients for surviving trees (i.e., growth release), thus influencing availability and vigor of suitable host trees (Fettig et al. 2007;Hessburg et al. 2005;Keeling and Sala 2012). In addition, occurrence of stand-replacing fires can reduce the likelihood of MPB attack until regenerating forests have attained a threshold diameter sufficient to attract beetles, especially when beetle populations are at relatively low levels (Kulakowski et al. 2003). Jentsch et al. (2021) found that fire helps dampen MPB outbreaks not only by removing host trees but also by altering the demographic structure of forest stands.
Although MPB was rarely considered a major contributor to post-fire pine mortality historically (Geiszler et al. 1980;McHugh and Kolb 2003;Six and Skov 2009), recent findings indicate that MPB may indeed be attracted to multiple fire-injured pine species and that the likelihood of mass attack and subsequent tree mortality was positively associated with several measures of fire injury, including crown scorch and cambial death (Davis et al. 2012;Lerch et al. 2016). Fires can therefore contribute to maintaining local MPB populations (Elkin and Reid 2004;Powell et al. 2012), but MPB response was found to pulse and recede within two years post-burn, thereby rarely generating sustainable outbreak populations with the potential for spread to adjacent stands (Davis et al. 2012;Powell et al. 2012).
Ecological change as a result of beetle outbreaks and wildfire also has important implications for wildlife populations (Saab et al. 2014;Saab and Powell 2005). Fireand bark beetle-caused tree mortality provides important habitat for wildlife dependent on snags and downed logs by providing nesting, roosting, and foraging substrate (e.g., several woodpecker species) . Disturbance-caused tree mortality, however, reduces occupancy of other species dependent on live foliage for nesting and foraging habitat (e.g., several songbird species) (Latif et al. 2021;Mosher et al. 2019).
Predicting landscape responses to MPB, wildfire, and climate is problematic because of the inherent complexity of their interactions (Loehman et al. 2020), the long time periods and broad spatial scales over which the disturbances co-occur (Archibald et al. 2013), and the important influences of historical and current land use on antecedent and future disturbance patterns and responses. Management actions, such as fire suppression and fuel treatments, confound projections of landscape responses to the interactions. Landscape simulation models provide an alternate method for assessing dynamics of large-scale interacting processes at appropriate temporal and spatial scales. Keane et al. (2011) developed the FireBGCv2 simulation modeling platform specifically to explore interactions of climate, vegetation, fire, and other disturbance agents. The original FireBGCv2 module that simulated MPB dynamics was empirical, simplistic, and overly generalized due to a general lack of comprehensive information to create a mechanistic simulation of MPB dynamics (Keane et al. 2011). Initiation of MPB outbreaks and outbreak severity were inferred from landscape conditions and a coarse evaluation of climate conditions, and MPB outbreaks were only held in check by winter temperature thresholds known to result in coldinduced MPB mortality. Tree mortality from MPB was simulated using a series of empirical regression equations dependent on tree diameter and developed over 50 years ago (Amman andCole 1983, Cole et al. 1985). Recently developed, mechanistic MPB models predict influences of thermal history on MPB cold tolerance (Régnière and Bentz 2007) and developmental timing and adult emergence synchrony Bentz 2009, 2014), both important factors in MPB population dynamics. We coupled these models with more up-to-date knowledge of MPB tree attack dynamics and stand susceptibility (Randall et al. 2011) within the FireBGCv2 platform to create a more complex, mechanistic MPB module.
In this study, our two major objectives were to (1) evaluate the importance MPB dynamics and its simulation rigor (no MPB, simple empirical simulation, and complex mechanistic simulation) on simulated landscape responses under future climates and land management actions and (2) assess the impacts of management actions on landscape-scale MPB, fire, and climate interactions. Simulations of these interactions over the entire western US are computationally intractable using the FireBGCv2 landscape model. Therefore, we selected three large, forested landscapes in the western US with varying fire regimes and densities of pine species that are known hosts to MPB. We then implemented a multifactorial simulation design that included three MPB factors (no MPB, simple empirical module, complex mechanistic module), two climate scenarios (observed historical weather, RCP 8.5 future climate), two wildfire suppression levels (no suppression, 96% suppression), and three levels of fuel treatment effort (no treatments, low and moderate intensity treatment).

Study areas
In our simulation experiment, we selected three western US study areas ( Fig. 1) that we felt represented a wide range of climate, vegetation, and fire regimes. The Ochoco watershed (OCH) in central Oregon is a 168,000 ha, dry mixed-conifer ecosystem consisting of low elevation ponderosa pine (Pinus Ponderosa Douglas Fig. 1 Map showing the location of the three landscapes (Ochoco-OCH, Elkhorns-ELK, Yellowstone Central Plateau-YCP) (in black) in the western US used in this simulation experiment. All three have pine forests with a history of MPB outbreaks. The Elkhorns have significantly less pines than the other landscapes. The ELK and OCH were initialized and parameterized specifically for this study while Clark et al. (2017) quantified the YCP landscape ex C. Lawson) and Douglas-fir Pseudotsuga menziesii (Mirbel) Franco) forests, middle elevation lodgepole pine (Pinus contorta Douglas) and subalpine fir (Abies lasiocarpa (Hooker) Nuttall) forests, and upper elevation whitebark pine (Pinus albicaulis Engelm.) forests, with historically frequent (25-to 100-year return intervals), low-severity fires. The Elkhorn watershed (ELK) in central Montana is a 140,000 ha, xeric mixed-conifer ecosystem of low elevation ponderosa pine and Douglasfir forests and high elevation lodgepole pine, whitebark pine, and subalpine fir forests, with moderately frequent (~ 20-50-year return intervals), low severity fires at lower elevations, and infrequent, high severity fires in the upper elevations. This landscape had the least area dominated by the MPB host pine species. The Yellowstone Central Plateau (YCP) in Yellowstone National Park, Wyoming, is a 175,000 ha high elevation landscape used in a previous simulation study (Clark et al. 2017) that consisted mainly of mature lodgepole pine forests with scattered late-successional subalpine fir and whitebark pine stands at the highest elevations. YCP mostly had an infrequent (> 150 year) fire return interval and stand-replacement fire regime.

FireBGCv2 model
FireBGCv2 is a mechanistic, individual-tree, gap model that is implemented in a spatial domain. The model was developed by integrating empirically derived deterministic functions with stochastically driven algorithms to approximate landscape and ecosystem behavior across time and space (Keane et al. 2011(Keane et al. , 1996 (Fig. 2). The executable FireBGCv2 program used for this project is available upon request and Keane et al. (2011) document model algorithms, assumptions, parameters, and validations.
The FireBGCv2 modeling platform integrates models of fire behavior and effects (FIRESUM (Keane et al. 1989), FOFEM (Reinhardt et al. 1997), and FARSITE (Finney 1998)) with the BIOME-BGC ecosystem model (Running and Hunt 1993;Thornton et al. 2002;White et al. 2000) and other modeled processes that describe tree and vegetation growth and mortality (Keane et al. 2011). The FireBGCv2 modules important in this study include simulation of stand-level management actions including timber harvest, tree planting, fuels treatments, and prescribed burning (Keane et al. , 2011; detailed fire algorithms that simulate natural ignitions based on Fig. 2 A schematic of the FireBGCv2 model. The right diagram shows important ecological processes and the scale that they are simulated including the mountain pine beetle (MPB), fire, climate, and fuel treatments (Land use). The middle diagram shows the five scales represented in FireBGCv2. Stand-level processes are simulated on a vignette of a stand called a plot. This plot is comprised of various ecosystem components (upper right panel). The flow of energy, carbon, water, and nitrogen is simulated across each of these components using mechanistic algorithms detailed in the Keane et al. (2011) document. Biophysical processes such as evaporation, transpiration, photosynthesis, and respiration are governed by five daily weather variables: precipitation, maximum and minimum temperature, radiation, and humidity weather, fuels, and topography; and spread fire across the simulation landscape based on vectors of wind and slope (Rothermel 1972), fire suppression, expressed as a proportion of naturally ignited fires that are suppressed at initial attack, an input parameter (Keane et al. 1996), and modeling of fire effects to vegetation based on fire intensity and soil heating (Reinhardt and Keane 1998).

Empirical mountain pine beetle FireBGCv2 module
The empirical MPB module initiated MPB population outbreaks based on percent landscape dominated by host pine species and two cold thresholds (− 20 and − 40 °C) that represent temperatures lethal to developing MPB. MPB-induced mortality of all pines larger than 20 cm DBH was then calculated using Amman and Cole (1983) algorithms as implemented in FVS (Crookston and Dixon 2005) for years with favorable, non-lethal freezing temperatures (see Keane et al. 2011 for further details).

Mechanistic mountain pine beetle module
The mechanistic MPB module implemented in FireB-GCv2 was based on a demographic model developed by Bentz (2009, 2014). This model simulated annual MPB population growth (R) as a function of temperature-dependent phenology as determined from daily maximum and minimum temperatures (Régnière et al. 2012). The demographic model incorporates the important role of tree defense in MPB success by connecting adult emergence timing (i.e., phenology) with population consequences based on emergence timing required for successful mass attacks on trees. High R values (R > 1) signify thermal conditions that support highly synchronized and seasonally appropriate adult emergence, and low R values indicate either unsynchronized emergence or emergence outside of a seasonally appropriate time window.
Cold temperature is considered a significant mortality factor in MPB (Weed et al. 2015). Acclimation of compounds that confer tolerance to extreme cold occurs over an extended time period as temperatures decrease in the fall prior to winter (Bentz and Mullins 1999). Here, we used a previously developed model to estimate annual probabilities of population survival as a function of daily maximum and minimum temperatures that influence larval cold tolerance from autumn to the following spring (Régnière and Bentz 2007).
Results from phenology, demographic, and cold tolerance models were linked by assuming that population growth would be most affected when a cold event results in extremely low offspring survival (also see Bentz et al. 2019). Several user input parameters for the mechanistic MPB model were required, including a threshold value at which cold mortality was introduced in a simulation year. If predicted survival from cold fell below the cold tolerance threshold, population growth predicted from the demographic model was reduced by the probability of cold mortality derived from the cold tolerance model. For FireBGCv2 results reported here, we assumed only extreme cold events would influence population growth. If predicted annual probability of survival from cold was ≤ 5%, population growth predicted from the module was multiplied by the probability of survival, thereby reducing population growth. No effect due to cold was applied if survival was > 5%.
The mechanistic MPB module does not predict initiation of an outbreak, but instead simulates population growth following initiation. To introduce MPB to a FireBGCv2 landscape, a background beetle mortality probability (set at 0.005) was used so every susceptible pine tree (diameter breast height (DBH) > 12.5 cm) in the stand had a 0.5% chance of dying from a random, endemic MPB attack in the current year. At the end of each simulation year, FireBGCv2 looped through all trees on a simulation plot and determined probabilistically whether each tree lives or dies due to beetles based on its probability of death. If a tree is chosen to be killed by beetles, based on this random draw, it was recorded as a "red" tree (i.e., dead tree with red needles) in the subsequent year and treated as a snag on the simulation plot, with its death due to MPB recorded.
Growth of MPB populations in a simulation year, once initiated, was also determined by the MPB mechanistic module, which combines effects of MPB phenology, demography, and cold tolerance models, as described above. Computed R-values represent ratios of new trees that will be killed in the stand in the current year based on the number of "red" trees killed the previous year, and daily maximum and minimum temperatures for that year. Daily weather was simulated at the site level; therefore, each site will have a different predicted R-value each year.
To implement tree mortality, stands with sufficient host tree species of adequate size were initially identified. To do this, we used a hazard rating system developed to estimate stand-level susceptibility of multiple pine species to MPB attack based on (1) average diameter of pines in the stand, (2) total basal area of all trees, and (3) proportion of the stand basal area comprised of pines (Randall et al. 2011). Based on these criteria, a single hazard rating is calculated on a scale from 0 to 27, with values 2 to 17 considered moderate, and values ≥ 18 considered high susceptibility. At the start of each simulation year, a "susceptibility" parameter was calculated using the hazard rating criteria for each stand. Another input parameter was the value above which a stand is considered susceptible. For the model results reported here, stands with a moderate rating (2-17) were considered susceptible to MPB attack for the current simulation year.
When the number of predicted red trees for a given year surpassed background beetle mortality threshold, the predicted number of red trees was then used as the population initiation for the next years' simulation. At low levels of MPB (< 10% of host trees attacked the previous year), only small host trees (12.5 < DBH < 18 cm) under the most stress (> 5 years of growth stress; see Keane et al. (2011)) were preferentially attacked. Tree stress in FireBGCv2 occurred when diameter growth dropped below a threshold value because of a combination of water stress, crowding, shading, and atypical temperatures. Each simulated tree had a counter recording the number of years that it has been under reduced growing conditions due to stress. At low MPB attack levels (< 10% of host trees attacked the previous year), the probability of death due to MPB assigned to each tree in the simulation plot for a susceptible stand was weighted by (1) total number of trees predicted to be killed according to the R-value and (2) stress level of each small tree as compared to all other available small trees on the stand. With increasing MPB population size (10-25% of total pines in the stand were attacked the previous year), stress was no longer considered an important factor in tree attack. At the highest levels of MPB presence (> 25% of total pines in the stand the previous year were attacked), only the largest host trees (DBH > 20 cm) were targeted. At the end of a simulation year, the fate of each tree on the simulation stand was determined by a random draw compared with its individual probability of death due to MPB. Due to the structure of FireBGCv2 (Fig. 2), probabilities of MPB death for individual trees were determined at the stand level, although random draws determine the fates of individual trees.

Simulation experiment
Our factorial simulation experiment consisted of replicated model runs for all combinations of four factors, with two or three levels for each factor (Table 1). Factors and levels included: (1) Mountain pine beetle (MPB)-no modeled MPB, MPB simulations using the simple empirical module, and MPB simulations using the complex mechanistic module as discussed in detail above. (2) Fire suppression level (FSL)-no suppression (historical fire regime with no fires suppressed) and full suppression (96% of fires suppressed which represents the current level for each of the three landscapes).
(3) Fuel treatment effort (FTE)-A "fuel treatment" is simulated by the removal of trees to reduce canopy fuels and a prescribed burn to reduce surface fuels (details discussed later in this section. There were three levels of fuel treatments-none, a passive, low-level fuels program (~ 1% landscape treated annually), and an aggressive, moderate-level fuels program (~ 5% landscape treated annually) (4) Climate (CLI)-an historical climate and a projected future climate under RCP8.5 scenario.
These last three factors (FSL, FTE, and CLI) were included to evaluate if differences in output from MPB modules were also evident for past and future climates, historical and future fire regimes, and after implementation of fuel treatments. The nine response variables in Table 1 were output at the end of every simulation decade over a 400-year simulation period for five replicates per factorial combination resulting in 200 observations per response variable per scenario. The short 400-year time span relative to some fire return intervals (e.g., YCP) was chosen because of a lack of adequate computing resources; it often took weeks to simulate landscape dynamics over this time span and simulation times for longer time spans increase exponentially.
Historical climate data were derived from a nearby weather station for each simulation landscape (see Table 1) using techniques from previous FireBGCv2 projects (Clark et al. 2017;Keane et al. 2019Keane et al. , 2011. Weather data included maximum and minimum temperatures and precipitation. Temperatures and precipitation, and also solar radiation and humidity, were extrapolated across each landscape using the model MTCLIM (Hungerford et al. 1989). Future climate was derived by offsetting daily, historical climate variables by offsets (Delta method) as described in (Clark et al. 2017;Loehman et al. 2010) (e.g., we added a temperature value computed from climate projections, called an offset, to each daily temperature value). Offsets were derived from CNRM-CM5 (National Centre of Meteorological Research, France) 30-arc-second resolution (~ 800-m) climate model, the highest ranked model overall for the three landscapes (Rupp et al. 2013), for the Coupled Model Intercomparison Project Phase 5. We used the 8.5 Representative Concentration Pathway (RCP85) data, an emissions scenario that predicts the highest average annual temperature increases among various possible emissions scenarios (about 6 °C above pre-industrial levels) because it represents an extreme climate scenario (Peters et al. 2012). The Delta method is commonly used when dealing with daily weather data because future daily projections poorly represent climate variability (Keane et al. 2015b).
For simulations using historical climate, we continually repeated the 60-100 + year weather record for the full 400-year simulation, and for future climate simulations, we also cycled the weather years but ramped weather so that all offsets were fully implemented after 100 years, and then kept offsets constant for the remaining simulation time (Holsinger et al. 2014).
A fuel treatment in this study was a combination of simulated silvicultural cuttings to remove all fire-sensitive, shade-tolerant species down to 4.0-cm diameter at breast height followed by a low intensity prescribed burn (~ 400 kW m −1 fireline intensity) to reduce slash and kill understory seedlings and saplings. FireBGCv2 simulated treatments only in stands that were (1) above a minimum basal area (10 m 2 ha −1 ); (2) contained fire-sensitive species (Table 1); and (3) had high fuel loadings (> 5 kg m −2 ). There were many simulation years when few stands met these criteria so treatment level targets were not reached .
Nine response variables (Table 1) were output every 10 years for a simulation time series of 41 observations per variable per run, and then performed five replicate Table 1 Details of the simulation experiment used in this study. We simulated four factors on three landscapes over 400 years with five replicates and output generated every 10 years (400/10 × 5 = 200 observations per factor combination). No other disturbances but fire and MPB were simulated. We started all simulations using the initial conditions of the landscapes circa 2010. While the model generated over 33 output variables, we used only nine in this study. Variables with asterisks are the focus variables in this paper for brevity   (Table 2), which dictates tree species available to inhabit or migrate into the landscape during simulations.

Statistical analysis
We compared responses among scenarios by fitting generalized linear mixed models that included the fixed effects of MPB, FSL, FTE, CLI, landscape, and year, and a first-order autoregressive covariance structure to account for temporal autocorrelation (Box et al. 1994;Pinheiro et al. 2007). All two-way interactions were included in the model. Separate models were fit for each landscape and response variable. We used the significance of each factor (P ≤ 0.05) and standard boxplots to assess differences between each MPB/FSL/ FTE/CLI factorial scenario. Of the nine response variables, we focused on the differences across scenarios using three ecological variables: (1) basal area of MPB host pine species (BAP), (2) proportion of landscape in MPB host species (HOS-TAREA), and (3) proportion of the simulation landscape burned annually (PLBURN) (Figs. 3, 4, and 5). We selected these response variables because we judged that they were the most important for understanding influences of MPB, fire, and climate interactions. The remaining six variables were used as explanatory variables and their simulated results are shown in supplementary material (Figures S1,S2,S3,S4,S5 and S6).
Recognizing that differences across the four factors and their levels might depend on combinations of all nine response variables ), we used principal component analysis (PCA) on the nine response variables in Table 1 to reduce the dimensionality of the data and to evaluate differences among scenarios and levels for each landscape. Only the scores of the first principal component (PC1) were evaluated in our PCA analysis to detect significant differences (P ≤ 0.05) using methods described in the previous paragraph (Fig. 6). The PCA analyses were performed in R using built in functions.

Sensitivity of ecological responses to the complexity of MPB sub-models
The significance of MPB module (none, empirical, mechanistic) on predictions of the three focus and six explanatory response variables depended on the landscape simulated and FSE, FTE, and CLI scenarios (Table 3, Figs. 3, 4, 5, and 6, and Figures S1,S2,S3,S4,S5 and S6). First and most importantly, MPB module type influenced predicted BAP, HOSTAREA, and PLBURN on the OCH landscape, BAP and HOSTAREA on the YCP, and only BAP on the ELK landscape (Figs. 3, 4, and 5). In both historical and future climates, predicted BAP and HOSTAREA were lower (> 30%) using the more complex mechanistic MPB module compared to simulation runs with no-MPB included (Figs. 3 and  4). Low values for BAP and HOSTAREA were direct consequences of greater predicted MPB-killed BA by the mechanistic MPB module, particularly on the YCP landscape (Fig. S5). Significant differences (> 50%) between the mechanistic and the other two MPB modules for the YCP landscape were supported by values of PC1 (Fig. 6).
Predictions of BAP and HOSTAREA, but not PLBURN, were also the most variable in historical climate simulations using the mechanistic MPB module (see width of box and whiskers in Figs. 3 and 4). In future climates, however, variability in predicted Table 2 The tree species (scientific and common names) that were included in the FireBGCv2 simulations; species names in bold (F) indicate that the species is fire-adapted while non-bold species are fire-sensitive. The asterisk indicates the species is a host of mountain pine beetles (MPB) Historical climate simulations with no-MPB resulted in greater BAP (> 20%) on both the OCH and ELK landscapes, and greater HOSTAREA on the OCH, but no differences between no-MPB and empirical MPB model simulations were observed on the YCP (Figs. 3 and 4). In future climate scenarios, however, the no-MPB simulations predicted greater BAP and HOSTAREA on all three landscapes, although variability in model outcomes was also greater, particularly for the YCP landscape. In the historical climate simulations, PLBURN was greater on the OCH landscape compared to the ELK and YCP, regardless of the MPB model type. This prediction was only found for no suppression treatments in historical climate simulations. In future climate scenarios, FSE, FTE, and MPB model type influenced predicted values of PLBURN.

Influence of management on MPB, fire, and landscape interactions
We found that the aggressiveness of fuel treatments (FTE) had relatively minor effects on most response variables across all three landscapes (Figs. 3, 4, 5, and 6), even though some differences were significant at the FTE factor level (Table 3). Moreover, FTE had little effect on MPB-killed basal area, especially when fires were allowed to burn on the landscape ( Figure S5). In future climate scenarios, PLBURN was predicted to be greater on the OCH than both other landscapes and was also greater in no suppression compared to suppression simulations. PLBURN was also greater in future climate scenarios for the YCP and ELK landscapes (~ 25% greater), but only in no suppression simulations (Fig. 5). Simulations with no-MPB and no suppression resulted in less PLBURN on the OCH landscape as the level of FTE became more aggressive (Fig. 5). Neither FSE, FTE, nor MPB model type influenced predicted PLBURN on the YCP and ELK landscapes. Relative to historical simulations, HOSTAREA and BAP were greater in future climate scenarios on the OCH for  Table 1). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range the mechanistic MPB model and no-MPB simulations, relative to the simulations with MPB included (Figs. 3 and  4). On the YCP landscape, variability in HOSTAREA and BAP increased in the future climate simulation runs that included MPB.
Forest fuels, represented by canopy bulk density (CBD; Figure S2), coarse woody debris (CWD; Figure S2), and fine woody debris (FWD; Figure S3) were similar across all three MPB model scenarios. FTE did not cause major reductions in the wildland fuel output variables. FSE had a much greater influence on fuels and fire dynamics than fuel treatment. In both historical and future climate scenarios, HOSTAREA on the ELK landscape was predicted to be very low, regardless of the MPB model type, FSE or FTE.

Discussion
The key finding from this simulation effort is the importance of including MPB dynamics in landscape models that simulate compositions and structures of current and future western US landscapes that contain MPB pine host species. Stand, landscape, and regional models that ignore impacts of this important disturbance are missing meaningful interactions of a major ecological disturbance (Kurz et al. 2008;Liu et al. 2011). We found significant differences and increased variability in simulated projections of landscape composition and basal area when either the empirical or mechanistic MPB module was used, compared to simulations with no-MPB. However, the mechanistic MPB module consistently predicted less BAP and HOSTAREA in historical and future climate scenarios relative to no-MPB but did not influence the third focus variable PLBURN. The mechanistic MPB module also predicted the greatest MPB-killed basal area. In addition, we found greater variability in predictions of BAP and HOSTAREA in historical simulations when the mechanistic MPB module was used. This high variability is a result of more ecological interactions represented in the mechanistic MPB module as it was integrated within FireBGCv2, including direct factors such as tree stress, climate, fire, and forest structure and composition (Keane et al. 2011). By contrast, variability in BAP and HOSTAREA was greatest in future climates when the empirical MPB  Table 1). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range module was used, especially for the YCP landscape. The likely reason for odd model behavior is that the empirical MPB module was largely based on statistical algorithms and cold temperature thresholds that were inappropriate because future weather was outside the bounds of data used to create the original equations and estimate important thresholds.
The mechanistic MPB module is likely to provide more realistic projections than the empirical model under future RCP8.5 climates because that module (1) includes an explicit representation of the influence of weather on MPB population success, (2) contains numerous interactions with other ecosystem variables (e.g., tree size and abundance, previous year MPB-caused tree mortality), (3) integrates several interactions in tree mortality algorithms, and (4) is not limited by the narrow ranges of empirical data used in the empirical MPB module. MPB is an important disturbance process in western US pine forests, and our results suggest it is the interactions of MPB population dynamics, climate, and their effects with other ecological processes that dictate landscape model projections (Bockino 2008;Jentsch et al. 2021;Keane et al. 2015a;Loehman et al. 2017).
Another important finding is that there are great differences among landscapes for most response variables. Simulation results for the ELK landscape, for example, differed from the OCH and YCP in terms of predicted proportion of pine on the landscape in both historical and future climates. This could be because the inventory data used to initialize FireBGCv2 for ELK contained mostly Douglas-fir and subalpine fir, and pine (MPB host) species were rare (ponderosa, lodgepole, and whitebark pine). As a result, differences across scenarios for the ELK landscape were difficult to detect because of low pine populations. Results from the YCP landscape often contained greater variation for all nine variables potentially due to the long fire return intervals (> 250 years) and the dominance of colder climates which favored lodgepole pine. Frequent fire regimes as found on the ELK and OCH landscapes tend to create shifting mosaics of vegetation communities that are surprisingly stable ). The OCH Fig. 5 Boxplots of FireBGCv2 simulation results of the annual percent of landscape burned (PLBURN %) showing the influence of the influence of mountain pine beetle (MPB) simulation (none = red, empirical = blue, mechanistic = green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP), treatment effort (none, passive, aggressive), and fire suppression level (none, suppression) under A current climate and B future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range landscape had a warmer environment and was dominated by ponderosa pine with higher annual % landscape burned, and a lower proportion of the historical landscape with fire adapted species compared to YCP.
The simulation of MPB dynamics, regardless of the module type, had significant impacts on burned area (PLBURN), but only for the OCH landscape and only under fuel treatments. The OCH landscape had a frequent fire regime that likely can be replicated by fuel treatments. In general, wildfire had a greater effect on wildland fuels in our simulations by decreasing CBD (Figure S2), CWD ( Figure S3), and FWD ( Figure S4) in both historical and future climates, especially for the OCH and ELK landscapes. The short simulation time span (only 400 years) was not long enough for the long fire return intervals for the YCP landscape (~ 250 years), resulting in fewer fires that probably masked impacts of wildfire on canopy and surface fuels. MPB-caused mortality often results in the deposition of canopy fuel from dead pines on the ground along with the dead canopy fuel killed by fire (Jenkins et al. 2012), but we found little change in simulated surface fuels as a consequence of MPB-caused tree mortality. These results are supported by Stalling et al. (2017) who found little change in forest floor fuels within a decade following MPB outbreaks. Under future climates, fire suppression was predicted to result in greater BAP on the YCP landscape and more MPB-killed basal area compared to no suppression actions in future climates. Interestingly, the opposite effect was seen in the OCH and ELK landscapes where greater MPB-killed basal area was predicted when there was no fire suppression. Fire suppression on the YCP allowed established pines to grow into trees of sufficient size to be vulnerable to MPB, thereby resulting in landscape homogenization favoring pine and, in turn, creating the potential for high severity outbreaks of MPB.
Another surprising result was minor influence of fuel treatments on most simulated response variables, especially their low impact on tree mortality due to MPB or fire. Aggressive fuel treatments reduced total and pine basal area and percent landscape burned on the OCH without MPB, but this was the only case. Keane et al. (2019) found that fuel treatments are only effective at  Table 1 in the PCA analysis of FireBGCv2 simulations showing the influence of the influence of mountain pine beetle (MPB) simulation (none = red, empirical = blue, mechanistic = green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP), fuel treatment effort (none, passive, aggressive), and fire suppression level (none, suppression) under A current climate and B future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range restoring landscapes if implemented under high fire suppression levels, which was validated in this effort using different landscapes. Wildland fire, whether as a controlled or uncontrolled wildfires, tends to do most of the work of fuel treatments but with significantly less resources . Fuel treatments, regardless of level of effort, didn't indicate any effect on MPBcaused tree mortality except for the OCH landscape Table 3 Abbreviated results from the generalized linear mixed model analyses that include fixed effects and two-way interactions of mountain pine beetle module type (none, empirical, mechanistic) (MPB), fire suppression level (FSL), fuel treatment effort (FTE), climate (CLI), and simulation output year (YEAR) (n = 200) (see Table 1 for factor details). Results are reported for all three landscapes (Fig. 1). The focus variables are basal area of pine species (BAP), percent of the landscape dominated by MPB host species (HOSTAREA), and percent of the landscape burned each year (PLBURN), while the explanatory variables are canopy bulk density (CBD), coarse woody debris (CWD), fine woody debris (FWD), basal area killed by MPB (BAMPB), and percent of the landscape in fire adapted species (PLFAS) (see Table 1). DF is degrees of freedom and values of each variable were output every 10 years for 400 years for five replications (n = 5 × 40 = 200). The p-value is reported for each response variable with numbers in bold indicating statistical significance (p < 0.05). NA indicates that fixed effect was not included to avoid overly dispersed models without MPB simulation. This is partially a result of wildfire effects mitigating the limited influence of fuel treatments over long time periods.

Limitations
The current implementation of MPB modules within the FireBGCv2 simulation platform did not incorporate dispersal of MPB outbreaks across stands because of the additional computational power required to integrate MPB dispersal models (Powell and Bentz 2014 Results of our statistical analysis of the simulation output also have shortcomings. Our analyzed output was generated every 10 years over the 400-year time span for five replications giving us a large 200 observations' data set. This large number of observations provides great statistical power, but as White et al. (2014) note, this statistical power is easily created in simulation projects by simply outputting more data. Because of this, it may be best to de-emphasize the importance of the p-values in Table 3 and put more confidence in the visual comparison of the box and whisker plots (Figs. 1, 2, 3, 4, 5, and 6 and S1, S2, S3, S4, S5 and S6). The variabilities of output variables over the 200 observations are also important. When employing mostly deterministic models, such as FireBGCv2 (the stochasticity in FireBGCv2 is only 3-5% of the variance) (Keane et al. 2011), the magnitude of the variance becomes an important indicator of the realism of simulated results. In general, increased variability often results from increased interactions causing non-linear model behaviors (Keane 2019). Obviously, the simulated variability is only a rough approximation of the actual ecological variability, but its comparison across scenarios provides an indication of possible realism in model results. For example, the same amount of variability across the box and whisker plots may indicate that the model may not be simulating alternate interaction pathways.

Summary and conclusions
Provided here are simulation results for three landscapes using FireBGCv2 integrated with empirical and mechanistic MPB modules that demonstrated the importance of including MPB when simulating landscape dynamics. As expected, adding MPB effects decreased total and pine basal area, regardless of fire suppression in both historical and a future climate, although there was little effect on the annual percent of landscape burned. Fire suppression increased large pine abundance and reduced annual percent of landscape burned but had varying effects on MPB killed basal area depending on the landscape and climate scenario. Simulations using integrated MPB models underscore the significant influence of MPB disturbances in these landscapes and highlights future stability in species composition with and without management. We anticipate integrating wildlife population models with these simulations to evaluate competing risks of MPB, wildfire, and fuels treatments on future habitat quality (based on forest structure and composition) and the distributions of selected wildlife species (e.g., Saab and Powell 2005).
Additional file 1: Figure S1. Basal area of all trees (BA; m2 ha-1) boxplots of FireBGCv2 results showing the influence of the detail of mountain pine beetle (MPB) simulation (none=red, empirical=blue, mechanistic =green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP) and fire suppression level (none, suppression) under (A) historical climate and (B) future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range. Figure S2. Canopy bulk density (CBD; kg m-3) boxplots of FireBGCv2 results showing the influence of the detail of mountain pine beetle (MPB) simulation (none=red, empirical=blue, mechanistic=green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP) and fire suppression level (none, suppression) under (A) historical climate and (B) future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range. Figure S3. Coarse woody debris loading (CWD; kg m-2) boxplots of FireBGCv2 results showing the influence of the detail of mountain pine beetle (MPB) simulation (none=red, empirical=blue, mechanistic=green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP) and fire suppression level (none, suppression) under (A) historical climate and (B) future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range. Figure S4. Fine woody debris loading (FWD; kg m-2) boxplots of FireBGCv2 results showing the influence of the detail of mountain pine beetle (MPB) simulation (none=red, empirical=blue, mechanistic=green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP) and fire suppression level (none, suppression) under (A) historical climate and (B) future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range. Figure S5. Basal area of all trees killed by mountain pine beetle (BAMPB; m2 ha-1) boxplots of FireBGCv2 results showing the influence of the detail of mountain pine beetle (MPB) simulation (empirical=blue, mechanistic=green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP) and fire suppression level (none, suppression) under (A) current climate and (B) future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range. Figure S6. Proportion of the landscape composed of fire adapted species (PLFAS; proportion) boxplots of FireBGCv2 results showing the influence of the detail of mountain pine beetle (MPB) simulation (none=red, empirical=blue, mechanistic=green) stratified by landscape (Ochoco-OCH, Elkhorns-ELK, and Yellowstone Central Plateau-YCP) and fire suppression level (none, suppression) under (A) historical climate and (B) future climate (RCP8.5). The horizontal bold line in each boxplot is the median, the box top and bottom are the 25th and 75th percentiles respectively, and the whiskers extend to the smallest or largest values but no further than 1.5*Interquartile range. Figure S7. PCA results for YCP landscape from FireBGCv2 simulations of historical climate, full fire suppression, mechanistic mountain pine beetle module, and no fuel treatments. Circle shows the range that contains 68% of the variation. Names with arrows indicate variable importance in PC1 and PC2 scores. Variable acronyms are defined in Table 1.