1. Introduction

There is a challenging dichotomy in the resource management community. On the one hand, there is a need to know and understand fundamental physical processes at a very specific and local level. On the other, there is also a need to know the larger context within which the system works. This is particularly true for integrating climate change into natural resource management practices. Ecological restoration and management actions require an understanding of local conditions (such as species competition, disturbance conditions, micro topography, and a given jurisdiction’s management options), but climate impacts are often measured and modeled at a coarser scale and imply the need to consider the wider landscape. Recent work has called for integrating various methods in order to overcome limitations of individual tools and the challenges associated with climate change adaptation [1].

Species distribution modeling (SDM) has proven useful for estimating the potential impacts of climate change and other abiotic factors on species distributions. Many SDMs are popular due to their ability to generate complex predictions without high computational requirements; however, there has been increasing recognition of the limitations of SDMs related to interspecific interactions, dispersal, equilibria, stochastic events, and fundamental niche space [2,3]. In addition to their inherent limitations, there are some specific questions that are not appropriate or possible within the framework of correlative modeling but may be of keen interest to natural resource managers. For example, what would happen to a species if the budget for a given restoration project was doubled or cut in half?

State-and-transition simulation models (STSMs) provide a spatial and quantitative framework to explore “what if” scenarios, which can be used to explore both management options and evaluate the sensitivity of the system to specific parameterizations or assumptions [4]. STSMs have been recognized as useful tools for incorporating the effects of multiple disturbances, biotic interactions, and management scenarios, but lack statistically robust techniques for relating climate data to species distributions. Together, STSMs and SDMs are well-equipped to account for the various facets of natural resource management challenges, at both regional and local scales, and combine them in a spatial simulation framework.

Our case study of whitebark pine in the Greater Yellowstone Ecosystem serves as a proof-of-concept for combining stochastic simulation models (STSMs) and correlative models (SDMs), and sets the stage for future exploration of management options. This research offers two novel contributions: 1) accounting for climate change and other dynamics through a combination of SDMs and STSMs; and 2) outlining an approach for validating STSMs based on recent developments in agent-based modeling, which remains an active area of exploration.

2. Background

2.1. Species distribution models

Species distribution models (SDMs) refer to a variety of empirical methods that associate species with the attributes of their preferred bioclimatic environment, also sometimes referred to as “habitat” [5,6]. SDMs evolved over the last several decades from our conceptual understanding of niche theory combined with the increased availability of geospatial datasets and geographic information systems with which to quantify bioclimatic attributes and their variation across space. SDMs typically represent individual species presence (and/or absence) at the spatial resolution of commonly available climate data (usually ~ 1 km or coarser) at temporal resolutions of decades or longer (often 30 year averages). More recently, interest in the ecological impacts of climate change has spurred rapid development and growth in application of SDMs because they are relatively easy to train on contemporary conditions and apply to projected future conditions. Future SDM projections often support commonly held assumptions about climate change and species extinctions and range shifts, although their temporal transferability (i.e., ability to be calibrated for one time period and make predictions for another) is rarely considered [7].

There are numerous challenges to producing actionable science using SDMs including that: they often represent unknown mechanistic relationships between species and their environments [8]; empirical methods do not always capture the direct effects of climate on species in ecologically meaningful ways [9]; there is often a mismatch of spatial and temporal resolution between SDMs and management action [10]; and SDMs rarely consider climate change impacts to communities of species and functional types whereas management is often charged with preserving assemblages of species and ecosystem function [11]. Numerous refinements of SDMs have been proposed [12]; however, the above challenges suggest that there remains a need to develop new methods that better produce scientific information that is useful to natural resource managers who are engaged in climate adaptation planning.

2.2. State-and-transition models

State-and-transition models originated as conceptual models that represented groups of vegetation communities and the shifts between them [13]. The definition of states often depends on the modeling objectives and data availability, but can generally be thought of as suites of vegetation communities that have distinct functional groups, ecosystem processes, and structure [14]. Transitions include natural events, management interventions, or a combination of both [15]. State-and-transition models are typically represented using box and arrow diagrams, in which boxes or nested boxes represent vegetation phases and states, and arrows represent the transitions between them.

These conceptual models remain at the heart of more recent quantitative computer-based state-and-transition models, or state-and-transition simulation models (STSMs, reviewed by Daniel and Frid [4]). In STSMs, transitions can be deterministic (e.g., growth, aging) or probabilistic (e.g., fire, invasion), and can be aspatial or spatially explicit. Spatially explicit STSMs are analogous to joint cellular automata-Markov models [16], hybrid Markov-cellular automaton models [17], and spatio-temporal Markov chains [18]. Bestelmeyer et al. [19] note that spatial context is important for conceptual state-and-transition models because spatial dynamics such as contagion, feedbacks between patches, spatial patterns in historical legacies, and variation in soils, topography, and climate can all affect the likelihood and location of transitions. These spatial dynamics can similarly be incorporated into spatially explicit STSMs.

STSMs can be used to compare and evaluate different resource management scenarios [20,21,22,23], and can incorporate climate effects [24,25,26,27,28]. Here we present a novel approach to incorporating the effects of climate on vegetation into STSMs that is spatially explicit and based on correlative models of habitat suitability.

2.3. Study area

Whitebark pine (Pinus albicaulis) is an iconic member of the subalpine forest community in the western North American mountains and is highly valued for the ecosystem services it provides. It is considered a keystone species because it has strong influence on ecosystem properties. Whitebark pine (WBP) is able to tolerate the harsh conditions typical of mountainous terrain and creates forest in high-elevation locations that would otherwise be shrubs and herb. By colonizing open subalpine patches, it provides safe sites for subalpine fir (Abies lasiocarpa) and other species resulting in greater forest cover. In addition to providing habitat for many wildlife species the unusual reproductive strategy of producing large crops of cones with large nutritious seeds every few years provides a rich food source to the threatened Grizzly bear and other wildlife species. In addition to these positive effects on high-elevation ecosystems, the species has strong social values. The unique umbrella shape, large size, and grizzled appearance developed over its centuries-long lifespan evokes strong emotional appeal to people visiting these mountain haunts.

The Greater Yellowstone Ecosystem (GYE), which includes Yellowstone National Park, Grand Teton National Park, and a number of state and federally managed forests, is a mid- to high-latitude region in the Northern Rocky Mountains of western North America. Conifers are dominant in the range, with forest types composed of Pinus contorta, Abies lasiocarpa, Pseudotsuga menziesii, Pinus albicaulis, Juniperus scopulorum, Pinus flexis and Picea engelmannii, although the deciduous hardwood Populus tremuloides, is also wide spread. Plateaus and lowlands are dominated by species of Artemisia tridentata and open grasslands of mixed composition. The GYE study area encompasses 150,700 km² with an elevational gradient from 522-4,206 m that represents 14 surrounding mountain ranges [29].

WBP is present in the GYE from below 2,100 m to nearly 3,300 m, an elevation range that spans the montane and subalpine zones [30]. WBP is subdominant to other conifer species in the lower portion of its distribution and is dominant in many locations at upper treeline [31]. The WBP population in the GYE has been particularly hard hit in recent years: in some areas whitebark pine mortality has exceeded 95% of cone bearing trees (DBH > 15 cm) [32] due to factors related to warming climate, mountain pine beetle (Dendroctonus ponderosae), and an exotic fungal pathogen (Cronartium ribicola) associated with white pine blister rust [33]. Beyond the current forest die-off, resource managers are concerned because the area of suitable habitat for WBP in the GYE is projected to decline dramatically in the coming century due to projected climate change [29,34,35,36]. Consequently, the US Fish and Wildlife Service listed the WBP on the U.S. candidate species list [37]. Yet questions remain as to how best to manage WBP under these various threats.

Figure 1. Study area, highlighting the approximate extent of the Greater Yellowstone Ecosystem.

DownLoad: Full-Size ImgPowerPoint

3. Methods

3.1. Species distribution modeling

For the SDM portion of the research, we used the VisTrails: SAHM software package to preprocess data, select predictor variables, compare model algorithms, produce model diagnostics, and capture data and workflow provenance. VisTrails is a free open source scientific workflow software [38] that has been customized for SDM using a set of add-on tools that comprise the Software for Assisted Habitat Modeling (SAHM) [39].

Field observations of adult tree presences and absences were compiled from the Forest Inventory and Analysis (FIA) program, Whitebark/Limber Pine Information System [40], and long term monitoring plots established by the National Park Service Greater Yellowstone Inventory and Monitoring Network [41]. “Adult” class WBP were selected for modeling based on a recorded diameter at breast height (DBH > 20 cm). WBP within central Montana are reported to reach 100 years of age at approximately 8-12 m in height with DBH 15-20 cm. Given previous silvicultural studies, it was assumed that 20 cm DBH for WBP represent adult class individuals for the GYE, with potential to reproduce. First, 2,545 WBP observations from the Forest Inventory and Analysis (FIA) [42] program were assembled. FIA plots are located on a regular gridded sampling design with one plot at approximately every 2, 500 forested hectares, with swapped and fuzzed exact plot locations within 1.6 km to protect privacy. Gibson et al. [43] found that model accuracy was not dramatically affected by data fuzzing, but to provide the most spatial accuracy, this study culled FIA field points where measured elevation were > 300 m different than a 30 m USGS digital elevation model. In order to generate probability of occurrence surfaces for competitor species, separate species distribution models for spruce-fir [1] and lodgepole pine were constructed using 2,489 FIA plots.

Probability surfaces were generated using a random forest algorithm. Random forest is an ensemble learning technique that generates independent random classification trees using a subset of the total predictor variables and classifies a bootstrap random subsample of the data. The predictor variables were 30-year means (1950‒1980) of Parameter-elevation Regressions on Independent Slopes Model (PRISM) variables and Thornthwaite-based dynamic water balance model variables [29]. We used a maximum correlation filter of 0.75 (all pairs of variables) to avoid collinearity issues, and physiologically relevant variables were determined by expert analysis and literature review [29]. A total of eight predictor variables were used to determine the probability of WBP occurrence: minimum temperature January, vapor pressure deficit March, precipitation April, snow water equivalent April, maximum temperature July, actual evapotranspiration July, potential evapotranspiration August, precipitation September. Atmospheric CO2 concentrations were not in the list of available predictors, but are noteworthy because they could affect tree physiology and response to warming [44].

Accuracy for the model was evaluated by calculating the receiver operator characteristic curve (AUC). Although AUC alone does not provide an explicit description of commission and omission error rates within a model, it does serve as an index of how likely a model can discriminate a presence versus an absence [45]. As a general rule of thumb, an AUC of 0.5 indicates performance no better than random and 1.0 indicates perfect model prediction. AUC measures above approximately 0.7 are generally considered to be good and above 0.9 excellent [46]. After fitting, our WBP model reported an AUC value of 0.94, displaying high specificity and sensitivity. In addition to examining the AUC for WBP, the out-of-bag (OOB) error estimation was also examined and found a rate of commission (13.1%) and omission (10.9%) for ensemble bagged trees, demonstrating low test error. Evaluation of AUC for the other species displayed good/excellent skill for present day climate conditions.

A post-hoc comparison of habitat niche fit with seedling class WBP (< 2.54 cm DBH, n = 497) was also constructed, displaying a spatial distribution for WBP analogous to adults (Figure 2). Binary classification of presence and absence based on a probability threshold, where sensitivity and specificity were equal, generated near equivalent distribution maps for adults and seedlings despite seedlings presenting lower occurrence probabilities as a result of lower empirical sample prevalence. Comparison of predicted presence distributions across predictor variables again verified similar environmental gradients for both adults and seedlings (Figure 3). This post-hoc analysis rationalizes use of the “adult” fit distribution probabilities under future climate for state and transition simulations as potential suitable habitats for dispersed recruit colonization (details below) due to its more complete empirical representation of the population on the landscape.

Figure 2. Modeled suitable bioclimatic envelope for adult WBP (938 presence, 1631 absence; accuracy rate of about 93%) (left panel), and for seedlings (497 presence, 2072 absence; accuracy rate of about 84%) (right panel).

DownLoad: Full-Size ImgPowerPoint

Figure 3. Binary classification of modeled presence for WBP adults (blue) and seedlings (red) display similar environmental gradient distributions.

DownLoad: Full-Size ImgPowerPoint

Following model fitting and applying the model to a historic climate period (1950‒1980 and 1980-2010, respectively), we projected our SDM under two global climate models (GCMs) and two carbon concentration scenarios. Using a Bias-Correction Spatial Disaggregation (BCSD) approach, an archive of 9 statistically downscaled CMIP5 climate projections for the conterminous United States at 30-arc-second spatial resolution was assembled by the NASA Center for Climate Simulation NEX-DCP30. For this analysis, two GCMs that represented the greatest and least change in total WBP habitat in the GYE by the year 2099 (CNRM-CM5 [46] and HadGEM2-AO [47], respectively), and two representative concentration pathway (RCP) scenarios were used to project WBP occurrence probabilities. RCP 4.5 was the first, representing increased radiative forcing until stabilization of greenhouse emissions between 2040 and 2050. RCP 8.5 was the second, representing the “business as usual” scenario, with uncontrolled radiative forcing with stabilization by 2099. Using this approach, we sought to demonstrate the “bookends” of range for projected WBP probable habitat under the high uncertainty and variability of global climate models [48]. Occurrence probability surfaces were constructed for the 2040, 2070, and 2099 climate projections.

3.2. State-and-transition modeling

We implemented the STSM in the ST-Sim modeling platform [49]. ST-Sim’s graphical user interface streamlines the process of defining distributions and probabilities, managing model inputs and outputs, viewing results, and creating and comparing alternative simulation scenarios.

Within ST-Sim, we defined state-and-transition pathways for each potential vegetation type (or stratum) in the model. Strata and state classes were defined based on vegetation dynamics models developed for the area as part of the Landscape Fire and Resource Management Planning Tools Project (LANDFIRE) [50]; specifically, we created three forest types (WBP, lodgepole, and spruce-fir) based on the dominant species in three LANDFIRE biophysical settings (BPS) (northern rocky mountain subalpine woodland and parkland, rocky mountain lodgepole pine forest, and rocky mountain subalpine mesic-wet spruce-fir forest and woodland). We also included an alpine stratum as a potential location for WBP colonization, and a shrub-herb stratum for forest that underwent a mortality event. The pathway diagrams for each forest stratum were then updated based on available information in the literature on species life histories.

We initialized the model with multiple sources of spatial data. The extent of the spruce-fir and lodgepole forest types were based on spatial data from LANDFIRE, and WBP was based on a previously published dataset that was generated with Landsat ETM+ imagery [51]. Cells that were classified as WBP by the LANDFIRE biophysical settings data, but not by the Landenburger et al. dataset were reclassified as spruce-fir because there was greater overlap in its BPS description than lodgepole. The location of state classes within each biophysical setting were spatially randomized in proportion to the LANDFIRE mapped state class distributions.

The entire simulated landscape covers just over 15 million hectares centered on the GYE, at a spatial resolution approximately 64 hectares per cell. The simulation was initiated in the year 1920 and run for 30 years to allow the model (specifically, tree populations) to stabilize. The historical period covers 1950-2010, and the model was then projected to 2100. The model included transitions for reproduction, aging, mortality, and competition. Modeled disturbances included replacement fire, mixed fire, disease, insect pathogens, and wind/weather. Parameters for these and other probabilistic transitions were derived from a variety of published sources (see Supplementary Material). Deterministic parameters (age-based transitions) were based on LANDFIRE pathway diagrams and life histories of the three main tree species. Fire size, seed dispersal distances, and mountain pine beetle outbreak size distributions were estimated based on the Monitoring Trends in Burn Severity (MTBS) dataset, U.S. Forest Service species information, and published literature [52].

The relationship between probabilities of mountain pine beetle infection and temperature was estimated using a modified logistic function:

where L is the maximum probability, μ is the 2-year moving average annual minimum temperature anomaly, Μ is the mean annual minimum temperature anomaly for 1999‒2009, and a is adjusted so the resulting probability matches observed beetle mortality for the same period. In order to calculate annual temporal multipliers, this probability was then divided by the historical rates of insect and disease infection from the LANDFIRE model (0.003). This relationship was used for both historical model runs (1950‒2010) and model projections (2010‒2100) in order to vary the probability of mountain pine beetle through time. Other probabilistic transitions were held constant through time.

Unfortunately, data was insufficient to parameterize blister rust in detail. There is a wide range of estimates of blister rust infection rates across locations within the GYE [31], and little information on the size of blister rust outbreaks. Moreover, the mortality of WBP with blister rust is often associated with other disturbances (mountain pine beetle, fire). Modeling blister rust is further complicated by the complex relationship between blister rust and climate (temperature and humidity both influence the spread of blister rust), and the fungus’ dependence on multiple hosts. As a result, we did not explore the influence of blister rust on WBP beyond including it in the model and ensuring that the modeled prevalence of infection and mortality matched observations from the GYE.

The establishment of new forest stands after disturbance was dependent on dispersal distances (from the literature), and on habitat suitability as defined by SDM output data. SDM output for each period (1950-1980, 1980-2010, 2010-2040, 2040-2070, 2070-2100) and each forest type (WBP, spruce-fir, lodgepole pine) were used as spatial multipliers to dictate probabilities of forest establishment (Figure 4). If establishment failed, the cell returned to the open state in the shrub-herb stratum.

Figure 4. Simplified pathway diagram of the integration of species distribution model output into the state-and-transition simulation model, including dispersal (blue arrow), depletion of the seedbank or mortality of juvenile trees (black arrow), successful establishment (whose probability is characterized by the SDM output, red arrow), and aging (green arrow). Other transitions (e.g., MPB infection, blister rust infection, fire) were omitted from the figure due to space considerations and because each forest stratum had a unique pathway diagram.

DownLoad: Full-Size ImgPowerPoint

We verified realistic model function and then evaluated model outputs to ensure they corresponded to real-world patterns (similar to techniques used in agent-based modeling [53,54]). In particular, verification involved model debugging and parameter calibration, and validation included sensitivity analysis and pattern evaluation. We calibrated the equation that dictated the relationship between mountain pine beetle and minimum annual temperature anomaly (Equation 1) to match observations (a similar approach to Frid et al. [55,56]). In particular, we compared observed levels of beetle and blister rust infection in GYE to results from simulations using four logistic curves with different maxima for the probability of beetle infection.

Evaluation of simulation models remains a significant challenge [57,58,59]. Our evaluation effort was limited by the availability of independent datasets, but we implemented sensitivity analysis and pattern evaluation [60] to check the validity of our model. We sought to reproduce the recent large-scale dynamics of WBP population and disturbance. We compared mountain pine beetle induced WBP mortality using different stratum shift (dead forest stands to open shrub-herb) probabilities to observations from the U.S. Forest Service aerial detection survey. Modeled results were subset by the area flown in each year of the aerial survey.

To test the robustness of our model to parameter uncertainty, we performed sensitivity analysis by running the model for the historical period with different values for parameters that might change substantially with climate or had lower certainty (were unknown or estimated, see Supplementary Material). In particular, we varied the probabilities of alpine colonization by WBP, replacement fire, and spruce-fir replacement of WBP stands by +/−50% [54,61]. We ran 40 Monte Carlo simulations (iterations or model runs) for each perturbation, and compared the means of our primary outcome variables (WBP population, and mountain pine beetle mortality and infection at year 2010) to those from the original model specification. The sensitivity index [61] was calculated as:

where p is the value of the independent variable, dp is the value for a change of p, x is the value of the dependent variable, and dx is the corresponding change in x in response to the change in p.

After model verification and validation we ran the STSM using a no climate change scenario and multiple scenarios that included combinations of two GCMs (CNRM-CM5 and HadGEM2-AO) and two RCPs (4.5 and 8.5) to capture the range of change in projected WBP habitat suitability. We compared temporal patterns and values of key outcome variables at 2100 for these four climate scenarios and a no climate change scenario. The no climate change scenario used SDM output from 1950‒1980 and the mean probability of mountain pine beetle infection for the same period.

4. Results

4.1. Verification and validation

Calibration of the relationship between annual minimum temperature anomaly and mountain pine beetle infection indicated that 0.15 was a reasonable estimate of the maximum probability of beetle infection. This value produced estimates of beetle infection and mortality, and blister rust infection that were comparable to observations from the literature (Figure 5); about half of sub-watersheds in the GYE had high to complete mortality from mountain pine beetle [62], and the proportion of infected WBP trees in 2010 was around 20% [63].

Figure 5. Results of model calibration using four different maxima for the probability of mountain pine beetle (MPB) infection (bars) as compared to observations (points). Columns represent the means of 40 Monte Carlo simulations and error bars the standard deviations.

DownLoad: Full-Size ImgPowerPoint

We also tested the sensitivity of key model outputs (total WBP, live mature WBP, beetle kill WBP, and beetle infected WBP) to changes in parameters with less certainty (alpine colonization, replacement fire, and spruce-fir competition) (Table 1). Perturbations of all of the three parameters produced significant differences in one of the three main model outputs. Reducing (by 50%) the time required for spruce-fir forest to replace WBP produced the largest sensitivity index values.

Table 1. Sensitivity analysis results

Outcome values are means of 40 model runs measured at year 2010. These values are compared to baseline model output using t tests; significant differences (p < 0.05 for the two-tailed distribution) are shown in bold, and sensitivity values are reported for statistically significant differences. Percent changes from the baseline model are in parentheses. (a)Alpine colonization could not be decreased further or increased by a smaller increment due to the limitations of the software (i.e., precision was limited to four decimal places) (b)The probability of replacement fire varied by state class; default and perturbation values for replacement fire are listed as ranges due to space considerations
Parameter	Default Value	Perturbation Value	Total WBP		Live Mature WBP		Beetle Kill WBP		Beetle Infected WBP
			Hectares (% Change)	Sensitivity	Hectares (% Change)	Sensitivity	Hectares (% Change)	Sensitivity	Hectares (% Change)	Sensitivity
Alpine colonization (probability)(a)	0.0001	0.0002	561, 060 (0.30)	0.30%	163, 524 (1.14)	1.14%	146, 815 (−0.05)	N/A	6, 838 (−0.68)	N/A
Replacement fire (probability)(b)	0.001-0.003	0.0015-0.0045	586, 759 (4.90)	9.79%	161, 757 (0.05)	N/A	142, 137 (−3.24)	−6.48%	6, 674 (−3.06)	N/A
		0.0005-0.0015	531, 778 (−4.93)	9.86%	162, 647 (0.60)	−1.19%	152, 048 (3.51)	−7.02%	6, 923 (0.56)	N/A
Spruce-Fir Competition (years)	350	175	459, 347 (−17.88)	35.76%	123, 146 (−23.83)	47.67%	103, 335 (−29.65)	59.31%	4, 818 (−30.01)	60.03%
		525	567, 269 (1.41)	2.83%	164, 816 (1.94)	3.88%	150, 577 (2.51)	5.01%	6, 537 (−5.05)	−10.10%

 | Show Table

DownLoad: CSV

In terms of pattern evaluation, the relationship between the probability of beetle infection and annual minimum temperature anomaly for the historical period were similar to simulated beetle survival published elsewhere [64, Figure 2A in 65]. Additionally, the final model specification (value of 0.1 for the probability of a stratum shift from dead forest to open shrub-herb) closely matched temporal patterns of mountain pine beetle induced WBP mortality from the aerial detections survey data, especially when compared to model specifications using other stratum shift values (0.05 and 0.0333) and a no climate change scenario (Figure 6). Model output for the total area of WBP in National Forest and National Park Service lands in 2010 (536,664 hectares) was also comparable to the observed area of WBP dominant stands (531,999 hectares) [31].

Figure 6. Comparison of mountain pine beetle mortality for four different model specifications (three different probabilities of the shift from dead forest to open shrub-herb, and a no climate change scenario).

DownLoad: Full-Size ImgPowerPoint

4.2. Projections

The area of WBP in 2100 was lower for all climate change projections as compared to the no climate change scenario (Figure 7), although this difference does not become pronounced until mid- to late-century (Figure 8). The no climate change scenario resulted in less spruce-fir and more lodgepole as compared to the GCM/RCP scenarios.

Figure 7. Mean areas of three forest types in 2100 under five different climate scenarios. Error bars represent the 95% percentile range for 40 Monte Carlo simulations.

DownLoad: Full-Size ImgPowerPoint

Figure 8. Projected area of WBP forest under five different climate scenarios. Lines represent the means of 40 Monte Carlo simulations and shaded regions represent the 95% percentile range.

DownLoad: Full-Size ImgPowerPoint

Despite the appearance of diminishing mountain pine beetle mortality over time and projections of low mountain pine beetle mortality in 2100 (Figure 9), beetle infected and killed WBP stands accounted for about 40% of susceptible WBP in all climate scenarios at the end of the century. Moreover, by 2100 late seral stage WBP was nearly absent from the landscape in all climate projections except for the no climate change scenario.

Figure 9. Projected area of WBP mortality from mountain pine beetle under five different climate scenarios. Lines represent the means of 40 Monte Carlo simulations and shaded regions represent the 95% percentile range.

DownLoad: Full-Size ImgPowerPoint

5. Discussion

Results of the sensitivity analysis suggest that some of the less certain model parameters had a substantial effect on model output. The time for WBP stands to transition to spruce-fir forest had the largest effect. More detailed information on this dynamic would likely improve model reliability. Estimates for seral (as opposed to climax) WBP communities for this parameter range from 100‒200 years [66]. Although the 350-year replacement time used in the model produced relatively conservative estimates of WBP persistence, a more nuanced approach could link spruce-fir competition to other factors (e.g., elevation). It is possible that this and other influential parameters, such as the rate of WBP colonization of the alpine zone and variation in the probability of stand replacing fires with climate, could be estimated based on fire, reforestation, and climatic patterns that have been reconstructed from fossil records [30,67].

Another limitation of this modeling effort was the use of output from SDMs of adult trees as input for seedling establishment. In other words, the spatial multiplier files used in our model represent the probability of the presence of adult trees, but dictated the probability of seedling establishment. However, these probability surfaces were based on 30-year moving averages of climate variables, so, depending on the age of the adult trees, these data may indeed represent the climate conditions during establishment. Moreover, our comparison of the presence and absence locations of mature and seedling WBP suggests that the distribution of seedlings is very similar to that of adults. In sum, without additional data points for the presence and absence of seedlings to parameterize more robust seedling SDMs, and given the overlap in WBP adult and seedling distributions it seems reasonable to have used the available SDM output for seedling establishment within the STSM.

Despite these limitations, the model produced results that closely matched historical point estimates and temporal patterns of WBP and mountain pine beetle activity in the GYE, indicating realistic model function. The counterintuitive finding of less spruce-fir under no climate change is likely the result of mountain pine beetle activity; the no climate change scenario had lower beetle infection probabilities, less beetle mortality, and thus fewer opportunities for spruce-fir to colonize cells vacated by beetle-killed lodgepole and WBP. Future work could incorporate SDM output of historic and projected mountain pine beetle habitat suitability in order to provide more nuanced and spatially explicit infection probabilities.

Output for key WBP variables were comparable across GCMs and RCPs. All scenarios showed a substantial decline in WBP forest, and especially late seral stage forest. The differential effects of climate projections may have been attenuated by similar levels of mountain pine beetle activity across climate scenarios. Across all climate projections, the probability of mountain pine beetle infection reaches a maximum by early- to mid-century.

Overall, this modeling effort integrated a wide range of available information on WBP ecology, reproduced observed patterns of WBP mortality in the GYE, and set the stage for future exploration of WBP management scenarios. The sensitivity analyses and findings also point to the importance of accounting for climate, disturbances, biotic interactions, and habitat suitability together. The approach outlined here for combining the SDMs and STSMs provides a strong foundation for exploring alternative approaches to managing resources under the combined pressures of climate change, insects, disease, and interspecific competition.

6. Conclusion

This research demonstrated the benefits of integrating correlative and stochastic simulation models. In particular, we combined SDM’s ability to produce statistically robust estimates of the relationship between climate and habitat with STSM’s ability to account for disturbances and biotic interactions. The resulting model reproduced observed patterns of WBP in the GYE.

Another important outcome from this research is the identification of data gaps and opportunities for the integration of additional datasets and modeling approaches. Our framework for model validation not only served to corroborate model function, but also identified important research needs, including the ability of WBP to colonize alpine zones, the rate of spruce-fir replacement of WBP across different parts of the landscape, and changes in fire regime. Looking to the historical and paleontological record could help address these shortcomings.

This model and the software platform offer a basis for exploring resource management scenarios such as thinning to protect high value trees from fire, prescribed fires, and the application of pesticides and pheromones to deter mountain pine beetles [31]. Ideally, the specific management scenarios would be carefully developed in conjunction with resource managers. Scenarios could then be implemented in ST-Sim by altering transition probabilities to further explore parameter space and model sensitivity, setting management targets for a given outcome (e.g., area of beetle mortality), or setting limits on expenditures for a given action (e.g., pesticide application). Moreover, these scenarios could be applied to specific management units (e.g., National Forests, National Parks), locations (e.g., < 10 km from roadways), and time periods. The costs of implementing management actions or achieving particular targets could also be tracked through time.

STSMs can integrate a variety of data types and sources to create spatially explicit, flexible, and verifiable representations of ecological dynamics. When paired with SDMs, they offer an especially powerful approach for anticipating the effects of climate change, and ultimately, exploring options for managing species in the face of an uncertain future.

Acknowledgments

We are grateful to Rick Lawrence for sharing data that were used to help initialize the STSM. We also thank Catherine Jarnevich and an anonymous reviewer for their helpful comments on previous drafts of this paper. This research was funded by the U.S. Geological Survey’s North Central Climate Science Center. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflict of Interest

All authors declare no conflicts of interest in this paper.

Appendix S1. Probabilistic transition parameter values, sources, and associated transition constraints

From Veg. Type	From Class	To Veg. Type	To Class	Transition Type	Probability	Age (years)	Source & Notes	Other Transition Constraints
All Forest	Propagule	All Forest	Juvenile	Establishment	From SDM	≥ 1	SDM output (habitat suitability grids)
All forest	Dead-Fire	Shrub-Herb	Open	Succession	Deterministic	≥ 2	Allows for seed dispersal to open cells after fire
Alpine	Open	WBP	Propagule	Colonization	0.0001	≥ 1	Estimated/unknown; tested with sensitivity analysis
Lodgepole	Late-Close	Spruce-Fir	Mature	Competition	Deterministic	350	LANDFIRE BPS 2110550 [68]	Min. time-since-transition (TST) (establishment) = 350 yr.
Lodgepole	Mid-Close	Lodgepole	Mid-Close	Competition/Maintenance	0.0020	20‒79
Lodgepole	MPB	Lodgepole	Dead- MPB	Death-MPB	Deterministic	≥ 20	Probability of death incorporated into probability of infection (as in LANDFIRE)	Min. TST (MPB infection) = 2 yr.
Lodgepole	Propagule	Shrub-Herb	Open	Establishment Failure	0.9500	≥ 1	Probability assumed equal for all forest types	Min. TST (dispersal) = 11 yr. (estimated based on MPB stands, Teste et al. [69])
Lodgepole	Mid-Close	Lodgepole	MPB	Infection-MPB	0.0030	20‒79	LANDFIRE BPS 2110550 [68]	Annual temporal multipliers vary this parameter according to Equation 1; outbreak size distribution estimated based on Aukema et al. [52]
Lodgepole	Late-Close	Lodgepole	MPB	Infection-MPB	0.0030	80‒350
Lodgepole	Late-Close	Lodgepole	Mid-Close	Insect/Disease- Other	0.0060	80‒350
Lodgepole	All	Lodgepole	Dead-Fire	Replacement Fire	0.0030	≥ 1	LANDFIRE BPS 2110550 [68]; tested with sensitivity analysis	Fire size distribution (MTBS [70])
Lodgepole	Dead-MPB	Shrub-Herb	Open	Succession	0.1000	> 20	Based on time for shrub-herb to early forest (10‒30yrs, Keane et al. [66]); used to track amount of dead forest	Min. TST (MPB mortality) = 2 yr. (loss of needles)
Shrub-Herb	Open	All Forest	Propagule	Succession	0.2500	≥ 1	Probability assumed equal for all forest types	Seed dispersal distances (USFS [71])
Spruce-Fir	Propagule	Shrub-Herb	Open	Establishment Failure	0.9500	≥ 1		Min. TST (dispersal) = 3 yr. (Johnson & Fryer [72])
Spruce-Fir	Late-Close	Spruce-Fir	Late-Open	Insect/Disease- Other	0.0010	≥ 150	LANDFIRE BPS 2110560 [68]
Spruce-Fir	Late-Close	Spruce-Fir	Late-Open	Mixed Fire	0.0050	≥ 150
Spruce-Fir	Early	Spruce-Fir	Dead-Fire	Replacement Fire	0.0020	1‒39	LANDFIRE BPS 2110560 [68]; tested with sensitivity analysis	Fire size distribution (MTBS [70])
Spruce-Fir	Mid-Close, -Open	Spruce-Fir	Dead-Fire	Replacement Fire	0.0020	40‒149
Spruce-Fir	Late-Open	Spruce-Fir	Dead-Fire	Replacement Fire	0.0025	≥ 150
Spruce-Fir	Early	Spruce-Fir	Mid-Open	Succession	0.0010	1‒39	LANDFIRE BPS 2110560 [68]
Spruce-Fir	Late-Open	Spruce-Fir	Late-Close	Succession	0.0010	≥ 150
Spruce-Fir	Mid-Close	Spruce-Fir	Mid-Open	Wind/Weather/ Stress	0.0010	40‒149
WBP	Late-Close	Spruce-Fir	Mid-Close	Competition	Deterministic	350	Based on LANDFIRE BPS 2110550 [68]; tested with sensitivity analysis	Min. TST (establishment) = 350 yr.
WBP	Mid-Close	WBP	Mid-Close	Competition/ Maintenance	0.0020	50‒129	LANDFIRE BPS 2110460 [68]
WBP	Early	WBP	Early	Competition/ Maintenance	0.0050	1‒49
WBP	MPB, Rust & MPB	WBP	Dead- MPB	Death-MPB	Deterministic	≥ 51	Probability of death combined with probability of infection (as in LANDFIRE)	Min. TST (MPB infection) = 2 yr.
WBP	Rust-Early	WBP	Dead-Rust	Death-Rust	0.5213	≥ 1	Keane et al. [66]	TST (rust infection) = 10‒30yrs (Hatala and Crabtree [73])
WBP	Rust-Mature	WBP	Dead-Rust	Death-Rust	0.5213	≥ 50
WBP	Propagule	Shrub-Herb	Open	Establishment Failure	0.9500	≥ 1	Probability assumed equal for all forest types	Min. TST (dispersal) = 4 yr. (USFS [71], Keane et al. [66])
WBP	Mid-Close	WBP	MPB	Infection-MPB	0.0030	50‒129	LANDFIRE BPS 2110460 [68]	Annual temporal multipliers vary this parameter according to Equation 1; outbreak size distribution estimated based on Aukema et al. [52]
WBP	Late-Close	WBP	MPB	Infection-MPB	0.0030	130‒350
WBP	Mid-Open	WBP	MPB	Infection-MPB	0.0020	50‒129
WBP	Late-Open	WBP	MPB	Infection-MPB	0.0020	≥ 130
WBP	Mid-Open, -Close	WBP	Rust- Mature	Infection-Rust	0.0067	50‒129	Logan et al. [65]; also similar to Keane et al. [66]	Initiated in 1970; outbreak size distribution (estimated)
WBP	Late-Open, -Close	WBP	Rust- LANDFIRE Mature	Infection-Rust	0.0067	≥ 130
WBP	Early	WBP	Rust-Early	Infection-Rust	0.0184	1‒49	Keane et al. [66]
WBP	Late-Close	WBP	Late-Open	Insect/Disease- Other	0.0020	130‒350	LANDFIRE BPS 2110460 [68]
WBP	Mid-Close	WBP	Mid-Open	Insect/Disease- Other	0.0030	50‒129
WBP	Late-Close	WBP	Late-Open	Mixed Fire	0.0020	130‒350
WBP	Mid-Open	WBP	Mid-Open	Mixed Fire	0.0070	50‒129
WBP	Late-Open	WBP	Late-Open	Mixed Fire	0.0070	≥ 130
WBP	Mid-Close	WBP	Mid-Open	Mixed Fire	0.0040	50‒129
WBP	Dead-MPB	WBP	Dead-Fire	Replacement Fire	0.0025	≥ 52	Estimated based on fire probabilities from BPS 2110460 (mixed evidence for influence of beetle/rust kill on fire); tested with sensitivity analysis	Fire size distribution (MTBS [70])
WBP	Dead-Rust	WBP	Dead-Fire	Replacement Fire	0.0025	≥ 2		Fire size distribution (MTBS [70])
WBP	Rust-Early	WBP	Dead-Fire	Replacement Fire	0.0025	1‒49
WBP	Rust-Mature	WBP	Dead-Fire	Replacement Fire	0.0025	≥ 50
WBP	MPB, Rust & MPB	WBP	Dead-Fire	Replacement Fire	0.0025	≥ 51
WBP	Early	WBP	Dead-Fire	Replacement Fire	0.0010	1‒49	LANDFIRE BPS 2110460 [68]; tested with sensitivity analysis
WBP	Mid-Close	WBP	Dead-Fire	Replacement Fire	0.0020	50‒129
WBP	Late-Close	WBP	Dead-Fire	Replacement Fire	0.0020	130‒350
WBP	Mid-Open	WBP	Dead-Fire	Replacement Fire	0.0030	50‒129
WBP	Late-Open	WBP	Dead-Fire	Replacement Fire	0.0030	≥ 130
WBP	Early	WBP	Mid-Open	Succession	0.0100	1‒49	LANDFIRE BPS 2110460 [68]
WBP	Mid-Open	WBP	Mid-Close	Succession	0.0500	50‒129
WBP	Late-Open	WBP	Late-Close	Succession	0.0500	≥ 130
WBP	Dead-MPB	Shrub-Herb	Open	Succession	0.1000	≥ 52	Based on time to transition from shrub-herb to early forest (10‒30yrs, Keane et al. [66]) in order to track amount of dead forest	Min. TST (MPB mortality) = 2 yr. (loss of needles)
WBP	Dead-Rust	Shrub-Herb	Open	Succession	0.1000	≥ 2

 | Show Table

DownLoad: CSV