Introduction: Cooperation among managers of protected areas and federal multiple use lands with private inholdings to increase restoration success and economies of scale creates ecological and regulatory complexity best studied with state-and-transition simulation models (STSM).

Objectives: Project partners asked whether (1) agency budgets are sufficient to lower the dissimilarity between current vegetation and the reference condition, (2) prescribed burning improves bighorn sheep habitat without excess reduction of older subalpine forest classes, (3) vegetation treatments reduce uncharacteristic fire activity, and (4) climate scenarios affect vegetation treatments and bighorn sheep habitat?

Methods: Spatial STSMs were run for a 161,569-ha remote-sensed vegetation map. Partners selected two management scenarios crossed with three 50-year climate projections applied to 22 focal ecological systems with defined objectives, budgets, treatments cost and success rates, and metrics for vegetation and bighorn sheep.

Results: Drier climate suppressed fire but its greater precipitation variability increased fire and avalanche frequency. Treatments decreased area burned at lower elevations where non-native annual fuels facilitated fires. Ecological departure (vegetation dissimilarity) from reference conditions was unchanged when lower elevation seedings included introduced species, but ecological departure decreased with native species seeding. Prescribed fire and increased avalanche frequency increased bighorn sheep habitat suitability.

Implications for Practice: Treatment effectiveness increases by focusing treatment resources on <15 ecological systems per landscape; otherwise no single system is sufficiently funded. Prescribed fire at subalpine elevation can be used to increase young forage for and habitat suitability of bighorn sheep while temporally causing slight departure from reference conditions. Reduction of heavy woody fuels and non-native annual fuels at lower elevation followed by plant seedings reduced fire frequency at lower elevations. Traditional ecological departure might not be the best metric of success when introduced (i.e., non-native) species seedings replace vegetation classes not found during pre-European settlement as both equally contribute to ecological departure. Three systems might not require treatments because they increased ecological departure.

Study Area. The south Snake Range Area of Interest (AOI) of 161,569 ha includes lands primarily managed by the National Park Service (NPS) and BLM, with marginal USFS land to the north, while private lands include inholdings and larger tracts at the periphery of the AOI. The AOI is in east-central Nevada adjacent to the Utah border (Wheeler Peak at 38^o59’25.80”N, 114^o20’09.48”W). The vegetation of the south Snake Range is highly diverse as the elevation gradient is considerable (from 1,219-m to over 3,900-m). Vegetation combines systems typical of the central Great Basin that co-occur with the western edge of the North American monsoon.

Vegetation Mapping. One 161,569 ha image from WorldView 2/3 pan-sharpened 2 m resolution multispectral satellite imagery (Maxar Space Systems) was captured on June 22, 2022. About 33,514 ha covered the 2010 mapping of GBNP and a private parcel adjacent to GBNP, labelled the Keyhole Property. An additional 128,055 ha covered BLM and USFS-managed lands, and private lands in the AOI. Remote sensing was conducted with the software ERDAS Imagine® from Hexagon AB applied to WorldView imagery.

Two methods were applied to the analysis of the WorldView imagery: 1) change detection measuring spectral differences between the 2007 QuickBird (used to create the 2010 map; owned by Digital Globe, Inc. and retired 2015) and the 2022 WorldView imagery for GBNP and the Keyhole Property and 2) new remote sensing analysis applied to the rest of the AOI.

Change detection. Imagery interpreted in 2010 (QuickBird imagery from 2007) and in 2022 (WorldView) were compared using a variety of standard algorithms of change detection analysis including image subtraction, vegetation index comparison, principal components analysis, and manual evaluation (Lillesand & Kiefer 2000). Image-to-image analysis and manual evaluation resulted in a data layer indicating areas of spectral difference for further field investigation.

Field visits to all spectral differences were used to describe: 1) the potential agent of change (e.g., fire, restoration treatment, climatic/weather impact); and 2) current land cover classification, ecological system, and vegetation class. The land cover classification, ecological system, and vegetation class assignment was selected from the same discrete classification scheme employed in the 2010 mapping effort but used updated ecological system descriptions. In the field, seasonal and precipitation-caused spectral differences (not a transition change) needed to be distinguished from spectral differences caused by disturbance events and/or succession to other vegetation classes (true change). Areas with no spectral difference were assumed unchanged from 2010. We visited 445 field locations in GBNP and the KeyHole property for change detection in 2022.

Remote sensing for new mapping. Ground-verification field work was conducted from 8-24 July and 17-24 October 2022 following an unsupervised classification of imagery to identify unique spectral combinations of red, blue, green and near infrared and image texture to specific ecological systems and vegetation classes. On BLM, USFS, and allowed private lands, 9,927 rapid observation points (9,043 on BLM, 285 on USFS, and 599 on private lands) were collected to ensure that a representative portion of the project area was visited. At each rapid observation location, the ecological system, vegetation class, and at least two georeferenced photographs were taken. 

State-and-Transition Simulation Models (STSMs). STSMs are stochastic models of landscape change used to forecast outcomes of what-if scenarios. The landscape in an STSM consists of a discrete set of simulation cells classified into a discrete set of states. Those discrete sets of states were obtained from the remote-sensed ecological systems and vegetation classes. Simulation cells change over discrete timesteps according to a set of possible transitions caused by either natural (e.g., wildfire) or anthropogenic (e.g., plant seeding) processes operating within temporal or spatial constraints. All simulations were conducted using the ST-Sim package (Version 3.3.11) in SyncroSim (Version 2.4.18; https://docs.stsim.net/).

STSMs were spatially simulated for 50 years (2023 to 2072 with initial conditions in 2022). In ST-Sim each pixel was assigned an initial condition state, which is the combination of an ecological system and a vegetation class obtained from remote sensing that can either (a) age one timestep and stay in the same class, (b) age one timestep into an older class (i.e., succession), or (c) experience a probabilistic disturbance and transition to ≥1 other states, including the originating state. Additionally, land ownership, categorized as GBNP, BLM, USFS, or private can constrain management actions and budgets.

Fire has been the dominant ecological disturbance in the south Snake Range, although fires have not been common or large since federal agencies began mapping fires in the 1980s. The Phillips Ranch Fire (in 2000 at 1,025 ha), Black Fire (in 2013 at 1,912 ha), and Strawberry Fire (in 2016 at 1,885 ha) are the largest recent fires. These sizes reflected fire suppression management.

The other dominant process was herbivory expressed as livestock grazing on BLM, USFS, and private lands (primarily cattle and domestic sheep grazing mostly in the southern part of the AOI), wild horse grazing limited to the southeast corner of the AOI, and native herbivory by elk (Cervus canadensis), mule deer (Odocoileus hemionus), pronghorn (Antilocapra americana), and lagomorphs wherever relevant. Cattle, domestic sheep, wild horse, and native grazing were modeled separately.

Disturbances other than fire and herbivory included plant mortality from severe drought, short droughts, avalanches, native conifer invasion of shrublands, flooding, non-native annual grass invasion, exotic forb invasion, insect and disease outbreaks, wet years, very wet years, aspen (Populus tremuloides Michx.) clone loss, natural recovery of degraded vegetation, and succession.

Management Scenarios. A management scenario was a group of treatments and specific projected temperature and precipitation time series that defined a simulation theme. Scenario development was determined by the general guiding objectives of GBNP, BLM, and NDOW that set the tone for land management priorities informing 10-15 years of future treatment implementation. Two workshops were held for partners to establish guiding objectives and treatments for the management scenarios. Managed systems were also selected at the workshops. Guiding objectives were: (a) Map potential and current vegetation, and ecological condition as expressed by ecological departure from reference conditions (formerly Fire Regime Condition Class defined as the dissimilarity of distributions of vegetation class proportions per ecological system between current observations and reference conditions); (b) maintain overall condition of and restore degraded native upland and wetted ecological systems to reference conditions or desired future conditions given temperature and precipitation change; (c) maintain and enhance bighorn sheep habitat given temperature and precipitation change; (d) treat Wildland-Urban Interface (WUI) areas and reduce fuel loads to help protect human settlements and cultural resources in and around the project area from wildfire; (e) meet wilderness area objective of maintaining or enhancing wilderness characteristics using wildland fire for resource benefit and/or targeted prescribed fire; and (f) help NPS and BLM, and other partners meet objectives specified in management plans.

Six management scenarios were simulated for 50 years: two levels of management were factorially crossed with three forecasted 50-year temperature and precipitation scenarios. Management level one was a status-quo (control) scenario (named “custodial”) that included the maintenance of current fire-suppression activities and, where applicable, livestock grazing practices, but did not use any mechanical or prescribed burning actions. Management level two was “preferred” management that represented an ambitious level of vegetation treatment implementation. While $0 per year was assigned as the expenditure for the custodial scenario, the preferred scenario was limited to a maximum average annual expenditure of about $200,000 on GBNP lands and $400,000-$600,000 on BLM lands from 2025 to 2039, with low level maintenance management thereafter. No management was simulated in private lands, so no funding was assigned to private lands, although fire and invasive species moved in and out of private lands to adjacent lands, or USFS-managed lands. Each action was assigned a cost per area and other implementation attributes embedded in the simulation library included success rate, vegetation class outcomes for success and failures, equipment use constraints based on slope, and mechanical treatment exclusion areas such as wilderness areas.

Three climate scenarios were statistically forecasted from time series of minimum and maximum temperatures and precipitation from 1945 to 2022 data obtained from the Parameter-elevation Relationships on Independent Slopes Model (PRISM; continuation of historic temperature and precipitation) and Localized Climate Analogs (LOCA) downloads of ACCESS1 (hotter and much drier than current conditions) and HadGEM2-ES (hotter and periodically wetter than current conditions). LOCAs were downloaded for location 38°56'46.2"N 114°15'26.8"W. Methods of statistical forecasting of replicated data and translation to Standard Precipitation-Evapotranspiration Index (SPEI; SPI initially published by Hayes et al. 1999) time series used in simulations. The six simulated scenarios were: custodial and PRISM, preferred and PRISM, custodial and ACCESS1, preferred and ACCESS1, custodial and HadGEM2-ES, preferred and HadGEM2-ES. The same actions were implemented at the same planned rates in each scenario; however, planned rates were not equal to realized rates as stochasticity caused differences in implementation among scenarios and replicates within each scenario.

Replication and Temperature and Precipitation Variability. The purpose of simulating future temperature and precipitation time series was to introduce temporal variability in dominant ecological processes. Each management scenario (custodial and preferred) used the same replicated time series to introduce temporal variability into the expression of ecological disturbances. Downloaded monthly precipitation and monthly minimum and maximum temperatures for PRISM, ACCESS1 and HadGEM2-ES were used with a stochastic weather generator (SWG) to statistically replicate precipitation, minimum temperature, and maximum temperature time series 10 times over 50 years. Variability directly affected processes in ST-Sim through temperature and precipitation, or indirectly mediated processes through the SPEI. Methodology of projected temperature and precipitation variability effects on ecological processes per replicate of SPEI was complex.

Metrics of Success. We limited results to three landscape level metrics: ecological departure (ED) from reference conditions, map of area burned frequency, and expert-driven bighorn sheep habitat suitability.

Ecological Departure. ED is the primary metric of the US LANDFIRE project to measure dissimilarity from ecological system reference conditions and prioritize regional fuels management funding (termed Vegetation Departure Condition by LANDFIRE, https://www.landfire.gov/). ED is the dissimilarity between the observed (O) distribution of vegetation class proportions from remote sensing, and the expected (E) distribution of vegetation class percentages in the reference condition per ecological system (i.e., generally assumed to be the pre-European settlement condition):

^^            ED = 1 - sum(min{O_i,E_i})                                                                                                      Eq. 1

where i = 1,…,N classes per ecological system. Lower values indicate that vegetation was closer to the reference condition and may not require intensive or any management interventions. ED is a powerful metric because the exact O_i’s causing departure can be identified by partial sums of Eq. 1.

Expert-driven Bighorn Sheep habitat Suitability. The habitat suitability index for the AOI’s Rocky Mountain bighorn sheep population was designed in R code based on two partner workshops supplemented with the literature and follow-up expert consultation. We used expert workshops because there was not yet a statistical habitat suitability model for the Snake Range population. Nine covariates were considered to build resource selection functions (RSF) partitioned between growing season and winter habitats that both shared escape terrain and distance to roads as RSFs (all equations and details in Supplement S6). A RSF is the mathematical or statistical equation relating an independent covariate to a partial measure of habitat suitability response.

Growing season suitability was defined by (1) distance to roads, (2) distance to escape habitat, and (2) growing season habitat attributes:

1.       Distance from the nearest classified road pixel (hiking trails were not included) was measured using the distance function in the ‘terra’ package. Distances greater than 400 m were given a value of 1.0 (i.e., no effect on bighorn sheep), but the RSF was increasingly smaller closer to roads.

2.       Escape terrain is defined by steep terrain above 30° (60% slope) calculated using the terrain function with an 8-pixel neighbor window in the ‘terra’ package. Pixels above 30° were retained and distance to those pixels to a location was estimated using the distance function in ‘terra’. Distances within 250 m were considered optimal.

3.       Growing season habitat was defined by forage, topographic roughness, and elevation. (a) Growing season forage was based on ecological systems and vegetation classes. Ecological systems used by bighorn sheep were identified by partners during workshops. Vegetation classes were assigned values of 1.0 (highest forage value for bighorn sheep), 0.75 (some forage value), or 0.0 (no forage value). Generally, early seral classes had the highest forage value. (b) Topographic Roughness Index (TRI) was calculated using the “Terrain” function in the ‘terra’ package. Higher TRI values were preferred by bighorn sheep. (c) Elevation was used to constrain growing season habitat, where elevations above 2,438 m were most preferred.

The non-adjusted growing season habitat (Habitat_s) was the average of forage, topographic roughness, and elevation RSFs. The adjusted and final growing season habitat suitability index was found by averaging Habitat_S with the distance to escape terrain RSF and distance from road RSF.  

Unadjusted winter habitat was estimated using four parameters: 1) aspect, 2) elevation, 3) winter forage, and 4) average snow depth:

1.       Aspect was used to identify areas where winter forage are more readily available and thermal conditions are preferred. The highest values for the aspect RSF were southerly facing estimated with the terrain function in the ‘terra’ package. The RSF was calculated as a deviation from a pure south aspect.

2.       The winter elevation RSF was characterized by the toe slopes of the project area with areas above ~1,830 m being preferred.

3.       The winter forage RSF was created in a similar process as the growing season forage RSF. Forage values assigned to each vegetation type were 1.0 (for highest forage value), 0.75 (for moderate forage value), and 0.0 (for no forage value).

4.       Partners identified that snow deeper than 0.9144 m was an important barrier for sheep to access forage. Average winter snow depth data between 2004 to 2022 were obtained from the National Operational Hydrologic Remote Sensing Center and interpolated to the project area.

As with the growing season habitat, winter habitat was calculated by averaging the four winter RSFs and then the final winter habitat suitability index was found by averaging the winter habitat value, distance to road RSF, and distance to escape terrain RSF.

Unique among metrics of success, the bighorn sheep habitat suitability R code was coupled with ST-Sim through the ‘rsyncrosim’ package and guided vegetation management actions. Each year, ST-Sim outputted vegetation rasters to the R script, which estimated suitability as a value from 0 to 1 per pixel that was returned to ST-Sim. Intermediate suitability values per pixel increased the probability of use of some restoration actions creating young forage that could cause habitat uplift. Already highly suitable pixels required no uplift, whereas poor suitability pixels would remain so for physical reasons, such as low elevation, unrelated to forage availability. All other treatments that did not contribute to bighorn sheep management were unconstrained by suitability and the treatment implementation algorithms of ST-Sim operated without adjustment for bighorn sheep suitability.