Developing an expert elicited simulation model to evaluate invasive species and fire management alternatives

Abstract

Invasive species can alter ecosystem properties and cause state shifts in landscapes. Resource managers charged with maintaining landscapes require tools to understand implications of alternative actions (or inactions) on landscape structure and function. Simulation models can serve as a virtual laboratory to explore these alternatives and their potential impacts on a landscape. To be useful, however, managers need to participate in model development to ensure that model structure can evaluate the response of key resources to plausible actions. Here, we detail development of a state‐and‐transition simulation model (STSM) to evaluate buffelgrass (Cenchrus ciliaris L. syn Pennisetum ciliare (L.) Link) in Saguaro National Park (SNP), Arizona, USA, through collaboration between managers and researchers. We integrate expert knowledge and research to create and parameterize a stochastic, spatially explicit STSM to evaluate specific management objectives. We also develop a dynamic link between the STSM and a fire behavior model to allow exploration of potential novel processes introduced to the ecosystem by buffelgrass invasion. Our projections show that buffelgrass can be expected to increase on the landscape and that the integration of fire into the model accelerates the projected rate of invasion and increases degradation of resources of management concern. We highlight the benefit of engaging end users in the modeling process so that the model is targeted to evaluate management objectives, in this case retention of saguaro cacti (Carnegiea gigantea (Engelm.) Britton & Rose) on the landscape. Being able to integrate an external model that can help address the unique characteristics of a problem such as the introduction of fire into the SNP desert ecosystem increases the ability of simulations to address complex ecological and management questions.

Introduction

Landscapes are dynamic, influenced by different processes interacting over space and through time, resulting in a heterogeneous pattern of landscape states. State and transition models (STM) characterize different states of landscape vegetation and the processes, or transitions, which occur to move a landscape between states. These models have a long history in the field of range management (Briske et al. ) and have typically been expressed as conceptual diagrams, describing hypotheses of vegetation and landscape change that can be used to inform management actions intended to move a landscape toward a desired state (or prevent a landscape from moving toward an undesired state). Simulation models are another class of models that can be used to forecast or hindcast landscape change through time. The state‐and‐transition simulation model (STSM) approach integrates STM conceptual models into a simulation model framework by assigning parameter values to model transitions and assigning each cell of a grid placed over the landscape to an initial state, allowing a quantitative implementation that can be used to project landscape conditions through time (Daniel et al. ). State‐and‐transition simulation models have captured many processes affecting landscapes including forest management (Daniel et al. ), invasive species dynamics (Jarnevich et al. ), carbon flows (Daniel et al. ), and rangeland management (Provencher et al. ).

Parameterizing an STSM can be challenging because it requires the assembly of detailed landscape data and the identification of parameter estimates for many uncertain variables. For some landscapes, historical information exists to derive a distribution of values for a specific parameter, such as fire frequency (Daniel et al. , Miller et al. ), but more often empirical data are lacking. In these cases, expert elicitation can be combined with available empirical data to parameterize models (Frid et al. , Miller et al. ). Global change dynamics, including climate change and biological invasions, often result in no‐analog communities for which historic information does not exist (Williams and Jackson ). When STSM representations are too simplistic or rely on historical information that is unavailable or that may no longer represent the system, it may be beneficial to couple an STSM with other types of simulation models. For example, some aspects of fire regimes can be simulated within an STSM based on historic information about fire cycle and fire sizes (Daniel et al. ). However, integrating a STSM with a fire model that is based on physical characteristics rather than historical information could alleviate issues arising from no‐analog conditions. Here, we specifically focus on the development and parameterization of a STSM linked to a fire growth and behavior model (FARSITE, Finney ). This study represents the first case of the STSM approach allowing for a dynamic linkage with other existing external models.

The STSM for this study is developed for Saguaro National Park (SNP) to simulate the invasion of buffelgrass (Cenchrus ciliaris L. syn Pennisetum ciliare (L.) Link), an invasive African bunchgrass that poses a novel fire risk to the Sonoran Desert. Saguaro National Park is located near Tucson, Arizona, USA, and was designated to preserve Sonoran Desert fauna and flora, including the giant saguaro cacti (Carnegiea gigantea (Engelm.) Britton & Rose) which are found only in a small portion of the United States. Saguaro National Park is comprised of two separate units to the east and west of Tucson. The 9726 ha Tucson Mountain District (TMD) on the west has elevations ranging from 664 to 1429 m and is almost entirely a desert ecosystem. The 27,279 ha Rincon Mountain District (RMD) on the east has greater topographic variation, with elevations between 814 and 2641 m, and contains a desert ecosystem, a desert grassland transition zone, and a higher elevation conifer forest. Portions of the desert ecosystem are susceptible to invasion by buffelgrass, which was first recorded in SNP in 1989 and has invaded an estimated 1000 ha despite an active management program across both units. A primary management concern is the potential for the formation of large patches of buffelgrass which have a much higher fuel load than native vegetation (McDonald and McPherson ). Although fires regularly occur at higher elevations in the conifer forests of SNP, the lower elevation native desert vegetation is not fire‐adapted, and there is concern that fires associated with large invaded patches could result in transformation to a system like an African savanna with the fire‐adapted buffelgrass taking over (Esque et al. ). Even with the absence of fire, buffelgrass can still transform the ecosystem through its presence and spread (Olsson et al. ). Our goal was to build an STSM to address SNP's and other nearby public land managers’ most pressing questions about the future of buffelgrass across the landscape. We implement the model to characterize the potential magnitude of the buffelgrass problem in SNP and build a tool that will enable SNP managers to evaluate alternative management actions and quantify the importance of ecological uncertainties.

Methods

Model overview

The process to develop the SNP buffelgrass STSM included the following: eliciting expert knowledge about the system from local managers and researchers via an expert workshop; obtaining spatial information about buffelgrass habitat suitability (Jarnevich et al. ); and gathering existing data on buffelgrass locations and spread, local climate, and LANDFIRE fuel loads (LANDFIRE ; Fig. ). These data were input into a dynamically linked landscape simulation model and fire behavior model. The landscape simulation model was built using the SyncroSim (version 2.0.21; www.syncrosim.com) software framework with the stsim‐farsite package (version 3.1.21) module (Daniel et al. ) linked to the fire growth and behavior model FARSITE (version 4.1.055; Finney ).

Saguaro National Park buffelgrass model development process. Orange boxes show model components, blue and purple boxes describe model input data, and green shows main model outputs.

Expert workshop

We convened a two‐day workshop of regional buffelgrass experts at SNP in January 2016, to solicit expert knowledge and guide the development and parameterization of the buffelgrass STSM. Workshop attendees included representatives from SNP, members of the National Park Service Exotic Plant Management Team, the Southern Arizona Buffelgrass Coordination Center, U.S. Forest Service, Tohono O'odham Nation, Sky Island Alliance, the University of Arizona, Pima County, and the U.S. Geological Survey. The experts were selected to include representatives from organizations working on buffelgrass management in the Tucson area, targeting individuals with field experience, and local researchers studying buffelgrass.

Problem framing

The workshop included facilitated discussions and exercises to identify primary buffelgrass management objectives and management scenarios to be addressed by the STSM. Workshop participants were divided into four groups, each with a facilitator, to discuss pre‐defined questions that could help with identification of management objectives, scenarios, and model structure (Appendix S1). The participants identified vegetation conversion, fire, and the loss of biodiversity and native species, especially the loss of the Park's iconic saguaro cacti, as primary threats from buffelgrass invasion. Identified secondary effects of vegetation conversion and fire included reductions in tourism, damage to historic and cultural resources, and loss of wildlife habitat. These threats and concerns guided the identification of three primary management objectives which provided sideboards for model development: minimizing loss of biodiversity, protecting and preserving saguaro cacti, and minimizing treatment costs. To measure treatment costs, the STSM would need to track management transitions made for each simulation scenario and pair these data with historic SNP buffelgrass treatment and inventory cost data. To measure biodiversity loss and saguaro cacti protection, participants would need to ensure that transitions in the conceptual STM could capture the major drivers of biodiversity and saguaro cacti loss and that states in the conceptual STM could provide proxy measures for biodiversity loss and saguaro loss. These considerations are discussed further in the conceptual STM section below.

The initial breakout discussions also identified management and ecological uncertainties that could be explored via simulations. Participants identified rainfall variability and the interaction of rainfall with the rate of patch cover increase as important ecological uncertainties. Another key source of ecological uncertainty was fire. Participants were specifically interested in exploring how different management strategies could affect the ability for the invaded desert ecosystem to carry and spread fire and in the possible transfer of fire between the lower and upper elevation ecosystems. They were also interested in how the frequency and severity of fires may change with buffelgrass invasion over time. We discuss how these uncertainties were incorporated into the conceptual model below.

Conceptual STM

Workshop participants modified a previous conceptual model developed by Jarnevich et al. (; Fig. A) to ensure that model outputs could provide proxy measures for identified management objectives, evaluate identified scenarios, and incorporate management concerns related to uncertainty in buffelgrass spread with fire and precipitation interactions (Fig. B). As with the previous conceptual model, each invaded state on the landscape is described as either detected or undetected. Transitions included natural transitions (buffelgrass patch cover increase, mortality, and fire) and management transitions (inventory and treatment). Inventory transitions were probabilistic, and either successful (solid lines in Fig. B moving to detected states) or unsuccessful (failure to detect established buffelgrass remaining as undetected; dashed lines in Fig. B). Only detected locations can undergo management.

Conceptual state and transition model (A) from Jarnevich et al. (), used to begin discussions for the Saguaro National Park (SNP) conceptual model, and (B) SNP conceptual model developed by workshop participants. The boxes represent potential landscape states related to buffelgrass cover classes (left to right) and detection (top to bottom). Arrows represent transitions between the states including increased buffelgrass patch cover, natural mortality of buffelgrass, detection of previously undetected patches (inventory), and management actions (treatment/retreatment). Dotted lines represent unsuccessful transitions. The SNP model modifies bin cover levels to reflect fuel loads and treatment thresholds and adds both fire transitions and treated and converted states.

Buffelgrass cover class bins were modified for the SNP model to reflect fuel loads and treatment thresholds. Workshop participants selected three cover thresholds as important for fire spread and management: <1% buffelgrass cover to represent a maintenance level, 10% buffelgrass cover required to carry a fire (or 0.5 ton/ac of fuel, estimated based on McDonald and McPherson ), and 50% buffelgrass cover, which is SNP's threshold for aerial application of herbicide identified in SNP's Restoration Management Plan and Environmental Assessment document (National Park Service ). At >50% cover, SNP considers the biodiversity of the localized invasion to be compromised and is willing to risk collateral damage to native plants.

Fire transitions were added to the conceptual model to reflect expected shifts in states anticipated by fires occurring in differing densities of buffelgrass cover (Fig. B). Because 10% buffelgrass cover is the minimum amount to carry a fire, fires that spark on cells with less than 10% cover would not shift the state of the cell. Participants agreed that fires in cells with >10–50% buffelgrass cover would transition to the denser >50% cover level, but that some saguaros and native vegetation would likely still be maintained in the cell. Due to fuel loads associated with dense patches of buffelgrass and the fire intolerance of the native desert plant community, participants agreed that fires occurring in dense buffelgrass (cells with >50% cover) would result in a converted landscape with a loss of nearly all saguaro cacti and native vegetation. This converted state would be persistent due to the inability of native plants to establish in buffelgrass patches (Morales‐Romero and Molina‐Freaner ). There was a discussion about whether the converted state could be reached without a fire (i.e., through continued buffelgrass infill over time, Olsson et al. ); the prevailing opinion was that it could not be reached within the 30‐yr time horizon of the simulation model. Participants therefore agreed that the conceptual model needed to distinguish between patches densely populated by buffelgrass (which would affect the objective to minimize losses to biodiversity but would be less likely to affect the objective to preserve saguaro populations), and patches that had been burned and thereby converted to buffelgrass grassland (which would affect both the biodiversity objective and saguaro objective). Changes in the initial state of the landscape compared to future simulated states of the landscape could be used as proxy measures for losses in biodiversity and losses to saguaro cacti populations.

Based on observed responses to treatment actions, participants agreed that dense patches of buffelgrass require consecutive years of treatment to reduce cover levels. There was consensus that buffelgrass treatment and retreatment would be best modeled using time‐dependent deterministic transitions. Participates agreed that patches with >50% cover would need to be treated two consecutive years to transition to the >10–50% cover class state; otherwise, these cells would transition back to the >50% cover class state. Participants agreed to add this as a deterministic rule within the conceptual model, and two new states were added to enable the tracking of consecutive year treatments for >50% cover level cells.

The Jarnevich et al. () model (Fig. A) was developed for a drier location and contained a buffelgrass natural mortality transition that was based on limited evidence of potential buffelgrass dieback associated with drought conditions; the group agreed that this transition could be removed from the SNP model due to SNP having a relatively wetter precipitation regime. Thus, the only way to reduce buffelgrass cover class levels in the SNP model is through treatment.

Model parameterization

Transition parameters

To parameterize model transitions, we used a combination of values from the literature and expert elicitation (O'Hagan et al. , Oakley and O'Hagan , Martin et al. ). Prior to the expert workshop, we conducted a literature review of the physiological limitations and life‐history traits of buffelgrass (Jarnevich et al. ; Appendix S1) and collected available information on each parameter included in previous STSMs and from regional research (Frid et al. , Jarnevich et al. ). We also asked participants to be prepared to discuss buffelgrass dispersal, patch cover increase, patch expansion, seedbank dynamics, and treatment and detection effectiveness. During the workshop, we utilized the Sheffield Elicitation Framework (SHELF) quartile method for identifying median, upper, and lower quartile estimates for select transition parameters (Oakley and O'Hagan ). We began the elicitation exercise with a training/practice problem to orient participants to the elicitation process. Following the orientation, we presented definitions for each parameter in the context of the STSM (Table ) and provided information gleaned from the literature and previous STSMs. Experts considered the available parameter data and collectively decided whether the presented values from the literature or previous STSMs were appropriate for use in the SNP model or whether alternative values were needed. We followed SHELF guidelines for eliciting values for parameters that participants agreed required updating. The process for each parameter started with a group discussion to come to a consensus on the plausible range (minimum and maximum) of values for the parameter under consideration. Individuals were then asked to anonymously record their estimates of the median and quartiles for the parameter within the plausible range. The set of individual estimates were shared with the group, and the group collectively discussed the parameter distribution to come to consensus median and quartile values for each elicited parameter. We elicited transition parameters for 10 × 10 m cells (100 m²) but built the STSM using 0.25‐ha cells (2500 m²). Therefore, we needed to scale elicited transition rates to 0.25 ha (Appendix S3).

Elicited and estimated transition parameters for the Saguaro National Park buffelgrass model

Transition parameter	Definition	1st quartile (yr)	Median (yr)	3rd quartile (yr)
Patch cover increase rate: <1% to 1–10%	Amount of time it takes for infilling of a buffelgrass patch to move from initially invaded (e.g., <1% cover) to 1–10% cover	1	3	5
Patch cover increase rate: 1–10% to >10–50%	Amount of time it takes for infilling of a buffelgrass patch to move from 1–10% cover to >10–50%	2	6	11
Patch cover increase rate: >10–50% to >50%	Amount of time it takes for infilling of a buffelgrass patch to move from >10–50% cover to >50%	3	5	10
Patch cover increase rate: <1% to >50%	Amount of time it takes for infilling of a buffelgrass patch to move from initially invaded (e.g., <1%) to >50% cover	6	11	17
Wet year time‐to‐transition divisor	Time‐to‐transition divisor to decrease the number of years it takes to move from one cover level to the next during a wet year	1.9	2.7	3.9
Probability of a wet year	Probability that a year has a monsoon length and cumulative precipitation above the 30‐yr medians	10%	27%	57%

Notes

Patch cover increase rates and wet year time‐to‐transition divisors were elicited during the January 2016 Saguaro National Park expert workshop using the Sheffield Elicitation Framework quartile method for a 10 × 10 m area with average environmental conditions (i.e., not a wet year or a drought year). The frequency of wet years distribution is based on weather station data from the Tucson International Airport (GHCND:USW00023160).

While we were able to elicit information about how wet years affect the rate of buffelgrass infill, the model also required an estimate of the probability of a wet year (i.e., wet enough to decrease the amount of time for buffelgrass to move between cover classes). To estimate this probability, we obtained daily climate data for 1987–2016 from the nearest long term weather station, located at Tucson International Airport (GHCND:USW00023160), and calculated the length of and cumulative precipitation during the monsoon season (July–September) in R (R Core Team ) using the seas (Toews et al. ) and tidyverse packages (Wickham ). Eight of 30 yr had both the monsoon length and cumulative precipitation above the 30‐yr medians, meaning the probability of a wet year was 27%. We used this as the frequency of wet years to be sampled for each time step in the STSM. We also calculated the upper and lower quartiles of wet year frequency similar to quartiles for expert elicited values that could be used for sensitivity analyses in the future. Lower and upper quartiles for the probability of a wet year were 10% and 57% (17 and 3 of 30 yr), respectively.

Spatial parameters

We used a 0.25‐ha cell size for spatial analyses in the STSM. A habitat suitability layer (Jarnevich et al. ) defined where buffelgrass was able to establish, partitioning the desert ecosystem into susceptible and unsusceptible strata (Appendix S2: Fig. S1). The habitat suitability layer was discretized into these two strata using a probability threshold value of 0.5. The habitat suitability model had a 10 × 10 m resolution, and we aggregated values to 0.25‐ha cell using the presence of one or more susceptible cells at the finer resolution to classify a cell as susceptible. Buffelgrass seeds landing in raster cells with the unsusceptible designation will fail to establish, so the unsusceptible strata only included uninvaded and seedbank states.

Saguaro National Park's STSM model is focused on the desert ecosystem; however, we also wanted to be able to simulate how increases in buffelgrass cover could affect fire spread from the desert through the grassland transition to the forest. A grassland ecosystem stratum was added, represented by a single static state and, based on SNP staff observations, defined as north‐facing slopes with elevation >1371.6 m (4500 ft) or south‐facing slopes with elevation >1524 m (5000 ft). The forest ecosystem (defined as any location >1828.8 m; 6000 ft) was not represented in the model.

We combined the most recent aerial survey data from 2012 and the most recent ground surveys from 2014 to 2015 to define initial landscape states (Appendix S2: Fig. S2). Because of the mapping size of most patches, we used a 10 × 10 m grid to rasterize the mapping data. To initialize undetected buffelgrass states, we used inventory detection probabilities provided by SNP staff to predict the amount of undetected buffelgrass within each state and then randomly distributed undetected buffelgrass to susceptible locations on the landscape based on this detection rate. We then scaled the detected and undetected patches from 10 × 10 m to 0.25 ha (50 × 50 m) by summing the area covered by buffelgrass in the 10 × 10 m cells being aggregated and then converting each summed value to a percent cover of the 0.25‐ha pixel. We ran the model for one‐time step and used locations classified as seedbank within susceptible desert at the end of the time step as the initial seedbank states.

We simulated spread of buffelgrass from invaded cells to neighboring cells using a probability distribution of annual spread distances (i.e., dispersal kernel) with a Weibull distribution, similar to that used in previous models (Frid et al. , Jarnevich et al. ). While SNP has mapped buffelgrass data from many years, there is active management occurring throughout the landscape, and therefore, no time series of unmanaged buffelgrass patch sizes was available to calibrate the dispersal kernel. Thus, instead of using SNP data, we calibrated the dispersal kernel for the STSM using a time series of spread data from unmanaged land in the nearby Santa Catalina Mountains (Olsson et al. ). We varied the Weibull scale parameter (α) between 0.02 and 0.26 and the Weibull shape parameter (β) between 0.1 and 1.5 and then selected the parameter combination that resulted in a spread distribution most similar to the Santa Catalina Mountain time series data (α 0.14 and β 1.0).

Saguaro National Park only treats buffelgrass within Park boundaries, but buffelgrass dispersal can occur from external sources. To account for seed dispersal into SNP from outside the boundaries, we included a buffer zone of one mile around the SNP boundary, matching the buffer distance defined in SNP's Fire Management Plan (National Park Service ). We did not allow buffelgrass patches to grow or have treatment applied within the buffer as we do not have information on control efforts outside SNP boundaries. We used EDDMaps records of buffelgrass locations to initialize the locations within the buffer (EDDMapS ). The area within the buffer was excluded from our analyses of simulation outputs.

Integrating fire model

We integrated fire into the STSM using the FARSITE fire area simulator software version 4.1.055 (Finney ). FARSITE was integrated dynamically into the ST‐Sim simulations using the built‐in ability of the underlying SyncroSim framework to link models, and a customized FARSITE package to ST‐Sim. At each simulation time step prior to applying model transitions, ST‐Sim pauses and invokes the command line version of FARSITE via a python script. This script supplies FARSITE with the fuel model attributes of each cell in the landscape at that point in time. As vegetation changes in the landscape due to buffelgrass invasion, so does the fuel model for each cell. In addition to fuel models, the FARSITE model also required non‐spatial parameters including conditions related to fuel moisture levels, ignition frequency, fire season length, fire duration, wind speed and direction, topography, precipitation, and temperature. These values were selected by reviewing historical weather and fire conditions in SNP (National Park Service ; Tucson International Airport weather station data [GHCND:USW00023160], Appendix S2: Table S1). Additional FARSITE parameters identified in Appendix S2: Table S1 were calibrated to match environmental conditions in which historic fire occurred and on recommendations by SNP's fire ecologist.

Based on the fuel model for each cell, the number of ignitions, and weather information, FARSITE produces fire perimeters which are passed back to ST‐Sim in the form of fire probabilities for each cell. Simulation cells within the fire perimeter that are burnable have their fire probability set to one. All other cells have their fire probability set to zero.

To enable the link between ST‐Sim and FARSITE, we modified the states in the STSM to capture buffelgrass cover classes specific to fire and fuel loads. These classifications were based on LANDFIRE surface and canopy fuel type designations with adjustments from SNP staff based on their knowledge of existing conditions. Since RMD covered two LANDFIRE fuel model geo‐area zones classified by separate subject matter experts, adjustments were made to decrease the fuel load for the western zone to more accurately reflect existing conditions observed by SNP staff. Vegetation classes included Short, Sparse Dry Climate Grass (GR1); Moderate Load, Dry Climate Grass‐Shrub (GS2); Low Load, Dry Climate Shrub (SH1); and Nonburnable Bare Ground (Scott and Burgan ). The uninvaded state was broken into four different states corresponding to existing LANDFIRE fuel type designations, and these changes were made to both the initial conditions (combine Appendix S2: Fig. S2 with Appendix S2: Fig. S3) and the conceptual model diagram (Fig. B expanded to Appendix S2: Fig. S4). The seedbank state fuel load matched the underlying fuel vegetation type (Appendix S2: Fig. S3). Buffelgrass cover classes were assigned to a fuel type based on equations relating percent cover of buffelgrass to fuel loads (McDonald ). This equation resulted in <1% and 1–10% cover levels assigned to Low Load, Dry Climate Shrub (SH1), >10–50% cover levels assigned to Low Load, Dry Climate Grass (GR2), and >50% cover levels assigned to Moderate Load, Dry Climate Grass (GR4).

Simulations

We conducted analyses for the TMD and RMD areas separately due to their geographic separation and different vegetation communities. We ran the simulations without management actions to establish a baseline for comparison with future analyses. We ran a scenario with fire and one without for 30 time steps (from 2014 to 2044), each with 20 Monte Carlo iterations. We analyzed both the state of the landscape at each time step and the invaded area (calculated by multiplying the area of each cell by the midpoint of its current percent cover class). We used burned cell information at the grassland/forest border to determine the number of fires in the grassland that could move into the forest ecosystem to determine potential impact to that system.

Results

Model files and results are available as a U.S. Geological Survey data release (Jarnevich ). See Appendix S2 for details on accessing the software and model.

Elicited parameters

Elicited values for the time to transition from the lowest (<1%) to the highest (>50%) cover class varied depending on whether workshop participants were asked to provide the total number of years it would take for a patch to transition from the lowest (<1%) to the highest (>50%) cover class or whether the total was calculated from participants being asked for the time to transition between each individual cover class (i.e., time from <1% to 1–10%, from 1–10% to >10–50%, and from >10–50% to >50%; Table  elicited values and Appendix S3: Table S1 transformed to 0.25‐ha cells). The sum of the median individual transition estimates is larger than the elicited median total years to transition from lowest to highest cover (14 vs. 11 yr, respectively). The range in quartile estimates between the two approaches was similar (lower quartile of 6 yr for both; upper quartile of 16 yr for the total approach or 17 yr for the summed approach). The difference in median times to transition across approaches, and the interquartile range of 10 or 11 yr, reflected the overall high level of uncertainty that participants had about the rate of buffelgrass patch infilling.

Experts also estimated the effect of precipitation stochasticity on buffelgrass cover class transition rates by defining a time‐to‐transition divisor that decreases the number of years it takes to move from one cover level to the next during a wet year. For example, with average wetness it takes 3 yr to move from the <1% to the 1–10% cover level; the median wet year time‐to‐transition divisor of 2.7 would cut this transition time to 1 yr. Participants agreed that the wet year time‐to‐transition divisor would be consistent across the cover class states.

Simulations

Modeled suitable habitat (i.e., areas potentially susceptible to invasion) in SNP was primarily located in lower elevations and on south‐facing slopes (Appendix S2: Fig. S1). Initial conditions represented the state of the landscape in 2014, the beginning of the model simulation time period (Table ; Appendix S2: Fig. S2). Rincon Mountain District was more invaded than TMD at the start of the simulation (235 ha in RMD; 53.5 ha in TMD), despite having a similar amount of susceptible habitat (6070 ha in RMD; 4775 ha in TMD). Both the number of cells classified as invaded by buffelgrass and the derived invaded hectares increased through time in both units (Figs. , ) and was still increasing at a steady or accelerating rate at the end of the simulation (Fig. ). With no management, the projected buffelgrass invasion after 30 yr reached an average extent of 6842 ha (RMD: 4323 ha; TMD: 2519 ha) without fire and 8547 ha (RMD: 5147 ha; TMD: 3400 ha) with fire (Fig. ).

Landscape area in different states initially (2014) and at the simulation end (2044) with and without fire for the Rincon Mountain District (RMD) and the Tucson Mountain District (TMD) of Saguaro National Park

Category	RMD 2014	RMD 2044 fire	RMD 2044 no fire	TMD 2014	TMD 2044 fire	TMD 2044 no fire
Susceptible desert (ha)	6070	6070	6070	4775	4775	4775
Seedbank (ha)	0.5	140	246	0.7	176	202
<1% cover (ha)	63	544	817	116	502	625
1–10% cover (ha)	593	412	340	225	412	312
>10–50% cover (ha)	417	259	359	86	247	291
>50% cover (ha)	102	1286	2807	20	1428	1291
Converted (ha)	0	2646	0	0	811	0
Invaded area, ha (Percent susceptible desert invaded)	235 (4)	3717 (61)	2238 (37)	53.5 (1)	1982 (42)	1071 (22)

Note

Invaded area is calculated by multiplying the mean of the cover class bin by the area classified as the cover class bin.

Average area of landscape cells assigned to each state class across Monte Carlo realizations, representing buffelgrass cover classes, at the start and end of 30 yr of simulation without management for Rincon Mountain District (RMD) and for Tucson Mountain District (TMD) with and without fire.

Average area invaded by buffelgrass through time, where area invaded is calculated by multiplying the mean of the cover class bin by the area classified as the cover class bin, including scenarios for the Rincon Mountain District (RMD) and for Tucson Mountain District (TMD) both with and without fire. The vertical gray bars at the last time step (2044) represent the range in values across the 20 Monte Carlo iterations.

Including fire resulted in both more buffelgrass cover and more invaded locations on the landscape (Figs. , ). The impact of fire on the area of invaded hectares begins to have a noticeable impact after 10 yr into the simulation (Fig. ). At this point, the presence of fire on the landscape results in much more area invaded by buffelgrass. Furthermore, the simulations that included fire resulted in an average of 2733 ha in RMD and 811 ha in TMD in the converted class by the end of the simulation (Fig. ).

On average, 15,418 ha of the landscape burned at least once during the 30‐yr simulation (RMD: 12,260 ha, model runs ranged from 9276 to 14,622 ha; TMD: 3158 ha, model runs ranged from 2103 to 4428 ha). There was a large variability in the amount burned each year, and this variability increased through time (Fig. ). The cumulative burned area of the simulation, calculated by summing each annual burned area, also varied across iterations, ranging from 19,796 to 39,902 ha with many locations burning more than once in the 30‐yr period. There was an average of 8 (range of 5 to 11) fires that burned from the grassland to the study area border with the higher elevation forest of RMD during the 30‐yr simulation. In RMD, 66% of years had an ignition; in TMD, 76% of years had an ignition. While we parameterized the model to have a 70% chance of an ignition each year, there were additional years that resulted in no burned area when an ignition occurred in locations with fuel loads that did not support fire. These ignition years without burning included an average of 4 yr in RMD (range in realizations from 1 to 9 yr) and an average of 7.4 yr in TMD (range in realizations from 2 to 12 yr). Ignition years without fire plus years without an ignition resulted in an average of 13.8 yr without burning in RMD (model runs ranged from 8 to 18 yr) and an average of 15.7 yr without burning in TMD (model runs ranged from 10 to 21).

Amount of the landscape that burned during each year of the 30‐yr simulation, where different colors represent different potential realizations (i.e., Monte Carlo iterations) for Rincon Mountain District (RMD) and for Tucson Mountain District (TMD). The solid black line represents average area burned, and the dashed black line represents average area invaded (from Fig. ) by buffelgrass over the 30‐yr simulation across the different potential realizations.

Discussion

We provide a detailed description of the process of designing, building, and parameterizing an STSM of an invasive grass with a dynamic link to a fire behavior model. While we currently present only basic information on the scope of the buffelgrass problem, this STSM tool could be used to explore management questions and simulate scenarios to inform buffelgrass management in SNP and the greater Sonoran Desert region, with the goal of providing insights into the possible scope and damages of the invasion, the required resources for control, effectiveness of alternative management strategies, and analyses of ecological uncertainties (i.e., buffelgrass class cover transition rate and climatic uncertainty). Furthermore, results from current analyses can be summarized and used for public education and awareness surrounding the buffelgrass invasion issue in SNP, communicating scientific results to the public to improve knowledge sharing (Shanley and López ). For example, the spatial results can be combined to develop a time‐lapse video of the landscape over the simulation period for a succinct visualization of the potential scope of the invasion, which we did for a single Monte Carlo using buffelgrass state classes at each time step, starting with initial conditions in 2014 (Appendix S2: Fig. S2) and ending with 2044 (Video S1). While a single Monte Carlo represents one possible realization of the model, which will vary due to stochasticity of the parameters, it provides a visual tool for educating people about how buffelgrass may spread in the absence of management actions.

Building the STSM was an iterative and collaborative process. We began by using previously developed models and then iterated through multiple revisions during ongoing communication with SNP managers concerning initial results, modifications due to model limitations, and modifications required to design the model to address identified key management questions. Engaging local experts and stakeholders in framing the problem, developing the conceptual model, estimating model parameters, and devising simulation scenarios was an important component to developing a model that could be used as a tool for decision‐making (Addison et al. ). Including practitioners in the development resulted in a model tailored to management objectives, including adding fire with no‐analog conditions and matching cover classes to fire and management activities. Furthermore, involving practitioners in the process also helps to facilitate the communication of model results (Walsh et al. ). Saguaro National Park staff who will be using the model results have a clear understanding of the model components and limitations because of their involvement in model development. This model co‐production has also facilitated model implementation and modeler/manager communications because of SNP staff's working understanding of the model. The importance of this co‐production has been highlighted in development of other STSMs (Miller et al. ), although the specific methodology used for engagement differed.

One of the primary management objectives expressed by the stakeholders is to minimize loss of the saguaro cacti, a keystone species (Drezner ). Measuring this objective requires the integration of fire into the STSM due to cacti's susceptibility to fire events (Rogers ). Red brome (Bromus madritensis), another invasive grass, which generally produces fuel loads orders of magnitudes less than buffelgrass, was the main carrier of a moderate severity fire in SNP in 1994 (Esque et al. ). Twenty‐five percent of saguaros in the burned area were dead six years post‐fire. Another study of fire effects on saguaro found mortality still occurring ten years post‐fire and greater impacts with higher severity fires (Narog et al. ). The converted class in our simulations represents this loss of saguaro cacti as a result of fire, whereas saguaro could still survive in unburned high cover classes of buffelgrass.

Fire, as simulated by FARSITE, was an important component of simulations and consistently resulted in a greater projected extent of buffelgrass invasion in SNP. These results suggest that fire may play an important role in the spread of buffelgrass, highlighting a positive feedback between fire and invasion that has been found in other studies both in Arizona (Esque et al. , McDonald and McPherson ) and elsewhere (Butler and Fairfax , Clarke et al. , Miller et al. ). For historic context, during the last 50 yr there have been only two fires that burned from low elevation into high elevation according to SNP's fire history records (National Park Service ), while our simulations resulted in an average of eight in 30 yr. Historically large wildfires in the Sonoran Desert were extremely rare, and fire behavior was relatively mild. Lack of continuous fuels meant that most fires were small (<1 ha; National Park Service ). Thus, buffelgrass may alter fire dynamics beyond the desert ecosystem of SNP. Additionally, both fire size and frequency may be expected to increase in SNP, impacting the non‐fire adaptive environment and SNP's namesake, the saguaro cacti.

There are limitations to using STSMs to address land management questions. Invasion by buffelgrass produces a novel state, with no historic analog (e.g., introduction of higher fire potential; McDonald and McPherson ); therefore, these simulations can only provide forecasts amidst considerable uncertainty. Further, the STSM development process is complex and time‐intensive, and the number of parameters can be large and rapidly increase as additional complexity is added, which also compounds the uncertainty in the model results. We attempted to limit the level of complexity; for example, when simulating fire for each year of the simulations, we used one set of conditions to represent the weather, fuel moisture, and other fire‐related inputs. Complexity in the STSM is also increased by involving multiple stakeholders and experts who may have a wide range of answers to elicited model parameters and differing opinions on how best to represent the system. The importance of some of these assumptions and uncertainties can be addressed with sensitivity analyses but also requires a careful consideration of simulation results. The workshop identified gaps in research (e.g., the rate of patch cover increase and patch expansion) with great uncertainty among experts. There are general estimates of patch expansion (i.e., spread) over time based on aerial photographs (Olsson et al. ), but previous research has not considered patch cover increase. This STSM could be used to determine whether this gap in information is important in answering management questions (Cullinane Thomas et al. ). We also observed that managers were challenged in confronting the simplification and abstraction required in representing an ecological system as a model. Although the flexibility and complexity of the STSM framework allow for intricate modeling of ecological dynamics and management activities, modeling landscape conditions with multiple ecological and anthropogenic interactions requires simplification of the system.

Issues related to the integration of ST‐Sim with FARSITE included the age of the FARSITE program and its limitations (e.g., the FARSITE program cannot be invoked multiple times in parallel, drastically increasing computational time for Monte Carlo realizations). Integrating the external FARSITE model also added complexity to the STSM by increasing the number of vegetation‐based state classes from 14 to 23 with the addition of fuel load to the uninvaded classes (Appendix S2: Fig. S4) and the need to calculate additional fire‐related parameters. However, calling FARSITE, a fire behavior model, rather than using the simple fire simulation available in the ST‐Sim software, allows landscape characteristics to drive fire rather than trying to match historic fire characteristics such as fire size.

In recent years, it has become easier to integrate an STSM with external models as part of the SyncroSim framework. SyncroSim now includes a scripting interface through its rsyncrosim package for R and upcoming python package. Scripting interfaces increase the speed of developing and running models while minimizing the potential for error. The ability to link multiple models may promote the ability to address complex questions and customize models to the unique characteristics of a problem. Our work presented here is the first example of how an STSM and external model linked via SyncroSim can be used to quantify and compare landscape scale scenarios and successfully answered one of the types of questions that can be addressed with landscape simulation models by evaluating the scope of the buffelgrass invasion problem in SNP. Future work on this landscape simulation model could evaluate other questions this tool is suited for, including evaluating alternative management actions to identify which management strategy is most effective or conducting sensitivity analyses of parameters identified in this manuscript to prioritize research needs that affect management decisions (Cullinane Thomas et al. ). This framework (coupled STSM‐external models) could be used to address many complex management and ecological questions for landscapes allowing the STSM to include strengths from other models.

Acknowledgments

Funding to support this work was provided by the Natural Resources Preservation Project (NRPP) and the U.S. Geological Survey Invasive Species Program. We thank Kathryn Thomas and two anonymous reviewers for comments on earlier versions of this manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.