Study area and ecological dynamics

The RGdNNM covers 1 255 km² (310 000 acres) in northern New Mexico and is managed by the BLM. In RGdNNM, elevation ranges from 1 830 m to 3 100 m. Precipitation ranges from 200 mm to 900 mm annually (averaged over 1981–2010) with higher precipitation trending south and on the volcanic cones (PRISM Climate Group 2004). Soil temperatures span frigid to mesic from north to south (Bauer 2011). The RGdNNM is situated within an intermountain valley. Prominent landscape features include rolling plains, volcanic cones, remnant basalt flows, playas, and a rift gorge. RGdNNM is geomorphically and hydrologically characterized by the predominance of closed basins, containing many small playas, with only two outlets into the adjacent Rio Grande (Johnson and Bauer 2012). The soils are young Aridisols formed from alluvial fill and striated with layers of fractured basalt at varying depths from Tertiary and Quaternary volcanic cones and vents (Johnson and Bauer 2012). Three distinct climatic-geologic units (major land resource areas [MLRAs]) converge within RGdNNM (Fig. 2). The convergence of MLRAs indicates high spatial variability in precipitation patterns and temperature regimes across RGdNNM (Salley 2016).

Figure 1.

Integration of an inductive workflow into ecological site description development. Our proposed multivariate workflow expands on the ecological site concept development stage of the ecological site description development process (after Moseley et al. 2010). The ecological site concepts developed with the multivariate workflow should be tested with additional data collection and stakeholder workshops.

A desert shrub-grassland mosaic dominates the plains with an understory of blue grama grass (Bouteloua gracilis Willd. ex Kunth), western wheatgrass (Pascopyrum smithii Rydb.), and other perennial grasses and forbs. Resprouting subshrubs broom snakeweed (Gutierrezia sarothrae Pursh) and Greene's rabbitbrush (Chrysothamnus greenei A. Gray) are common in the southern and northern monument, respectively. Wyoming big sagebrush is the dominant shrub in the southern RGdNNM. Efforts to quantify relationships among environmental controls on Wyoming big sagebrush plant communities have been met with varying success, which is attributed to the broad ecological amplitudes and ecotypic variation of Artemisia-steppe species (Davies et al. 2007). Wyoming big sagebrush tolerates the driest soils and lowest elevations of the big sagebrush subspecies and is also considered to be the least resilient (i.e., less able to regain structure and function post disturbance) and resistant (i.e., less able to retain structure and function despite disturbance) of the big sagebrush subspecies (Chambers et al. 2014; Miller et al. 2014). Winterfat (Krascheninnikovia lanata Pursh) is the dominant woody plant in the northern RGdNNM. Winterfat grows in dry soils, typically with accumulations of calcium carbonate or other salts (Woodmansee and Potter 1971). Twoneedle piñon (Pinus edulis Engelm.), oneseed juniper (Juniperus monosperma Engelm.), and Rocky Mountain juniper (Juniperus scopulorum Sarg.) trees occupy patches on the plains and are common on the volcanic cones. Pinyon-juniper (PJ) communities can comprise woodlands, typically on shallow, coarse soils with a sparse herbaceous understory; savannas, with a low to moderate density of trees and continuous herbaceous understory, typically on moderately deep to deep soils; and wooded shrublands, where tree density is variable and the shrub stratum is well developed, occurring on a variety of substrates (Romme et al. 2009). Mixed conifer forest is present at higher elevations on the volcanic cones and is composed of white fir (Abies concolor Gord. & Glend.); Rocky Mountain juniper, ponderosa pine (Pinus ponderosa Lawson & C. Lawson); and Douglas fir (Pseudotsuga menziesii Mirb.). In MLRA 48A, which occupies a relatively small area of RGdNNM, black sagebrush (Artemisia nova A. Nelson) is the dominant woody plant species, with an understory of perennial grasses, broom snakeweed, and Greene's rabbitbrush.

As in other arid and semiarid systems, spatial variability in plant communities is determined by patterns of soil water availability, as influenced by landscape position, soil texture, and subsurface soil properties (Costantini et al. 2016). In shrub-grass mosaic landscapes like that of the RGdNNM, loamy soil textures and the presence of a clay-rich argillic horizon near the soil surface are often associated with increased herbaceous cover due to increased availability of shallow soil water that favors shallow-rooted grasses and forbs (Sala et al. 1997; Davies et al. 2007; Pennington et al. 2017). Deep-rooted woody shrubs and trees are able to access moisture deeper in the soil profile, as well as moisture that gathers in fractured bedrock in shallow soils (Duniway et al. 2010) Bulgamaa et al., 2020. Numerous, and often interacting, drivers including climatic trends, drought episodes, altered fire regime, and livestock grazing can cause persistent increases in woody species in shrub-grass mosaics, as the deep roots of woody species are able to access soil moisture reserves at a depth that perennial grasses cannot (Miller 2005; Archer et al. 2017; Bestelmeyer et al. 2018). A review of field notes from late 1800s' cadastral surveys indicated that the plant species and groups dominant in the landscape today (i.e., big sagebrush, bunchgrasses, and PJ) were also dominant at the time of survey (Bureau of Land Management, accessed January 2019). Cadastral survey field notes and land managers also indicate that livestock grazing has been ongoing throughout RGdNNM since the late 1800s, though long-term data regarding the intensity and spatial distribution of grazing are unavailable (Bureau of Land Management, accessed January 2019).

Figure 2.

Monitoring plots (n = 213) within the Rio Grande del Norte National Monument (RGdNNM). Major land resource area boundaries are shown in red; the monument boundary is in green. Average annual precipitation is indicated by colored shading.

Records are available for chemical, mechanical, and fire vegetation treatments in RGdNNM between 2002 and 2014 that were intended to reduce the cover and density of Wyoming big sagebrush or PJ trees and to increase herbaceous cover (Bureau of Land Management 2004; Bureau of Land Management 2012). A review of RGdNNM vegetation treatment responses showed trends consistent with responses in Wyoming big sagebrush communities across the western United States (Traynor et al. 2020). Mechanical shrub removal treatments in Wyoming big sagebrush ecosystems can increase annual plant cover and produce communities dominated by a single plant functional group (Prevéy et al. 2010; Ripplinger et al. 2015; Davies et al. 2016). This is a particular concern when the herbaceous, perennial understory of sagebrush steppe has already been depleted through overgrazing or another disturbance and along warmer/drier environmental gradients where community resilience and resistance to invasion are lowest (Prevéy et al. 2010; Davies et al. 2012; Chambers et al. 2014). Chemical treatments tend to increase perennial grass cover and total foliar cover when environmental conditions are favorable for grass growth in Wyoming big sagebrush ecosystems, although these treatments often do not permanently reduce sagebrush cover (Sneva 1972; McDaniel et al. 2005). In RGdNNM, treatments to remove PJ trees were more difficult to interpret (Traynor et al. 2020). Invasive plant cover is thought to increase in piñon-juniper removal treatments, particularly when native herbaceous cover is below 20% before treatment (Chambers et al. 2014).

Monitoring data

Data were acquired from 213 monitoring plots sampled within RGdNNM between May and October from 2011 through 2018. Sampling was conducted by ecological monitoring programs implemented by the NRCS (Landscape Monitoring Framework [LMF]) and the BLM (Assessment, Inventory, and Monitoring [AIM]). We refer to these plots collectively as “core methods plots.” More than 27 000 core methods plots have been described nationally across the AIM and LMF programs, with a subset revisited on a 5-year rotation. At core methods plots, data were collected following standardized methods of Herrick et al. (2018) that are “core” to the monitoring programs: line-point intercept, vegetation height, canopy gap intercept, soil stability, species inventory, and plot characterization. Plot locations were generated from spatially weighted stratified random sample designs (Herrick et al. 2010; Toevs et al. 2011). Line-point intercept, vegetation height, canopy gap intercept, and soil stability were collected along three transects 25–50 m in length. Species inventory was collected across entire plots. A soil pit was excavated at each plot to a depth of 70 cm or until a restrictive layer was reached, and horizons were numbered but taxonomic horizon designations were not assigned. Horizon depth, soil texture as determined by hand texturing in the field, soil structure, soil color, percent clay, rock fragment content, and effervescence class in response to 1 M HCl were recorded for each soil horizon. Photo points were taken for each transect and for the soil pit. Available water holding capacity was calculated with monitoring plot soil pit data (J. Williamson, personal communication; Saxton et al. 1986).

Supporting data

Additional soil data were acquired from SoilGrids (Hengl et al. 2017). SoilGrids predicts physical and chemical soil properties at multiple depths and was generated with machine learning algorithms and environmental covariates derived from digital elevation models, remote sensing data, and land cover classifications. Soil property predictions are presented at a 250-m resolution. Mean soil chemical and physical properties (soil pH, soil bulk density, percent soil clay, percent soil sand, soil organic carbon) were extracted for plot locations at depths of 0 cm, 5 cm, 15 cm, 30 cm, and 60 cm to supplement the core methods soil pit data. A geographic information system (ArcMap, version 10.2.2) was used to calculate elevation and landform variables from the US Geological Survey National Elevation Dataset 10-m digital elevation model (NRCS Geospatial Data Gateway, accessed 2019) and average annual precipitation from the PRISM 30-yr-average dataset (PRISM Climate Group 2004) at the plot level using ArcMap.

To validate the ecological site and vegetation state concepts derived from the monitoring plot data, a low-intensity traverse was completed in the study area from May to August of 2019. The low-intensity traverse is used in the initial phase of ecological site development to observe broad environmental gradients and general vegetation-soil-landform relationships (Caudle 2013). Rapid characterizations of plant communities, slope, landform, slope shape, and soils are recorded. Observation points are often selected opportunistically while traversing major roadways through the study area. During the low-intensity traverse, qualitative observations of landforms, MLRA boundaries, soils, and vegetation communities were recorded. These observations, in addition to literature review as cited, were used to describe the general ecological dynamics of RGdNNM (section 2.1). A complete list of variables included in the analysis and their sources are provided in Table 1.

Polygons delineating vegetation treatments implemented by the BLM between 2002 and 2014 were used to assign core methods plots to the following categories: no treatment (n =135), fire treatment (n = 23), chemical treatment (n = 8), or mechanical treatment (n = 47). Fire treatments targeted piñon-juniper trees through thin/burn and broadcast burn approaches. Chemical treatments targeted big sagebrush with Tebuthiuron through aerial application. Mechanical treatments targeted big sagebrush with disking, brushhogging, shaving, or drill seeding treatments, with some areas treated with multiple mechanical treatment types. The time in years between the most recent treatment and data collection was calculated for each plot in a vegetation treatment. Observations made during low-intensity traverse sampling indicated additional, suspected treatment areas that were not described by the treatment history polygons.

Analysis

Statistical analyses were run in R 3.5.0 (R Core Team 2018) and JMP (JMP, version 13 1989–2021). Analytical steps were ordered to first identify vegetation communities based on functional group abundance and plot structural attributes and then describe important abiotic conditions distinguishing vegetation communities. First, plant species were classified into functional groups within genera, based on longevity, structure, and photosynthetic pathway when applicable (e.g., C3 vs. C4 grasses; Gondard et al., 2003). Fuzzy cluster analysis was run on functional group abundance, percent bare soil, and proportion of the plot with large canopy gaps (i.e., canopy interspaces > 100 cm that have been associated with perennial species loss and increased erosion risk; Derner and Whitman 2009; Webb et al. 2014). Percent bare soil and proportion of the plot with large canopy gaps were included in the vegetation community analysis due to their importance as indicators of attributes of ecosystem function in rangelands (Pyke et al. 2002; Pellant et al. 2020). Fuzzy clusters were visualized in a principal coordinate analysis (PCoA) ordination. The monitoring plots were subset so that clusters included only those plots with a membership value of > 0.6, to create clusters with strong within-cluster homogeneity. Modal concepts of vegetation communities were defined by summarizing the minima, maxima, and mean for functional group abundances and structural attributes within clusters. Modal concepts were validated with plot photos. Plots with membership values lower than 0.6 to any cluster were assigned manually to existing modal concepts or to new modal concepts when they represented a novel vegetation community.

BIOENV analysis was run on the functional group abundance dissimilarity matrix and the set of abiotic variables to identify the abiotic variables with the highest correlation to the functional group abundance data (Clark and Ainsworth 1993; Spiegal et al. 2014). A classification tree was generated with the JMP Partition platform (JMP, version 13) using vegetation modal concept class as the response variable and the abiotic variables identified through BIOENV analysis (see Table 1) as predictor variables to generate ecological site concepts for sets of vegetation communities. A tree was first generated using K-fold cross-validation and then interactively pruned. The variables and breakpoints selected by the final classification tree were recorded and used to subset the monitoring plots into ecological site concept groups. The ecological site concepts identified through the classification tree were stratified by MLRA. Kruskal-Wallis (95% CI) and post-hoc Wilcoxon tests were used to evaluate significant differences in structural attributes of vegetation communities and abiotic variables among ecological site concepts. Within treatments, functional groups and structural indicators were summarized across treatment type and by ecological site concept. A Bray-Curtis dissimilarity matrix, which calculates shared abundance divided by total abundance, was used for all functional group abundance and structural attribute analyses and was selected for its compatibility with community ecology analyses (McCune et al. 2002). AIM and LMF monitoring data were accessed from the national TerrADat database using the R package “terradactyl” (McCord and Stauffer 2020).  Appendix A (mmc1.docx) provides additional details regarding methods, including setup and implementation of BIOENV and the classification tree.

Table 1

Variables used in analysis, with descriptions and sources. All soil and site variable were included in the BIOENV analysis. Variables that were selected for use in the final workflow (to define vegetation communities or in selection by classification tree analysis) are indicated by an “X.” Soil and site variables not selected by BIOENV analysis were summarized for each ecological site concept.

Table 2

Vegetation community clusters identified through fuzzy cluster analysis (clusters 1–16) and qualitative assessment of “fuzzy” plots (clusters 17–20). Vegetation communities are named for dominant functional groups and structural attributes. Column two lists functional groups that occurred on 100% of plots and are listed in descending order by mean percent foliar cover.

Vegetation communities

The analysis of functional group abundance and plot structural attributes resulted in a total of 20 vegetation communities (Table 2). Fuzzy cluster analysis produced 16 clusters, and 139 plots had membership values > 0.6 to the clusters (Fig. 3). When “fuzzy” plots (membership values < 0.6; n = 74) were reviewed, 56 were assigned to the existing clusters (clusters 1–16) and 14 plots were assigned to one of four new vegetation communities (clusters 17–20) that were underrepresented by the plot data. Four fuzzy plots were removed from subsequent analysis because they spanned multiple vegetation communities, as determined through plot photo validation. Of the final vegetation communities, three were dominated by Wyoming big sagebrush, two were dominated by winterfat, two were dominated by trees, four were dominated by western wheatgrass, two were dominated by black sagebrush, five were dominated by perennial bunchgrasses, and two were dominated by equal proportions of mixed shrubs, grasses, and subshrubs (see Table 2). When the vegetation communities were grouped by their dominant plant functional groups (e.g., mean structural attributes tested across big sagebrush/C4 bunchgrass, big sagebrush/piñon-juniper, and big sagebrush/bare soil), mean plot structural attributes varied significantly for all functional group assemblages except the forested communities, indicating structural vegetation differences even when functional group compositions are similar (Table 3).

Ecological site concepts

BIOENV analysis identified elevation, average annual precipitation, soil depth, predicted soil pH at a depth of 30 cm, soil depth to the maximum effervescence response, and maximum effervescence response class as the abiotic variables with the strongest correlation (0.51) to the functional group abundance dissimilarity matrix. These variables were used as predictors in a classification tree analysis with vegetation modal concept as the response variable. K-fold cross-validation (K = 10) produced a tree with 16 splits and an R² of 0.512. The tree was pruned interactively to 10 splits and 11 ecological site concepts (R² = 0.429; Fig. 4). Pruning was based maintaining plant community assemblages (terminal nodes) that matched the vegetation community groups as described in section 2.1, defined by low-intensity traverse observations and literature review, and consolidating splits that produced similar assemblages.  Appendices A (mmc1.docx) and  B (mmc2.docx) describe this approach in further detail. Subsequently, two concepts (terminal nodes C1, C8; see Fig. 4) that spanned across MLRAs 36, 51, and 48A were stratified by MLRA boundaries. Stratification by MLRA was important to ecological site concept development, as ecological sites are constrained by a single MLRA. The Mesic Loamy, Mesic Alkaline, Frigid Loamy, and Frigid Alkaline plains were then consolidated where they had been split by the classification tree. An ecological site concept was added to represent the Rock Outcrop site identified through the qualitative review of fuzzy plots, which showed a unique plant community and stretches of exposed bedrock or boulders. Each ecological site concept was predictive of a suite of vegetation communities, with the exception of Rock Outcrop, which was represented by a single vegetation community. A qualitative spatial analysis of the ecological site concepts revealed that the concepts were constrained by discrete landforms, which were used to name the concepts: one for the gorge rim, four for the intermountain plains, one for rock outcrops within the intermountain plains, one for fan remnants, one for alluvial fans, two for the volcanic cone footslopes, and one for the volcanic cone slopes (Table 4).

Figure 3.

PCoA ordination and fuzzy cluster analysis of core methods plots (n = 213) by abundance of functional groups. The plots are represented by colored wedges. The color of the wedge indicates fuzzy cluster (1–16) and the size of the wedge indicates membership value (0–1). Plots represented by multiple wedges have membership to multiple clusters, with membership values proportional to the size of the wedges. The ordination plot axes (Dim. 1 and Dim. 2) are labeled with their strongest predictors, based on permutations (see  Appendix A (mmc1.docx) for methods and  Appendix B (mmc2.docx) for results). Descriptions of clusters are given in Table 2.

Table 3

Means and standard deviations of structural indicators (percent bare soil and proportion of canopy gaps > 100 cm) for vegetation communities. Letters in superscript denote vegetation communities dominated by the same functional group with structural indicators that differ significantly as established by Kruskal-Wallis and post-hoc Wilcoxon testing. Where no lettering is applied, differences are not significant.

Though MLRA was not included as a predictive variable in the final classification tree analysis, it was used as a predictive variable by subsetting the ecological site concepts when they spanned across MLRAs. When the ecological site concepts were divided by MLRAs, the resulting suites of vegetation communities that appeared similar at the broad functional group level were differentiated by species. For example, the Alkaline Plains ecological site concepts across the entire study area were represented by a shrub-grass matrix. In MLRA 36 (Mesic Alkaline Plains) this is a Wyoming big sagebrush-grass matrix; in MLRA 51 (Frigid Alkaline Plains) this is a winterfat-grass matrix. Black sagebrush is the dominant shrub species in MLRA 48A. Spatial patterns of dominant plant functional groups suggest that the boundary between MLRAs 36 and 51 could be shifted northward, and the boundary of MLRA 48 could be shifted eastward to better represent ecologically meaningful boundaries. All of the concepts representing the volcanic cone landforms occurred in MLRA 36 with the exception of one Deep Footslopes plot in MLRA 51; this is likely an effect of minimal sampling of the volcanic cone landforms in MLRAs 51 and 48 by the monitoring data.

Figure 4.

Results from classification tree analysis with vegetation modal concepts as the response variable and the abiotic variables identified through BIOENV analysis as predictor variables (predicted soil pH at a depth of 30 cm, soil pit depth (cm), elevation (m), maximum effervescence response class in soil pit, and depth to maximum effervescence response class in soil pit (cm). The numbers (e.g., C7) representing the terminal nodes correspond to the vegetation community concept with the highest proportion of plots represented. The terminal nodes are annotated with the ecological site concept that they contribute to.

Table 4

Mean values for elevation, average annual precipitation, slope, soil pit depth, soil depth to maximum effervescence response, predicted soil pH, and maximum effervescence class across ecological site concepts. Numbers in superscript denote ecological sites that differ significantly.

Vegetation community–ecological site concept relationships

Suites of vegetation communities were summarized across the ecological site concepts where they occurred (Table 5). Wyoming big sagebrush–dominated communities were present in 5 of the 10 ecological site concepts, across a variety of landforms: Gorge Rim; Mesic Alkaline Plains; Fan Remnant; Alluvial Fan; and Deep Footslopes. No foliar cover of Wyoming big sagebrush was recorded in the cold, dry basin that occupies most of the northern RGdNNM in MLRA 51. Winterfat-dominated communities occurred on the Frigid Alkaline Plains site concept, which was defined by limy, calcic soils, either where calcium carbonate was recorded within 25 cm of the soil surface or predicted soil pH at a depth of 30 cm was > 7.9. Black sagebrush–dominated communities occurred on the Fan Remnant ecological site concept where temperatures were relatively cool and average precipitation was relatively high. Forest vegetation communities occurred on alluvial fan, volcanic cone footslopes, and volcanic cone slopes. Forest communities co-occurred with perennial bunchgrasses where soils were shallow (Shallow Footslopes; mean soil depth 26 cm) and with predominately cool-season rhizomatous grasses where soils were deep (Deep Footslopes; mean soil depth 69 cm). There were no shrub-dominated communities supported on the Loamy Plains ecological site concepts, where maximum calcium carbonate accumulation occurred at a soil depth > 25 cm. The plots classified to the Rock outcrop ecological site were dominated by warm and cool season bunchgrasses, subshrubs, and annual forbs.

Vegetation treatments

The vegetation communities associated with known vegetation treatments were summarized across the study area and within ecological site concepts (Fig. 5; Fig. B.1). Across all ecological site concepts, plots in chemical shrub removal treatments (n = 8) had high perennial grass cover (68% ± 19%), low percent bare ground (15% ± 9%), and a low proportion of large (> 100 cm) canopy gaps (9% ± 6%). They were classified into three vegetation communities (see Fig. 5). Plots within mechanical shrub-removal treatments (n = 47) were classified into nine vegetation communities, five of which were characterized by either presence of invasive species, high percent bare ground, or high proportion of large canopy gaps. Within all mechanical shrub removal treatments, mean perennial grass cover was 34% (± 16%), mean percent bare ground was 30% (± 12%), and mean proportion of large canopy gaps was 30% (± 19%). Cover of broad functional groups and structural attributes are summarized in  Appendix B Table B.1) (mmc2.docx)(. Plots in fire treatments that occurred on volcanic cone landforms (n = 17) were classified into tree-dominated vegetation communities with the exception of plots monitored less than 3 yr post treatment. In these cases, the Post-Fire/Shrub/Subshrub/Grass vegetation community or the Perennial Bunchgrass/Subshrub community was recorded, neither of which had a notable tree component. Plots in fire treatments occurring on the Fan Remnant site concept (n = 6) were classified into communities dominated by perennial grasses, shrubs, and subshrubs. Within all fire treatments, mean perennial grass cover was 23% (± 20%), mean percent bare ground was 18% (± 8%), and mean proportion of large canopy gaps was 29% (± 16%) (Table 6).  Appendix B (mmc2.docx) provides additional details regarding results.

Differentiation of vegetation community assemblages

Our workflow produced ecological site concepts and associated assemblages of potential vegetation states that differ from one another in topography, geomorphology, climate, and soil properties. Plant community composition is often used to define states in STM applications but may not capture important functional attributes that differentiate alternative states. Biotic thresholds are often crossed before abiotic thresholds, and states can comprise multiple community phases that may have substantial variation in plant composition but not in ecological process function (Stringham et al. 2003; Briske et al. 2005). Including measurements of percent bare ground and proportion of large intercanopy gaps provided additional information about the influence of plants on surface hydrology, soils, and nutrient cycling that plant composition alone may not represent (Bestelmeyer et al. 2009). As ecosystem processes such as nutrient cycling, energy flow, and recovery mechanisms are difficult to measure directly, indicators representing attributes of ecosystem function, such as percent bare ground and intercanopy gap sizes, are used in rangeland ecosystems to interpret soil/site stability, hydrologic function, and biotic integrity (Pyke et al. 2002; Pellant et al. 2020). Statistically significant differences in percent bare ground and proportion of large intercanopy gaps among vegetation communities dominated by the same plant functional group (e.g., Wyoming big sagebrush/C4 grass and Wyoming big sagebrush/bare soil) suggest that these communities may represent true alternative states with functionally significant differences. The suites of plant communities within ecological site concepts identified through our workflow can be validated by their similarity to previously described plant communities and alternative states. Two different PJ community types described by Romme et al. (2009) for the western United States, including the intermountain west, southwest, and southern Rocky Mountain regions, are supported on lowland and upland ecological site concepts. The Gorge Rim, Mesic Alkaline Plains, and Alluvial Fan support plant communities similar to the wooded PJ shrubland (Romme et al. 2009), wherein the community contains a well-developed shrub (Wyoming big sagebrush) stratum with a tree component that increases and decreases with climatic conditions and disturbances. Also consistent with Romme et al. (2009), the Shallow Footslopes, Deep Footslopes, and Volcanic Cones ecological site concepts support persistent PJ woodlands with sparse understories on upland sites that inherently favor tree growth. The assemblages of vegetation communities assigned to the Gorge Rim and Mesic Alkaline Plains closely mirror alternative states posited in existing STMs for Wyoming Big Sagebrush (8–12" PPT) and Big Sagebrush (12–14" PPT) by Chambers et al. (2014) and developed with data from Utah, Nevada, California, and Oregon in sites with Mesic/Aridic to Xeric soil temperature/moisture regimes, though not all of the states in those existing STMs were represented by our data.

Typically, inductive approaches to ecological site concept development have first identified ecological site units based on abiotic characteristics and potential vegetation states have been identified second, within hypothesized ecological site units (e.g., Spiegal et al., 2014; Ratcliff et al. 2018; Svejcar et al. 2018). While this approach was tested in the development of our workflow, we found that identifying vegetation communities first and subsequently testing their distribution across abiotic units (ecological sites) produced more interpretable results. This could be due to the strong effect of land use and management history on plant communities and lack of detailed spatial and temporal data documenting management history before 2002; when classification tree analysis on plant communities included treatment type as a predictor variable, it was selected as the first break (most predictive variable), dividing monitoring plots in known mechanical treatments from plots in fire, chemical, or untreated areas and overruling all other climatic, physiographic, and soil variables. The importance of incorporating land use legacies into STMs as drivers of state change, combined with lack of management history data as in RGdNNM, reinforces the necessity of pairing inductive and deductive methodologies in ecological site development (Knapp et al. 2011; Kachergis et al. 2013).

Table 5

Ecological site concepts that are distinguished by landforms, precipitation regimes, soil properties, major land resource areas (MLRAs), and vegetation community dynamics.

Figure 5.

Vegetation communities visualized by treatment type: chemical (n = 8), fire (n = 23), or mechanical (n = 47).

Table 6

Means and standard deviations of broad functional group foliar cover (%) and cover of structural indicators (%), reported by treatment type.

Vegetation treatment evaluation

Vegetation treatments in RGdNNM have been implemented to mitigate woody plant encroachment since at least the 1950s (Pilliod and Welty 2019; Traynor et al. 2020). Our analysis suggests that woody plant encroachment may be occurring by Wyoming big sagebrush into historically grass-dominated sites and by PJ trees into historically shrub-dominated sites. The Mesic Loamy Plains and Frigid Loamy Plains ecological sites, which are characterized by relatively greater depth to soil carbonates in the soil profile (> 26 cm) and relatively low slope (< 10%), support perennial grass-dominated communities, one of which has a shrub component (see Table 5). Due to the predominance of closed basins and playas in RGdNNM, areas of low slope likely have an accumulation of fine soil particles and increased moisture from run-in that favor perennial grass growth (Duniway et al. 2010; Johnson and Bauer 2012; Williams et al. 2016). The occurrence of the western wheat/big sagebrush community on these ecological sites may be an example of infilling of sagebrush into a perennial grass-dominated reference state due to drought, overgrazing, and feedbacks (e.g., reduction in herbaceous growth and decreased infiltration) perpetuated by altered resource availability and vegetation patterns (Miller 2005; Archer et al. 2017; Bestelmeyer et al. 2018). The Mesic Alkaline Plains, Gorge Rim, and Fan Remnant ecological sites are characterized by calcic soils and balanced grass/shrub matrices. Each of these sites also supports the big sagebrush/tree vegetation community, which may be an example of piñon-juniper encroachment in sites where the reference state is a sagebrush-perennial grass steppe (e.g., Miller 2005; Chambers et al. 2014). There are multiple drivers of PJ encroachment, including grazing, climatic variability, fire exclusion, and rising atmospheric CO₂ concentrations (Romme et al. 2009; Miller et al. 2019; Reinhardt et al. 2020). Our hypothesized shrub-steppe reference states reinforce observations of favorable shrub conditions occurring with shallow (< 26 cm) or increased soil carbonates (Bestelmeyer et al. 2009; Svejcar et al. 2018).

Our classification of plots to vegetation communities within treatment types is consistent with previous studies of the treatments in RGdNNM using monitoring data (Traynor et al. 2020), as well as with treatment results reported for Wyoming big sagebrush communities across the western United States (Sneva 1972; McDaniel et al. 2005; Davies et al. 2016). Chemical and mechanical shrub-removal treatments have no natural analogue and thus produce longer-term (decades- or centuries-long) legacy effects on plant community structure and function when compared with natural disturbance regimes (Ripplinger et al. 2015). Concurrent with studies of shrub removal treatments in sagebrush-steppe ecosystems, particularly Wyoming big sagebrush-steppe, we found mechanical treatments in RGdNNM were associated with greater foliar cover of invasive plant species when compared with other treatment types (Prevey et al. 2010; Ripplinger et al. 2015; Davies et al. 2016). Of the two vegetation communities that were characterized by substantial presence of invasive species, Invaded Grassland (n = 8) and Invaded Bunchgrass (n = 15), 75% and 73% of the plots, respectively, were in recorded mechanical treatments. Invasibility of sagebrush steppe ecosystems increases with drier and warmer conditions, the loss of perennial herbaceous understory, and the removal of big sagebrush, the presence of which can facilitate native bunchgrasses and prevent establishment of exotic forbs through resource competition (Prevey et al. 2010; Reisner 2010; Chambers et al. 2014). Although no pretreatment data were available for RGdNNM, it follows that if the perennial herbaceous understory was depleted before mechanical treatment, conversion to an invasive species–dominated vegetation community would be likely. Spatial data documenting treatments in RGdNNM were available from 2002, although treatments have been recorded in the area since the late 1950s (Pilliod and Welty 2019). It is possible that many more of the monitoring plots fell within historic vegetation treatments that produced long-lasting legacy effects and influenced the vegetation communities resolved by our analysis. However, recent, known mechanical treatments across ecological site types produced similar vegetation communities. When treatment type was included as a predictor along with environmental variables in the classification tree analysis, the first split of the tree separated plots in mechanical treatments from plots with no treatment, fire treatment, or chemical treatment. Treatment type was not used in the final classification tree analysis, as ecological site concepts are based on the abiotic variables that determine spatial heterogeneity rather than land use history.

Benefits of fuzzy clustering

By using fuzzy clustering to catalog the vegetation communities occurring within RGdNNM, and then examining how they were distributed across and constrained by abiotic variables, we were able to develop interpretable, ecologically consistent groups with a data-driven, inductive approach. Plant functional group compositions differed among ecological sites with variation in topography, climate, and soil properties. Plant functional group compositions and distributions of structural attributes also differed within ecological sites due to management and disturbance history, indicating differences in ecosystem processes and subsequently, alternative states (Allen-Diaz and Bartolome 1998; Stringham et al. 2003; Kachergis et al. 2012). We found fuzzy clustering to be a superior method when compared with hierarchical or crisp k-means clustering in its ability to represent noisy and complex vegetation data, consistent with other studies that have used fuzzy cluster analysis for plant community classification (Equihua 1990; Salski 2007; De Cáceres et al. 2010; Miller et al. 2011). As illustrated in Figure 3, many plots are represented by multiple colors (overlapping color wedges). This indicates the overlap among vegetation community clusters due to similarities in plant functional group compositions and/or the distributions of structural indicators. The use of hierarchical or crisp k-means clustering would have forced boundaries around these clusters and made the overlap invisible. Instead, we used the overlap to identify the plots that did not fit neatly into a cluster (the fuzziest plots). The use of fuzzy clustering required the qualitative classification of the fuzziest plots to create a grouping schema in which every plot was assigned to a modal concept. Plots with high membership values and clear identities were used to establish cohesive modal vegetation community concepts. The fuzziest plots were used to identify outliers to well-represented vegetation communities and additional vegetation communities that were underrepresented by the dataset. The exercise of landscape classification necessitates drawing boundaries around units and labeling them as discrete, although landscapes are typically composed of gradients and edge cases. Analytical routines that can maintain these gradients, such as fuzzy clustering, may be able to portray a more accurate picture of the landscape to inform management.

Management implications

Using relationships among plant communities, environmental attributes, and disturbance history, we have hypothesized ecological sites and vegetation states for the RGdNNM. Our landscape classification includes links between environmental heterogeneity and plant communities, as well as previous vegetation treatments where possible. The associations among environmental heterogeneity, vegetation communities, and previous vegetation treatments posited by our analysis can be used to inform ecological site development and guide management in RGdNNM. For example, quantitative summaries of plant communities can be compared with one another, as well as with existing quantitative structural and functional ecosystem thresholds, to estimate resilience, ecosystem service provisioning, and risk of degradation (Brown and MacLeod 2011; Oliver et al. 2015; Pellant et al. 2020). Three vegetation communities identified in the dataset have ranges of percent bare ground and large canopy gaps that fall below 20% and 35%, respectively (see Table 3); beyond these thresholds, fluvial erosion and aeolian sediment flux increase exponentially as bare ground and large canopy gaps increase (Webb et al. 2014). Land managers may seek to maintain or restore these vegetation communities on the landscape to reduce erosion potential. In addition to comparing quantitative indicator values with known quantitative thresholds, vegetation community concepts can be used to identify management thresholds (e.g., for which vegetation communities that are candidates for chemical shrub-removal treatment) (Anthony et al. 2021). Within ecological sites, plant communities can be compared with one another, as well as with existing quantitative structural, functional, and management thresholds, to estimate ecosystem service provisioning, vulnerability, and areas of concern (Suding et al. 2004; Monaco et al. 2012; Webb et al. 2020).

Both quantitative and qualitative data sources are essential to ESD development (Knapp et al. 2011a). The workflow presented here is not intended to produce ESDs and STMs with an automated, exclusively data-driven approach but rather is intended to use inductive reasoning and multivariate analyses to posit ecological site and state concepts that can be further deductively refined (see Fig. 1; Moseley et al. 2010). The core methods monitoring datasets include additional quantitative data, such as percent volume rock fragments by size class, soil aggregate stability ratings, and qualitative rangeland health observations, that were not used in our workflow. These data were not included in our analyses due to the incompatibility of their data structure with the multivariate methods selected to make best use of the variables listed in Table 1. However, quantitative indicators of soil erosion and qualitative rangeland health observations present an opportunity to better understand alternative states and should be explored as ecological site and state concepts are refined. Ecological site and state concepts in RGdNNM should be verified through workshops with land managers and users (e.g., Knapp et al. 2011b; Kachergis et al. 2013). Although quantitative data describing legacy land use and disturbance histories were not available for RGdNNM, local expertise could provide qualitative assessments of how these dynamics have shaped vegetation communities to determine transitions between states, as well as management thresholds for initiating restoration or change in land use. Ongoing data collection in RGdNNM in collaboration with the BLM and NRCS will serve to verify geographic boundaries of ecological sites and to solidify relationships between plant communities and ecological sites. Despite the large number of monitoring plots within RGdNNM, and across landscapes on a national scale, it can be assumed that monitoring data for any landscape may not represent all ecological sites and states, or the full ranges of variability in plant communities contained therein. Collecting additional data in ecosystems that are not spatially extensive will fill gaps where existing ecological sites and states may not have been adequately sampled by the monitoring dataset. Collecting pretreatment and post-treatment monitoring data within future vegetation treatments stratified by ecological sites will, over time, inform data-driven STMs.

The established workflow may be applied across other landscapes lacking ESDs where monitoring data exist. The workflow may also be applied where ESDs exist and need to be updated to reflect current states or establish quantitative ranges of biotic and abiotic indicators for use in developing monitoring benchmarks. An important consideration for using monitoring data to develop ecological site concepts is the trade-off between using monitoring data for model development and for assessment purposes. If the same data are used to generate ESD concepts and to subsequently evaluate the condition of resources against benchmarks taken from those ESDs, there is a risk of introducing circular reasoning into management decisions based on these evaluations (Webb et al. 2020). Conversely, deriving benchmarks solely from monitoring data (e.g., based on percentiles of indicator values), without ecological site and state context, presents the potential risk of assessing land health against shifting baselines (away from desired conditions) and/or managing for conditions that put sites at risk of undesirable state transitions (Soga and Gaston 2018; Webb et al. 2020). Thorough documentation of data use in developing benchmarks can help by providing date limits for the selection of monitoring data for evaluation purposes (e.g., evaluate with monitoring data collected after 2019). Depending on the size of the dataset, the data can be subset before developing ecological site and state concepts to reserve data points to be used exclusively for monitoring purposes (e.g., to evaluate grazing allotments for permit renewals). Including multiple lines of evidence in the development of benchmarks is important where data will be used for both establishing benchmark values and evaluating resource condition to validate that quantitatively derived benchmarks align with desired conditions.