Preparing ecological site and soil data

This study focused on ESDs from MLRA 35, 36, 34A, and 34B (USDA-NRCS 2006) and producing maps for roughly the areas of these MLRAs that fall within the Colorado River Basin. We queried ecological sites from a 2018 snapshot of the gridded Soil Survey Geographic Database (gSSURGO) (Soil Survey Staff 2018) in the study area and tabulated them by unique ecological site code concatenated to the associated ecological site name (see Fig. 1a). This resulted in > 1,200 entries, but many were not documented further than the record in gSSURGO or were duplicates due to slight typographical errors. We were able to link 826 unique entries to sites in the Ecosystem Dynamics Interpretation Tool (EDIT, Bestelmeyer et al. 2016) database ( https://edit.jornada.nmsu.edu/). Of these, further analysis included 1) 353 ESDs that had complete tabular data available documenting reference vegetation state communities, the most commonly available supporting data for ESDs; and 2) 252 ESDs that had published STMs. Experts characterize reference communities for each ESD by collecting and summarizing vegetation production using double sampling (Wilm et al. 1944) and other data in the field in areas deemed to represent reference conditions.

Many MLRAs have dichotomous keys for helping to identify ecological sites (USDA-NRCS 2019a). These keys reflect important structural information about how abiotic properties distinguish ESDs. We acquired keys for all study area MLRAs through email from NRCS staff (see Fig. 1a) and then encoded them into a digital format for analysis (explained later in ‘Defining soil geomorphic units’ section). The EDIT website also provides many ESD keys.

To quantify representative soil and geomorphic values for each ESD, we queried properties from gSSURGO for map unit component polygon delineations linked to each ESD. We weighted representative property estimates more to commonly correlated soils for each ESD by averaging values across all polygons ESDs are linked to. We summarized depth-specific soil properties (e.g., sand content) for surface soil as a depth-weighted average from 0 to 30 cm and for subsurface soil from 30 to 100 cm (or until bedrock). Soil and geomorphic properties examined from SSURGO include electrical conductivity (EC) in saturated paste to measure salinity; flooding frequency; sodium adsorption ratio; gypsum content (% weight in < 20 mm fraction); slope gradient (%); soil depth to bedrock contact (lithic, paralithic, or densic); volumetric rock content (%); clay content (% weight); and sand content (% weight). An important theme in ecological sites is delineating areas receiving extra moisture beyond ambient precipitation due to either runoff or subsurface moisture to define ‘Bottoms’ sites. We defined the extra moisture criteria as either an ESD that explicitly states that it receives extra moisture or if the SSURGO component flooding frequency was “rare” or more frequent.

Using reference state data from EDIT, we averaged reported production of specific important species, groups of species, and functional groups across all reference state communities for each ESD (Table 1). We found many errors in functional group attribution in the ESDs (e.g., obvious grass listed as a shrub), so we created a custom look-up table to re-correlate the species-level data to appropriate functional groups (  Table S1 (mmc1.xlsx), available online at …; n = 1 573 species). We also correlated species-level data to plant duration (e.g., perennial vs. annual) and photosynthetic pathway for all graminoids (  Table S2 (mmc2.xlsx), available online at …) following recent studies (Bruhl and Wilson 2007; Roalson et al. 2010; Osborne et al. 2014; Sage 2016). We defined vegetation parameters based on commonly important vegetation described in ESDs and from input in meetings with rangeland conservationists from the Bureau of Land Management (BLM) in western Colorado.

Figure 1.

This flowchart shows the variety of data sources and analyses employed for both creation and mapping of soil geomorphic units (SGUs) and ecological site groups (ESGs). a, Development started by preparing data from ecological site descriptions (ESDs) in the Ecosystem Dynamics Interpretative Tool (EDIT), ESD keys, National Soil Information System pedons, and the Soil Survey Geographic (SSURGO) database. We expand previous work (Duniway et al. 2016, b) to initially create an SGU key by optimizing reference production differentiation through evaluating different candidate keys using nonparametric multivariate analysis of variance (NPMANOVA) models (c). d, The developed key was used to classify National Soil Information System soil pedon field data for use in producing a digital soil map of SGUs. e, Final ESGs were determined by testing combinations of thresholds in gridded climate data (Aridity Index [AI], Max and Minimum Temperatures of hottest/coldest months [MAXT, MINT], and monsoon precipitation ratio [PPTRT]) combined with SGU classification of ESDs in SSURGO to then test for optimal breaks using a NPMANOVA model optimization process that considered both reference production data and documented ecological states from a generalized state-and-transition model table (f). g, The chosen combination of climate thresholds with SGU classes was then used in combining the SGU map with gridded climate data to render a final ESG map. h, Finally, rangeland monitoring vegetation cover data from the Bureau of Land Management Assessment Inventory and Monitoring and Natural Resources Conservation Service Natural Resource Inventory were overlaid on the maps of SGUs and ESGs to corroborate and assess the final classification using plots and NPMANOVA marginal effects models (i).

Table 1

Vegetation indicator variables used for subsequent analysis of both reference production and monitoring cover data. Species names are based on the US Department of Agriculture Plants database (USDA-NRCS 2019). If multiple species are listed, that category includes a sum of reference production from all species or “any-hit” summaries (% of transect points where any of those species were observed) for cover data. Short names are used heretofore in text and figures for brevity.

We tabulated possible ecological states for each ESD from published STMs in EDIT (see Fig. 1a). We categorized states first at a higher-level grouping as reference, herbaceous invasion, woody encroachment, eroded, bare ground, or perennial species loss. Where possible, we further distinguished more detailed levels of these groups (e.g., tamarisk vs. greasewood woody encroachment). This resulted in a binary table recording, of which 39 different types of unique ecological states occurred in each ESD (  Table S3 (mmc3.xlsx), available online at …). We chose not to explicitly include transition drivers (e.g., grazing, fire, drought, management actions) in analysis because of recent concerns about validity of ESD transition drivers due to the lack of longitudinal studies available needed to establish causality, among other concerns (Twidwell et al. 2013; Bestelmeyer 2015). Although some important processes (e.g., disturbance) that result in alternative states were not explicitly considered, the inclusion of alternative states that are characterized by altered processes (e.g., accelerated erosion) in our STM state data do account for these influences in the analysis to some extent. Ecological states essentially represent site behavior as varying propensities of different sites to manifest certain states (e.g., dune states in sandy sites), even if causal processes of an alternative state might not be fully characterized.

Assessing multivariate ecological responses with nonparametric multivariate analysis of variance

We used nonparametric multivariate analysis of variance (NPMANOVA) tests as implemented in the vegan R package (Oksanen et al. 2007) at various points in the study to understand how a set of responses vary with one or several independent variables. NPMANOVA tests can analyze multivariate response data with non-normal and binary distributions using the Bray-Curtis distance metric (Anderson 2001, 2017; Anderson et al. 2008; Anderson and Walsh 2013). For this study NPMANOVA assessed how well abiotic independent variables (soils, geomorphology, and climate) explained variation in vegetation and ecological state response variables. For models with one independent variable, we used the “adonis” function to run NPMANOVA tests. Since variable effects are evaluated sequentially in NPMANOVA with multiple independent variables, we used the “adonis2” function for models with more than one independent variable. With “adonis2,” repeated model building allows a variable of interest to be evaluated for the iteration where it is added last to the model (i.e., marginal effect). We ran 999 model permutations to calculate all reported P values. NPMANOVA test results (mainly in P values) reflect both difference and dispersion effects for unbalanced data. In most of the model runs presented here, models had weak to moderate (P values 0.1–0.01) dispersion effects, but dispersal box-plots showed such dispersion effects in mostly one or two classes. Although our goal was to maximize differences between groups, differing amounts of variation within groups (dispersal) may indicate differing ecological behavior as this may distinguish ecological units with differing diversity. However, all models used for interpretation were also examined with bivariate plots to verify that differences between groups were the dominant drivers in models.

Defining soil geomorphic units

We adopted groupings from recent regional expert ecological site generalization workshops as a priori “strawman” concepts (Bestelmeyer et al. 2016; Duniway et al. 2016) and initially classified all ESDs in our study to match these groups (see Fig. 1b). These workshops focused on grouping ESDs with similar STMs by intuitive soil and geomorphic concepts that did not explicitly address climate gradients, and we refer to this type of thematic grouping hereafter as “soil geomorphic units” (SGUs). Our work began by classifying ESDs into previously defined strawman SGUs using the published ESDs to qualitatively assign an SGU, but during this process we decided that the addition of new SGUs, adjustments to previous concepts, and explicit definitions were necessary. Our initial strawman SGU classification was then evaluated using MLRA ESD keys acquired from ESD specialists and reference community production data (see Table 1). Keys were translated into data.tree objects (Glur 2019) in the R statistical language (R Core Development Team 2017) for analysis of how frequently variables were used to split ESDs at various hierarchical levels. We assumed that variables used at higher hierarchical levels of keys and/or more frequently were more important for SGU differentiation. With this approach we modeled our own SGU dichotomous key off of the structure of existing keys to better translate our qualitative strawman classes to quantitative dichotomous definitions consistent with regional expertise. After translation of SGU definitions to a key format, we ran boxplot evaluations and a series of NPMANOVA tests to maximize variation in reference state production explained by the SGUs (see Fig. 1c). Tests included adjusting SGU definitions focusing on 1) potential inclusion of newly identified SGUs and 2) adjustments in dichotomous key structure and soil/geomorphic value thresholds. We chose the SGU key resulting in the NPMANOVA model representing the highest percentage of variation explained among our trials.

Mapping soil geomorphic units using digital soil mapping

We based the mapping of SGUs on a digital soil mapping (DSM) workflow (McBratney et al. 2003) relating raster datasets representing soil forming factors to a set of field soil descriptions classified to SGU to create an interpolative prediction model (see Fig. 1d). Model training sites (n = 12,320, see Fig. 1a and 1d) came from a 2019 snapshot of the internal USDA-NRCS National Soil Information System (NASIS). Raster covariates included 59 layers (Table A.1) derived from the 30-m National Elevation Dataset (Gesch et al. 2002), 2011 National Land Cover Dataset (Homer et al. 2015), a terrestrial vegetation inventory (Gergely and McKerrow 2013), regional landforms (Iwahashi et al. 2018), gamma radiometric layers (Hill et al. 2009), and Landsat 5 TM and Landsat 8 Collection 1 Tier 1 top of atmosphere reflectance products prepared in Google Earth Engine (Gorelick et al. 2017) using calibrations from Chandler et al. (2009). We calculated terrain metrics (see Table A.1) using tools in SAGA GIS (Conrad and Wichmann 2011) and ArcGIS (ESRI 2014). We produced geomorphons to represent local landforms (Jasiewicz and Stepinski 2013). Flow accumulation routing for relevant terrain layers followed the “Dinf” approach (Tarboton 1997). We reprojected all non-terrain-based rasters to the ∼30-m resolution of the elevation grid using bilinear resampling. We performed all raster modeling with the R raster package (Hijmans et al. 2016).

Riparian and rock outcrop ESDs are not consistently available, and we did not analyze them but recognized their importance for future work. We mapped both as separate ESGs using pseudotraining points identified in Google Earth Pro (40 points in each group for each MLRA picked to represent the diversity of both groups). Rock outcrops points were chosen to be pixels with at least ∼75% bedrock outcrop.

We classified each soil training site (excluding riparian and outcrop pseudotraining points) to a nearest soil series and then linked to the nearest SSURGO soil component in its 1) home polygon or 2) an adjacent polygon following a workflow from several recent DSM studies to query soils property estimates from SSURGO for a location (Maynard et al. 2019; Nauman and Duniway 2020). We then classified these sites to an SGU using soil parameters and the correlated ESD from SSURGO. We built a random forest model using the randomForest package (Breiman 2001; Liaw and Wiener 2002) to predict SGUs for all 30-m pixels in the study area. Since our SGU dataset had an unbalanced class distribution, we used a synthetic minority oversampling technique shown to help ameliorate problems with minority class accuracy in tree-based models (Branco et al. 2016; Taghizadeh-Mehrjardi et al. 2019). We iteratively adjusted custom weights for all classes to improve expert visual evaluation and class-wise validation accuracy. We evaluated performance with class-wise user and producer accuracies, overall accuracy, and the kappa coefficient (Cohen 1960) calculated on a randomly withheld 20% of training instances for validation. We also mapped the probability of the most probable class and Scaled Shannon's Entropy (Shannon 1948; Maynard et al. 2019) to provide relative metrics of uncertainty for users. We limited mapping to below 2,750 m in elevation due to a lack of consistent data above that elevation for both soil pedons and ESDs.

Preparation of climate data

Raster climate data included aridity index (AI), maximum temperature of the warmest month (MAXT), minimum temperature of the coldest month (MINT), and a monsoon precipitation ratio (PPTRT) of summer (June–September) to annual precipitation over 30-yr normalized datasets (see Fig. 1e). We downloaded the AI data from the second version of a global AI compilation that uses WORLDCLIM 2 as inputs (Trabucco and Zomer 2018). AI is calculated as mean annual precipitation divided by mean annual potential evapotranspiration as derived by the precipitation and temperature-based Hargreaves model (Hargreaves and Allen 2003). We downloaded MINT and MAXT maps from the 1-km resolution 1970–2000 WorldClim 2 Bioclim dataset (Fick and Hijmans 2017). We derived the PPTRT map from the 1981–2010 800-m resolution PRISM dataset (PRISM Climate Group 2010). To enable representative climate characterization of each ESD, we averaged climate datasets for all gSSURGO polygon delineations with components linked to each ESD in the study area. We reprojected climate layers to the ∼30-m elevation grid using bilinear resampling only for the final ESG map rendering, but users should view derived climate zone boundaries as approximate.

Combining climate zones and soil geomorphic units into ecological site groups

Final ESG definitions combined climate thresholds with SGUs in three broad steps: 1) identification of climate threshold candidates, 2) testing SGU-climate groupings using reference vegetation production data in NPMANOVA models to narrow down options, and 3) final ESG determinations using STM data in NPMANOVA models in addition to reference production data models for final ranking (see Fig. 1f). An information theoretic approach ranked NPMANOVA models built with differing combinations of climate thresholds and SGUs to identify the best combinations for final ESG definitions.

Identifying climate thresholds

We created candidate climate thresholds using a combination of quantile breaks and AI values equivalent to common PRISM-based mean annual precipitation (MAP) cutoffs documented in ESD keys. We calculated equivalent AI values to the MAP breaks by querying the mean of all AI pixels that overlaid PRISM MAP pixels within ± 0.01″ of the 9,″ 13,″ and 16″ breaks from the Utah MLRA 35 ESD key, resulting in AI values of 0.1260, 0.1749, and 0.2179 to include as candidate thresholds. We added these values to additional quantile-based breaks identified to represent other parts of the AI distribution (i.e., 7.5%, 20%, 40%, 50%, 60%, and 92.5% quantiles). We used slightly different quantile breaks, 7.5%, 20%, 35%, 50%, 65%, 80%, and 92.5%, for the other climate parameters since no NRCS-based equivalents existed (Fig. 2).

Akaike's information criterion ranking

To evaluate combinations of climate breaks with SGUs, a small sample size implementation of Akaike's information criterion (AICc) ranked potential SGU-climate groupings in NPMANOVA models. AICc incorporates sum of squares statistics in order to balance fit and complexity for model selection (Burnham and Anderson 2001; Anderson et al. 2008). To better examine relative model strengths, we also calculated Akaike weights and relative likelihood ratios (Burnham and Anderson 2004) from the AICc distributions using the qpcR package (Ritz and Spiess 2008).

Testing SGU-climate combinations with reference production data

In combining climate breaks with SGUs to create ESGs, we did not assume that all SGUs should be kept separate across the regional climate gradient. Instead, we tested two grouping combination strategies that allowed for both 1) testing possible SGU combinations nested within climate zones (SGU groupings could be different in different climate zones) and 2) simultaneously testing SGU aggregation schemes and climate zonation breaks (SGU groupings the same across climate zones). The rationale with simultaneous grouping is that the effects of SGUs and climate are independent, whereas the nested scheme reflects potential interactions. In trials, we considered up to four AI breaks and up to three breaks for the other climate parameters in addition to up to four SGU aggregations in all combinations that resulted in fewer than 60 final classes. Nested combinations were evaluated first to help constrain potential SGU aggregations allowed in simultaneous groupings to ease computational burden.

To narrow down the large number of possible nested SGU-climate combinations, we implemented an initial reference data NPMANOVA AICc ranking of climate breaks without any SGU aggregation. This resulted in four sets of climate zones identified out of 3,145 unique trials based on Akaike weights. Because there were still > 1 x 10¹⁰ potential combinations when all SGU combinations were considered within each of the identified potential sets of climate breaks, we also chose an initial set of allowed SGU pairing combinations based on our judgment of similarity (25 pairs documented at  https://github.com/usgs/Predictive-Soil-Mapping/blob/master/ESGs/ClimateOptimization_v2b_nested.R commit 06ea86a, line 670).

The four climate zone schemes identified were then used to further narrow down possible combinations with SGUs in nested configurations. Out of the four climate zonation schemes identified for nested SGU trials, the one with the lowest AICc had only two climate zones, while the others had three zones. Since potential SGU aggregation schemes within three climate zones still numbered in the billions, we made a pragmatic assumption that the SGU combination pairs identified in the best two-zone nested models were more likely to be the best candidates for three zone models. Thus, we further narrowed down SGU combination pairs based on two-zone nested climate-SGU models by selecting only SGU pairings that were identified in at least 5% of NPMANOVA models with AICc relative likelihoods > 0.5. We limited all further trials to these potential SGU aggregation possibilities to keep numbers of possible combinations to under 1 x 10⁸ for practical computational reasons, and because the SGU pairs identified seemed reasonable. We analyzed reference production NPMANOVA trials for the other three nested climate configurations using AICc, Akaike weights, and relative likelihoods to find the best combinations possible within each set of climate zones. We retained trials for final STM state occurrence analysis based on picking combinations from NPMANOVA reference production models within each nested climate zone scheme with relative likelihoods > 0.5 due to a lack of distinctively better Akaike weights among best-performing models.

Figure 2.

Plots showing distribution of climate parameters means for the 353 ecological site descriptions used for optimizing climate zones for ecological site grouping. The red triangles indicate candidate climate thresholds to be tested as selected by quantile values and aridity index equivalents to representative NRCS precipitation breaks. Candidate thresholds were used to define climate zones boundaries to then combine with soil geomorphic units to create final groups.

For simultaneous SGU-climate grouping trials, we opted to use the limited set of potential SGU aggregations as determined from the two-zone nested climate configuration models to help constrain grouping possibilities to ease computation burden, which resulted in ∼2.7 million trial possibilities. The AICc analysis showed that top models were distinctly better than others among the simultaneous model trials, so models with relative likelihoods > 0.1 were retained to narrow down possibilities for final evaluation with STM state occurrence data.

Testing SGU-climate combinations with ecological state data

For final ESG model selection (climate-SGU aggregation scheme), we examined NPMANOVA model trials using binary state existence data from STM models for each underlying ESD (see   Table S3 (mmc3.xlsx)) as response variables to further evaluate retained candidate climate-SGU combinations from the initial reference production trials. We then separately assessed all retained climate-SGU combinations between STM state occurrence and reference production response datasets using AICc ranking of NPMANOVA models. We selected the final climate-SGU schema for ESG definition to have the lowest average AICc rank between both the set of STM state occurrence NPMANOVA models and the set of reference production data NPMANOVA models (see Fig. 1f and 1g).

Corroborating concepts with independent monitoring data

To further corroborate ESG concepts, we compiled line point intercept vegetation cover data from the NRCS National Resources Inventory (NRI) program (Nusser and Goebel 1997; Herrick et al. 2010) and the BLM Assessment Inventory and Monitoring Landscape Management Framework (AIM-LMF) monitoring program (MacKinnon et al. 2011; Toevs et al. 2011, see Fig. 1h). These programs were designed around consistent protocols (Herrick et al. 2017) and provide broad statistical inference over private lands (NRI) and BLM lands (AIM-LMF). There were 2,155 AIM-LMF plots measured from 2012 to 2018 and 1,569 NRI plots measured from 2004 to 2016 in our study area. All plot locations were sampled once. We calculated “any-hit” plot vegetation cover (i.e., a species or functional group present at any level in the canopy is counted) using the groupings in Table 1. We linked monitoring data to soil geomorphic units and ecological site groups by intersecting plot locations with predictive maps.

We generated bivariate plots and NPMANOVA models to examine relative strengths of climate variables, SGU classes, and ESG classes for explaining the variation in the three different vegetation datasets—reference production data by functional groups, STM state occurrence, and vegetation cover monitoring data (see Fig. 1i). To understand broader controls on the system, we compared total variation explained in NPMANOVA models with only SGUs, only ESGs, SGUs plus climate parameters, and SGUs plus climate parameters and climate-SGU interactions. We also examined marginal variable effects in NPMANOVA models including SGU and climate variables to look at relative influences of SGUs and climate variables. Data generated during this study are available from the USGS ScienceBase-Catalog (Nauman 2021).

Soil geomorphic unit definitions

Development of SGUs resulted in 16 classes defined in a dichotomous key (Fig. 3) hierarchically separating sites by 1) areas with moisture beyond ambient precipitation; 2) salinity; 3) depth, slope, and rock content; and 4) soil texture. These classes included several new units not documented in previous grouping efforts: Very Shallow, Loamy Uplands, Clay Uplands, and Riparian. The SGU key structure was mainly based on a pooled analysis of all MLRA ESD keys by hierarchical level showing that annual precipitation, soil depth (to bedrock or restriction), run-on versus runoff geomorphic positions, and salinity were the most common splits in the higher levels (first 1–3 or 1–4 levels in each key depending on key size; keys had 9–14 total levels). All MLRA keys had higher-level climate groupings mostly defined by annual precipitation, but sometimes by elevation or important geographic features that compartmentalize climate variation. After precipitation, MLRA keys generally split out areas with moisture available beyond ambient precipitation due to run-on or shallow water tables separating bottoms and riparian sites from uplands areas. After extra moisture, sites were generally split by salinity or soil depth. At lower levels, rules often included soil texture, slope gradient, landform, and dominant vegetation. Since current vegetation often reflects land use history, we focused on the abiotic concepts of previous keys for formalizing groups.

Figure 3.

Dichotomous key diagram defining soil geomorphic units (SGUs) starting from top. When only one branch from a node is defined, it is assumed that instances in the other branch of that node include remaining sites that do not meet the criteria in the opposing branch. Nodes with no box represent intermediate pools of sites that need further rules to differentiate to an SGU. Note that there are two paths for a site to be classified as Saline Hills by either having a high sodium adsorption ratio (SAR) or high electrical conductivity (EC) if there is not significant gypsum and low SAR. EC is based on the saturated paste method.

Soil descriptions in ESD documents were not always consistent with properties queried for soil components in SSURGO linked to the ESDs. This was particularly noticeable for the important variables of soil depth, salinity, and texture. This is not surprising given that both the previous ecological site groupings (Duniway et al. 2016) and the actual ESDs often have somewhat qualitative and overlapping (i.e., not mutually exclusive) abiotic definitions. Iterative testing of SGU key structure and thresholds resulted in explained variation of reference community production rising from 29% to 34% in NPMANOVA models. Table A.2 provides more description of each SGU and general ecological community, and Figure 3 displays the final dichotomous key.

Map of soil geomorphic units

The final predictive map (Fig. 4) of the SGUs had an overall validation accuracy of 57.4% and Kappa coefficient of 0.53 for the first most probable class with overall accuracy jumping to 73.6% and 82.3% if second and third most probable class matches were considered, respectively. Many of the misclassifications of the most probable class were among similar SGUs (e.g., Shallow as Very Shallow, Table 2). Producer's accuracies (PAs) averaged 0.52 (standard deviation [SD] = 0.15) and User's Accuracies (UAs) averaged 0.57 (SD = 0.18) with most ranging from 0.45 to 0.7. Sandy Bottoms had notably low-class accuracy with a UA of 0.21 and PA of 0.15. This may be due to a lack of training data or to these units occurring in small areas in rugged drainages where they may not always cover a full pixel. Sandy Bottoms tend to be confused with Sandy Uplands, Riparian, and Saline Bottoms SGUs in validation. As many landscapes in the region are quite heterogeneous, a portion of misclassifications are also likely due to mixed pixels, where multiple SGUs may be present in one 30-m pixel.

Figure 4.

a, Predicted soil geomorphic unit map of study area showing general distribution of dominant units in the southwestern United States and b, a closer look at the area near Hanksville, Utah, showing contrasts among shallow sandstone-derived soils to the east, saline soils derived from Mancos Shale to the northwest, and rocky soils coming off the Henry Mountains to the southwest.

The SGU map shows large areas of shallow soils (pinker areas, see Fig. 4a) common in the Upper Colorado River Basin, along with distinctive saline soils associated with various regional marine shale and evaporite formations (bluer areas, see Fig. 4). Other distinct features include the sandy soils associated with the weathering and subsequent transport of sediments sourced from the famous scenic Triassic and Jurassic sandstones of southern Utah and northern Arizona (yellow areas in Fig. 4a). In closer visualization of the Hanksville, Utah area, the SGU map reflects landform and geologic details, with rocky colluvial and alluvial soils coming off the Henry Mountains in the southwest, shallow and sandy soils in the canyons of the Dirty Devil River system on the east, and saline Mancos Shale soils of the Factory Butte area in the northwest (see Fig. 4b). All maps are available online (Nauman 2021 at  https://doi.org/10.5066/P92OPRMV).

Combining Climate Zones and SGUs into ESGs

Reference production data

Models with climate-nested SGU aggregation. We identified four candidate climate zonation schemes out of 3,145 possibilities to use in nested SGU aggregation trials as part of final ESG determination. Analysis and assessment of the top 10 models identified with AICc revealed several common climate breaks. The top two models were the most informative with larger jumps in Akaike weights than subsequent models (Table 3). AICc weights did not increase substantially among the models ranked third and higher, but examination of other models in the top 10 led us to include the third and fourth best models for further testing due to moderate AICc weights and contrasting climate breaks. All top 10 models included an AI break, with most including an AI of 0.144, followed by three that included a break of 0.104 and three that included an AI of 0.134. Only three of the top 10 models included a MINT break (all –17.07 °C), four models included a MAXT break (25.02 °C for three and 28.17 °C for one), and none of the models included a monsoon precipitation ratio break. The top model included just one AI break of 0.144; the second-best model included AI breaks of 0.104 and 0.144. The third and fourth best models had AI breaks of 0.144 in addition to a MINT break of –17.07 °C and a MAXT break of 25.02 °C, respectively.

Analysis of the two-zone climate scheme (top model from the 3 145 previous zonation trials) in nested trials identified 10 SGU pairs to constrain potential SGU lumping in further trials with the other more complex climate zonation schemes (pairs present in at least 5% of the tested two-zone models with relative AICc likelihoods > 0.5; n = 561 models). This narrowed the potential SGU pairings in later analysis to Breaks-Gypsum, Breaks-Very Shallow, Gypsum-Very Shallow, Saline Uplands-Loamy Uplands, Shallow-Deep Rocky, Sandy Uplands-Loamy Uplands, Loamy Uplands-Finer Uplands, Finer Uplands-Clay Uplands, Sandy Bottoms-Saline Bottoms, and Saline Bottoms-Bottoms. All nested NPMANOVA trials exhibited a similar lack of distinction in model AICc values with no clearly better candidates. So, we narrowed down nested trials from each climate scheme to the set with relative likelihoods > 0.5. We retained a total of 13,151 potential models from nested trials for testing with the more limited STM state occurrence data.

Table 2

Confusion matrix with Producer's and User's accuracy rates for each SGU class based on the 20% withheld validation set of NRCS soil pedons. Boxes on the diagonal represent the counts of pedons correctly classified.

Table 3

Climate break model performance and AICc criteria for top 10 NPMANOVA trials on reference production data including Akaike weights (AICcWt) and Akaike relative likelihood (AICcrelLL).

Incorporating STM state models for final ESG determination

In final ESG analysis, the AICc rankings varied considerably between the STM state occurrence and reference-based NPMANOVAs (Fig. 5). Trials for nested climate schemes 2 and 4 along with the simultaneous grouping trials had the best ranks for the reference production NPMANOVA models, whereas nested approaches 1, 3, and 4 were top ranked with the STM NPMANOVAs. Two models tied for the lowest average rank between reference production and STM models with both averaging a rank of 337.5 out of the 13,170 models (see Fig. 5). The next closest model had an average rank of 365.5. To break the tie, we chose the model with the higher average relative likelihood between the response sets. The final chosen model ranked 478th among all reference production models and 197th among the STM state occurrence models.

Figure 5.

Akaike's information criterion (AICc) rank comparisons for soil geomorphic units–climate combinations in NPMANOVA models using a, reference production data and b, state-and-transition model state binary data using the boxplot base function in the R statistical language. Model selection involved minimizing the AICc rank average between reference production and STM-based models. Ranking is done within all modeling approaches for each set of response variables (i.e., models in part a are ranked separately from models in part b).

Description of ESGs and corroboration with monitoring data

The final selected ESG classification (i.e., best NPMANOVA results) included three climate zones and nested SGU aggregation that varied among the different zones (Table 4). The first zone includes areas with AI < 0.144 and maximum temperature of the warmest month (MAXT) > 25.02 C, which we labeled as “Arid Warm.” The second zone has AI > 0.144 and MAXT > 25.02 C, which we labeled as “Semiarid Warm.” The third zone has AI > 0.144 and MAXT < 25.02 C, which we labeled as “Semiarid Cool” (see Fig. 4). Production and cover generally increase in analogous ESGs from Arid Warm to Semiarid Warm to Semiarid Cool (see Tables 4 and A.3) but are also far higher in the ESGs that include hydrologically enhanced SGUs (Riparian, Bottoms, Saline Bottoms, and Sandy Bottoms). Sandy Uplands and Loamy Uplands were aggregated in all climate zones. We combined Saline Bottoms and Bottoms, as well as Finer Uplands and Clay Uplands in the Arid Warm climate zone. In the Semiarid Cool zone, we combined Saline Uplands, Sandy Uplands, Loamy Uplands, and Finer Uplands into one ESG. Additionally, in the Semiarid Warm zone, we combined Sandy Bottoms with Bottoms and Shallow with Deep Rocky. Overall, we documented 29 ESGs with the reference production data, excluding Riparian and Outcrops. However, in mapping overlays, there were 35 groups due to some SGUs extending into the Semiarid Cools zone outside of ESD climate averages. This disparity in ESG classes between the mapped areas and the ESDs is probably due to climate variability not captured when summarizing the ESDs to mean map unit values, but it may also indicate areas that still need further characterization. However, cover data available from BLM, NRCS, National Park Service (NPS), and high values from relevant ESD production tables could supplement this analysis.

Tree production is limited in the Arid Warm groups and mostly occurs in the Very Shallow, Shallow, Deep Rocky, and Breaks SGUs. Cover values from the monitoring data corroborate that tree cover is highest in Shallow, Very Shallow, Deep Rocky, and Breaks in the Warm climate zones. Cover data also reveal that tree cover is lower in shallow soils in the Semiarid Cool zone but still high in Breaks (see Table A.3). Cover data also show elevated tree cover in Arid Warm Sandy Bottoms in contrast to the reference data, suggesting that the random forest model incorrectly mapped some Riparian areas (high tree cover) as Sandy Bottoms. Pinyon and juniper species are similar in their affinity for SGUs that are rocky, shallow, and coarser in texture. Juniper has higher production than pinyon in the Arid Warm zone, both are similar in the Semiarid Warm zone, and both have lower production in the Semiarid Cool zone in the SGUs with available reference data. Ponderosa has its highest production in the Semiarid Cool Shallow ESG but also has small production values for some other ESGs. Ponderosa is also found in limited areas of Finer Uplands, Loamy Uplands, Deep Rocky, and Shallow SGUs in both the Semiarid Warm and Semiarid Cool climate zones. Aspen production is only recorded in the Semiarid Cool Bottoms. Cover data also show aspens in Semiarid Cool ESGs including Shallow, Deep Rocky, and the Saline-Sandy-Loamy-Finer Uplands ESG. Cover data are similar to reference production results for pinyon and juniper.

Shrubs are significant plant community components throughout the climate zones but tend to be more dominant in the Semiarid Cool climate zone and in ESGs that include Breaks, Deep Rocky, Shallow, and Very Shallow SGUs. Shrub cover in the Semiarid Cool ESGs is also pointedly higher than the other climate zones (see Table A.3). Saltbush species tended to be most prolific in the saline SGUs and in finer textured soils in the Arid Warm climate zone. Greasewood only has significant production in the various bottoms run-on SGUs, particularly in Saline Bottoms (Table A.4). Most sagebrush subspecies are not very productive in the Arid Warm zone. Basin sage is only a significant producer in the Sandy Bottoms and Bottoms of the Semiarid Warm zone. Mountain sage has its highest production in Semiarid Cool Deep Rocky or the Semiarid Warm Sandy and Loamy Uplands. Big sage is in nonsaline deeper soils of the Semiarid Cool and Semiarid Warm climate zones with particularly high production in Finer Uplands and the Semiarid Cool Deep Rocky ESG. Blackbrush production is mostly constrained to the Arid Warm zone with preferences for Shallow, Very Shallow, and Deep Rocky SGUs but with some production on Breaks, Gypsum, and Sandy/Loamy Uplands. Monitoring plot cover data are generally consistent for overall shrub-related parameters but differ slightly in how various big sage subspecies are distributed. Basin sage cover is more prominent in Loamy Uplands, Finer Uplands, Breaks, and Deep Rocky SGUs in the Semiarid Cool zone, while patterns in big sage cover data share similarity to both the basin sage and big sage reference production. This suggests that basin sage is not being identified well in the field relative to Wyoming big sage (much of the data are not identified to subspecies). Mountain sage cover seems to reflect reference production data across all big sage subspecies including cover in several of the Arid Warm ESGs, which also suggests some misidentification in the field.

Table 4

Mean reference production (kg/ha) by important vegetation indicators for each ecological site group (row). Each column is heat mapped (lightest to darkest) to highlight the ecological site group with the highest (darkest) production for each category. Vegetation categories follow short names in Table 1.

Grass production and cover vary substantially by both climate zone and SGU. Grass production dominates ESGs with Bottoms and Saline Bottoms SGUs but also in groups with deeper and less rocky soils. Alkali sacaton is a key producer in ESGs with Saline Bottoms and Bottoms SGUs. In general, C4 perennial grasses produce more in Arid Warm ESGs relative to C3 perennial grasses. The C3 and C4 perennial grasses produce more evenly in Semiarid Warm ESGs, and C3 perennial grasses dominate in Semiarid Cool ESGs. Grass cover from the monitoring plot data is relatively consistent with production trends except in Arid and Semiarid Warm Bottoms ESGs. In these areas, grass cover from field plots is much lower than seen in patterns in reference production data. Clay Uplands have much higher grass cover in all climate zones relative to trends in the reference data, but particularly so for the Semiarid Cool Clay Uplands ESG, which only has two supporting ESDs. In the monitoring plot data, alkali sacaton is also almost completely constrained to the Arid Warm and Semiarid Warm Saline Bottoms and Bottoms ESGs, in contrast to ESD reference production data, which shows it present throughout the climate zones.

Patterns in STM state occurrence within ESGs and SGUs differed in more subtle ways. All SGUs had an herbaceous invasion state of some kind (  Table S4 (mmc4.xlsx)), but ESGs in the Semiarid Cool climate zone had less overall relative occurrence of invaded states (  Table S5 (mmc5.xlsx)). Bromus spp. and Salsola spp. were the most common invaders but became less common in the Semiarid Warm and Semiarid Cool climate zones. Wheatgrass species (e.g., Pascopyrum smithii [Rydb.] Á. Löve and Elymus lanceolatus [Scribn. & J.G. Sm.] Gould) were listed as common invaders in Semiarid Warm areas. All SGUs also have some kind of woody encroachment, particularly of shrubs, but Breaks had much less documentation of encroachment. Greasewood, tamarisk (Tamarix L. spp.), and Russian olive (Elaeagnus angustifolia L.) encroachment are generally limited to ESGs with a bottoms SGU. Snakeweed (Gutierrezia spp.) and rabbitbrush (Chrysothamnus L. spp. and Ericameria nauseosa [Pall. ex Pursh] G.L. Nesom & Baird) encroach in all SGUs, but snake-weed is not explicitly documented in Semiarid Cool ESGs. Pinyon is primarily recorded as encroaching in Semiarid Warm ESGs of Clay Uplands, Finer Uplands, and Shallow/Deep Rocky. Juniper encroachment was broader, including the Arid Warm ESGs of Deep Rocky and Sandy/Loamy Uplands and even more so in Semiarid Warm ESGs including Clay Uplands, Finer Uplands, Sandy/Loamy Uplands, Shallow/Deep Rocky, and Very Shallow. Perennial species loss states were also ubiquitous across ESGs excepting those with the Gypsum SGU. The Semiarid Cool climate zone also tended to have higher relative occurrence of perennial species loss states. Loss of a grass species is the most common perennial loss with C3 species recorded as lost more often than C4. Breaks, Deep Rocky, Shallow, and Very Shallow SGUs had less propensity for perennial species loss states than other SGUs. Big sage subspecies loss is not recorded in Arid Warm ESGs and concentrated mostly in the Semiarid Cool ESGs, particularly in Clay Uplands, but is also recorded in ESGs with Deep Rocky, Sandy Uplands, Saline Uplands, Loamy Uplands, Finer Uplands, Shallow, and Very Shallow SGUs.

Erosion and erosion-prone states (with elevated bare ground) follow more specific patterns than other documented alternative states. Aeolian mobilization states (Dunes and Coppicing) are specific to Sandy and Loamy Uplands in the Arid Warm and Semiarid Warm climate zones. Sandy Uplands also have some ESDs where dunes are listed as a potential reference state (mostly in Arid Warm zone). Sites with rill or gully erosion states are concentrated primarily in Bottoms and Saline Bottoms SGUs and in the Arid Warm and Semiarid Warm climate zones but are also documented in Finer Uplands, Saline Uplands, and Sandy Uplands SGUs. More general “eroded” states are also documented for Bottoms, Saline Bottoms, Finer Uplands, Clay Uplands, Saline Hills, Shallow, Very Shallow, and Sandy Uplands SGUs. Elevated bare ground states are least frequently observed in Gypsum, Breaks, and Saline Bottoms SGUs and are most common in Finer Uplands, Clay Uplands, Sandy Bottoms, and Sandy Uplands across climate zones. Of the documented erosion-related states, all are in the Arid Warm or Semiarid Warm climate zone with the exception of one Clay Uplands ESD in the Semiarid Cool zone.

Relative compositional patterns are similar between reference production and cover datasets across ESGs with some exceptions (e.g., Fig. 6). More general indicators like total foliar cover had more differences with analogous metrics in total production (see Figs. 6 a and 6b), whereas more specific indicators like tree production showed more congruent trends to tree cover among ESGs (see Fig. 6c and 6d). Total foliar cover and production showed similar ESGs with lower values (e.g., Saline Hills, Gypsum, Very Shallow in Fig. 6) but diverged in trends for more productive groups, particularly for ESGs with various bottoms SGUs.

A soil-geomorphic framework for organizing ecosystem behavior

Distribution of the new ESGs is closely tied to bedrock geology; anticlinal and faulting features; parent material (e.g., eolian, alluvial, residual); and relative topographic position (Fig. 8), reflecting patterns in similar studies (McAuliffe 1994; Bestelmeyer et al. 2011; Michaud et al. 2013). In mapped SGUs, lower drainages like the Dolores and San Miguel Rivers in western Colorado provide clear examples of the more productive riparian and bottoms units. In the illustrated area, the bottoms tend to be saline due to the prevalence of evaporite and marine shale geologic formations like the Mancos Shale, a highly erodible and salt rich geologic formation that generally forms hilly badlands in drier areas (Cadaret et al. 2016; McGwire et al. 2020). A number of saline shales are exposed regionally (Fillmore 2010; Morrison et al. 2012), and studies have tied them to surface water salinity loading (e.g., Miller et al. 2017; Nauman et al. 2019) and dust sources due to sand bombardment from upwind sandstone structural benches (Flagg et al. 2014; Nauman et al. 2018). The saline and usually finer-textured soils derived from shales (blue classes in Figs. 3, 6, 7) have profound implications for the structure and composition of local ecosystems with low productivity and salt-adapted species (see Tables 4 and A.3). Both Saline Bottoms and Bottoms are also unique in their susceptibility to rill and gully water erosion (see   Table S4 (mmc4.xlsx)), which has been reported by others for the Colorado Plateau region as mostly related to Holocene and modern climate change but also to land use (e.g., grazing) in more limited cases (e.g., Bailey 1935; Balling Jr and Wells 1990; McFadden and McAuliffe 1997). However, our review of STMs in this study supports linkages between grazing impacts on water erosion in saline SGUs (e.g., R034BY101UT: Desert Alkali Bench ESD) similar to what others have reported in more recent studies (Duniway et al. 2018; Fick et al. 2020).

Shallow and sandier eolian soils tend to be associated with canyons that dissect structural bench and mesa systems dominated by sandstones (e.g., see Fig. 8). On top of mesas and structural benches, deeper soils occur that vary in texture with elevation, with Sandy Uplands and Loamy Uplands in the drier areas and Finer Uplands in wetter areas or on benches with finer-grained source lithology. Local anticlinal collapses and faulting are also associated with larger valleys that tend to have mixes of sandy eolian and alluvium deposits that are often coarser Sandy Uplands along the lower footslope and toeslope positions and grade into Loamy Uplands on the valley floor and Bottoms or Sandy Bottoms if there is enough runoff to facilitate water accumulation. Deeper and coarser soils are associated with grassy systems with lesser woody components (see Tables 4 and A.3) that can increase with overgrazing, drought, and/or temperature increases (see   Tables S4 (mmc4.xlsx) and  S5 (mmc5.xlsx)) (Gremer et al., 2018; Munson et al., 2016; Winkler et al., 2019). Sandier areas are also more prone to dune mobilization and associated invasion of exotic herbaceous species tolerant to saltating sands (e.g., Salsola spp.) (see   Tables S4 (mmc4.xlsx) and  S5 (mmc5.xlsx), available online at 10.5066/P92OPRMV) (Draut et al., 2013; Ellwein et al. 2018; Miller et al. 2011; Thomas and Redsteer 2016).

Figure 8.

3-Dimensional model showing predicted soil geomorphic units (SGUs) looking north across the Paradox salt collapse valley in western Colorado. Some notable geologic and landform features closely related to SGUs regionally are noted. A vertical exaggeration of 6x was employed to better highlight topography.

Higher exposed positions dominated by shallower soils and higher rock content (Outcrops, Deep Rocky, Breaks, Shallow, and Very Shallow SGUs) favor more woody species (see Fig. 8, see Tables 4 and A.4) and had less documentation of alternative states, particularly erosional states, than uplands areas with deeper soils and less rock (see  Table S5 (mmc5.xlsx)). Shallow and rocky areas are often rugged (e.g., Fig. 8), with ecological communities adapted to extremes, and usually have less livestock use (Duniway and Palmquist 2020) and overall lower land use impacts due to inaccessibility. However, these soils offer little water storage capacity, and even the drought-tolerant pinyon and juniper tree species dominant in these areas have experienced recent die-offs due to drought and drought-related pests that need better characterization in STMs (Flake and Weisberg 2019; Hartsell et al. 2020; Kannenberg et al. 2021). Additionally, these woodier communities often have documented fire-driven state transitions, particularly in Semiarid Warm and Semiarid Cool ESGs, which need careful future assessment in STMs due to uncertainty in how changing fire regimes are interacting with evolving management practices, pests, and climate change in ways likely to exacerbate future community changes (e.g., Floyd et al. 2004; Romme et al. 2009; Haffey et al. 2018; Hartsell et al. 2020).

In contrast, deeper upland soils that occupy less rugged and more accessible areas have many more documented invasion, encroachment, perennial species loss, and bare-ground related states (see   Tables S4 (mmc4.xlsx) and  S5 (mmc5.xlsx)). Considerable work has examined the degradation of these ESGs by grazing (e.g., Cole et al. 1997; Neff et al. 2005; Belnap et al. 2009) and climate change (e.g., Munson et al. 2011a; Witwicki et al. 2016; Hoover et al. 2021). These studies and our work here also suggest a broad trend toward transitions to states with less biological soil crust, more bare ground gaps, and more shrub-dominated ecosystems that are more susceptible to wind erosion (Okin et al. 2006; Munson et al. 2011a; Webb et al. 2014).

Big sage (Artemisia tridentata), a widespread and dominant shrub species in the western United States (Davies et al. 2011), has seen dramatic losses across its range due to interactions of invasion by annual grasses, fire, and land use (Knick et al. 2003; Schroeder et al. 2004) and is the focus of much research and policy due to sagebrush-obligate grouse species (greater sage-grouse, Centrocercus urophasianus; and Gunnison sage-grouse, Centrocercus minimus; e.g., Schroeder et al. 2004; Arkle et al. 2014; Chambers et al. 2017). We found that big sage subspecies were present in all SGUs but most prevalent in the Semiarid Warm and Semiarid Cool climate zone in Finer Uplands, Clay Uplands, Loamy Uplands, and Saline Uplands in both reference production analysis (see Table A.4) and plot cover data (Table A.5). Interestingly, substantial reference production of big sage subspecies occurred in Sandy Uplands, Sandy Bottoms, Deep Rocky, and Bottoms, but the cover data did not support prevalence of big sagebrush in these SGUs, particularly in Sandy Uplands and Sandy Bottoms. Sage loss in STMs is listed most prominently in Sandy Bottoms and to a lesser degree among seven other SGUs in the semiarid climate zones in regional STMs (  Tables S4 (mmc4.xlsx) and  S5 (mmc5.xlsx)). Although big sage systems have been studied extensively as habitat for sage-grouse management (e.g., Davies et al. 2011; Kachergis et al. 2012; Arkle et al. 2014), broader studies of big sage distribution have focused primarily on climate and ecohydrological relationships (e.g., Schaepfer et al. 2012a, 2012b). Our results largely corroborate the climate niches previously reported for big sagebrush but also help to better establish specific salinity and soil depth constraints that may help to better prioritize management efforts for these ecosystems (Davies et al. 2011; Chambers et al. 2017). Our results do contrast from these restoration frameworks, as our STM analysis indicates that loss of sage is more frequently documented in the Semiarid Cool climate zone (see  Table S5 (mmc5.xlsx)) relative to the Semiarid Warm zone, whereas other studies point toward wetter and cooler areas as priorities for restoration due to better conditions for adapting to ongoing changes in climate and disturbance (Monroe et al. 2020).

Data limitations and uncertainty

ESGs and underlying environmental factors are not able to explain a substantial portion of variability in vegetation and STM data (see Figs. 6 and 7), reflecting some combination of 1) error in input data, new maps, and analysis; 2) potential for omission of other important explanatory data; and 3) some expectation of overlap in characteristics between classes (e.g., we do not expect all grass to be mutually exclusive to one ESG). Clearly, users must consider within-ESG uncertainty in what vegetation and ecological states can be expected for applications using ESGs (see Fig. 6). We emphasize the need to corroborate mapped classes with field confirmation for site-specific management.

The differences in relative dominance among species and functional types between the observed cover (AIM and NRI field data) and the reference conditions documented in the ESD production data are likely due to a combination of methodological differences (fractional cover vs. pounds per acre) and land conditions being departed from the proposed reference state. For example, Saline Bottoms and Bottoms differed notably in trends between cover and reference production, potentially due to their propensity to transition to erosional states with lower water tables that favor loss of productive grass species, possibly reflected in lower relative cover values (e.g., Bailey 1935; McFadden and McAuliffe 1997; Nauman et al. 2018; see   Table S4 (mmc4.xlsx)). Cattle also preferentially use bottoms areas that are more productive and closer to water, especially in saline landscapes (Duniway et al. 2018). Regional packrat midden studies have indicated that livestock grazing-induced state changes have caused widespread regional vegetation change and often homogenize areas (Cole et al. 1997; Fisher et al. 2009). Homogenization of vegetation across soil units potentially contributes to some of the similarity in cover between ESGs observed here. Other work has shown that land use impacts and climate change have also caused regional shifts in vegetation away from presumed reference states (Herrick et al. 2010; Munson et al. 2011b; Gremer et al. 2015; Duniway et al. 2018), which may explain the differences in our observed results. However, despite selected differences, cover data largely corroborated general ecological trends separating ESGs.

Finally, it is notable that the local expert-based characterization of ESDs and soil survey result in data that carry the risk of inconsistent assumptions and in many cases lack direct traceability to primary data or relevant studies. Researchers have challenged how processes documented to cause ecological state changes have been attributed to STMs in ESDs (Twidwell et al. 2013; Bestelmeyer 2015) causing us to focus on directly observable current abiotic properties and ecological states rather than processes attributed in ESDs for analysis. However, understanding ecological processes is the ultimate goal of ecological classification, and the documented ecological states are results of the process(es) driving change (e.g., erosional, woody encroachment, and annual invasion states). Some processes like geomorphic mediation of runoff (e.g., defines bottoms sites) and aeolian transport (e.g., important in sandy groups) are closely tied to what defines the ESGs presented, but others are more subtle. We think that clear definition and explicit mapping of relevant abiotic features across broader scales will help facilitate better links to the relevant studies of causal drivers and processes of ecosystem state change (e.g., restorative management actions, climate change, fire, grazing) that will enable better informed STMs. We have provided examples of many of these process linkages in this discussion, but there is much uncertainty over how ecological systems will behave as climate change and land use create more novel situations (Winkler et al. 2018; Lynch et al. 2021). Clearly, additional primary field data and longitudinal studies will help test and refine ESGs and underlying ESDs and should be pursued in future work.

Applicability to other landscapes and future climates

Although our analysis is specific to the arid and semiarid portions of the Upper Colorado River Basin under recent climate normals, similar soil-geomorphic-climate optimization could be expanded to other areas. The separation of climate from SGUs in our framework can accommodate nonstationary climate trends by allowing for updates of climate thresholds with periodic analysis of new data as conditions change. The SGU concepts built here are likely to mediate ecological processes in other regions and could be assessed using similar analyses. Some SGUs linked to local geologic anomalies may not be present in some regions (e.g., Gypsum), and SGUs may need to be added in some areas, such as petrocalcics in the Basin and Range Province (McAuliffe 1994; Duniway et al. 2007), podzol soils in Appalachia (Nauman et al. 2015a; Nauman et al. 2015b), or Andic soils in the Pacific Northwest (McDaniel et al. 2005), among others.

We suggest that, given the challenges that ESD development has faced, incorporating a more generalized and quantitative level of classification combining local expert knowledge derived from field efforts (Bestelmeyer et al. 2016; Duniway et al. 2016) with quantitative optimization done on a broader regional scale can streamline efforts and provide valuable quantitative support when conflicts arise. Our framework challenges existing soil survey land resource hierarchies with results showing that soils and geomorphology explain variation in productivity and ecological behavior at much broader regional scales than current classification reflects (Salley et al. 2016). We hope that the general framework of SGUs can serve as a more static land unit that can be combined with more dynamic parameters (e.g., climate metrics) to create flexible groupings that can accommodate changing uses and climate trends that will continue to impact ecosystems (Seager et al. 2007; Ault et al. 2016).