A Quantitative Soil-Geomorphic Framework for Developing and Mapping Ecological Site Groups

ABSTRACT Land management decisions need context about how landscapes will respond to different circumstances or actions. As ecologists' understanding of nonlinear ecological dynamics has evolved into state-and-transition models (STMs), they have put more emphasis on defining and mapping the soil, geomorphological, and climate parameters that mediate these dynamics. The US Department of Agriculture Natural Resources Conservation Service ecological site descriptions (ESDs) have become the foremost system in classifying lands into ecological units based on STMs. However, an exhaustive inventory of ESDs has proved challenging to complete in the United States, and there have been questions about the consistency of detail in areas completed and the ability to objectively support some assertions made in existing ESDs. To address these issues, this study examines ESDs in the diverse Upper Colorado River region, where ESDs are only partially complete, to look at quantitative approaches to generalizing ecological site concepts based on unifying underlying soil, geomorphology, and climate patterns. Using existing ESDs and vegetation monitoring plot data, results show that a simple hierarchical soil geomorphic unit (SGU) framework based on topographic mediation of moisture, soil salinity, soil depth, slope, rock content, and soil texture can represent much of the ecological dynamics cataloged in ESDs. Analyses of reference plant production data, ecological state attribution, and regional monitoring data show that the new SGUs represent more variation than common climate parameters. This study also included predictively mapping SGUs at 30-m resolution (Kappa of 0.53, 74% agreement with top two predictions in validation). An optimized combination of SGUs with climate zones derived from an aridity index and maximum temperature of the hottest month resulted in an ecological site group framework that condensed over 826 unique ecological site records at various stages of completeness in the regional soil survey down to 35 intuitive and mappable ecological site groups.


Introduction
Understanding how both ecological potential and behavior of lands vary spatially is critical for land management and policy decisions. In recent decades, frameworks for grouping rangelands into ecologically and managerially meaningful units have focused on classifying areas by abiotic features that control ecological potential and behavior under the paradigm of state-transition dynamics driven by nonlinear thresholds ( Winthers et al. 2005 ;Bestelmeyer 20 06 ;Herrick et al. 20 06 ;Bestelmeyer et al. 2009 ;Twidwell et  2013 ). Similar developments in fire science have recognized the importance of natural potential but place less emphasis on soil and geomorphic drivers ( Rollins 2009 ). Of the major land potential frameworks, the ecological site description (ESD) inventory, managed by the US Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS), has been endorsed broadly by government agencies for supporting land management policy ( Caudle et al. 2013 ). Central to ESDs is the development of reference conditions that reflect ecological potential and are indicative of desirable site conditions. The concept of a reference has evolved from an estimate of precolonial historic reference vegetation community to more realistic sets of functional site parameter ranges that maintain ecological integrity while considering pragmatic goals appropriate to various management objectives ( Pyke et al. 2002 ;Brown and Herrick 2016 ;Herrick et al. 2019  erence concepts allow for broad interpretations of land condition with monitoring programs (e.g., Miller 2008 ;Herrick et al. 2010 ).
ESDs also include a state-and-transition model (STM) that catalogs potential alternative ecological states and transition pathways (or drivers of change) and provides managers a roadmap for decision making in a given landscape setting Bestelmeyer et al. 2009 ). Despite 20 + years of active development, there are challenges to producing ESDs across the United States. A major challenge is the tremendous number of individual ESDs needed for full US coverage, which is driven by a combination of the relatively fine conceptual and spatial scale of most ESDs. Production of ESDs occurs separately within each NRCS major land resource area (MLRA; n > 20 0; USDA-NRCS 20 06 ). The large number of ESDs and quickly evolving ecological science make incorporation of supporting evidence for each individual ESD and STM difficult-a perpetual work in progress ( Bestelmeyer 2015 ). Characterization of transitions between states often lacks available longitudinal studies needed to fully understand drivers and subsequent threshold behavior in the changing earth system Bestelmeyer and Briske 2012 ;Twidwell et al. 2013 ). A historic focus on livestock has led to more attention on grazing relative to fire, climate, and other change drivers in published STM transitions ( Twidwell et al. 2013 ). In addition to MLRA, state boundaries also have dictated administration of ESDs in some cases. The political and geographic divides that dictate ESD development can create added cost and administrative burden for both producers and users. Several groups have concluded that more generalized ecological management units with less redundancy in STMs result in useful groupings that can be more easily mapped and linked with relevant studies for STM development Bestelmeyer et al. 2016 ;Duniway et al. 2016 ;Salley et al., 2016 ;Stringham et al. 2016 ;Maynard et al. 2019 ).
Another concern with the current ESD inventory is spatial scale. Soil survey map units that provide the current spatial representation for ESDs often include multiple soil components linked to different ESDs in one polygon that can be quite large (from hundreds up to tens of thousands of hectares). Many ecological processes characterized in STMs happen at a different range of scales (i.e., submeter to hundreds of meters, Bestelmeyer et al. 2011 ), which often creates a mismatch in scale with ecological sites as mapped in traditional soil survey. Several recent efforts have attempted to bridge this scale gap using digital soil mapping approaches ( Maynard and Karl 2015 ;Nauman and Duniway 2016 ;Maynard and Karl 2017 ;Maynard et al. 2019 ) that can also be integrated better with new remote sensing vegetation datasets (e.g., Zhou et al. 2020 ). These effort s have shown potential for improving ecological site mapping but are still limited by the scope and consistency limitations of the ESDs themselves. This highlights a need for a more generalized soil-geomorphic-climate template (e.g., building on Duniway et al. 2010 ;Nauman and Duniway 2016 ) that can be consistently mapped at finer resolution across geopolitical boundaries.
Our objective in this paper is to produce a quantitative soil-geomorphic framework mappable across broad extents with landscape-scale detail as the foundation for a generalized ecological site framework. To achieve this, we employed a variety of data and information available for all ESDs in our study area using multivariate statistical techniques to test and optimize decisions made to define new ecological site groups (ESGs) based on climate, soil, and geomorphic properties ( Fig. 1 ). We undertook this for the diverse Upper Colorado River region, where resource degradation is well documented at a broad scale ( Miller 2008 ;Herrick et al. 2010 ;Miller et al. 2011 ;Munson et al. 2011b ), but ESD development is not fully complete. We used reference state vegetation data and documented ecological states from published STMs as response variables for optimizing group classification. We further corroborated resulting ESGs using field-based monitoring data. To facilitate adaptability, we separated climate classification from soil-geomorphic mapping (suggested in Salley et al. 2016 ), allowing for easier updates to both as conditions in the future and knowledge improve. The outcome is a hierarchical soil geomorphic unit (SGU) framework for the region that can be mapped with available soils and geographic data that provide a robust starting point for similar effort s in other regions.

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. 1 a). 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. 1 a) 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 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.

Short name Species included as identified in ESDs
Production/total cover Total site production of all species/total foliar cover of all species Grass All grasses Shrub All shrubs (includes subshrubs) Forb All a custom look-up table to re-correlate the species-level data to appropriate functional groups ( Table S1 , 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 , 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. We tabulated possible ecological states for each ESD from published STMs in EDIT (see Fig. 1 a). 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 , 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( Anderson , 2017Anderson 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 boxplots 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 Duniway et al. 2016 ) and initially classified all ESDs in our study to match these groups (see Fig. 1 b). 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. 1 c). 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. 1 ( 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 ∼30m 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. 1 e). 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 −20 0 0 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. 1 f). 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 NRCSbased 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 SGUclimate 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 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 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.
of possible combinations to under 1 x 10 8 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.
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 ) 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. 1 f and 1 g).

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. 1 h). 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. 1 i). 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 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. 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.
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  ) 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. The SGU map shows large areas of shallow soils (pinker areas, see Fig. 4 a) 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. 4 a). 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. 4 b). All maps are available online ( Nauman 2021 at https://doi.org/10.50 6 6/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 Bot- 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 toms, 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.

Models with SGUs aggregated simultaneously across climate zones.
Trials combining climate breaks and SGUs aggregations without nesting included ∼2.7 million possible models. Analysis of the AICc distribution of models showed a more concentrated pool of bet-ter candidate models than the nested schemas with relative likelihood dropping below 0.5 for the third best model. On the basis of this higher concentration of stronger model candidates in the AICc weights, we retained models with a relative likelihood of > 0.1 for STM state-based trials ( n = 48).

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 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). 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.

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 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 .

ESD indicates ecological site description.
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 pat-terns 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.
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 ), but ESGs in the Semiarid Cool climate zone had less overall relative occurrence of invaded states ( Table S5 ). 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 snakeweed 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 pro-duction showed more congruent trends to tree cover among ESGs (see Fig. 6 c and 6 d). 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.

Relative explanatory power of SGUs, climate, and interactions
NPMANOVA models explaining reference production data, ecological state occurrence, and plot cover data showed that soil geomorphic units generally explain more variation than climate parameters ( Fig. 7 ). Models explained 11 −56% of variation, but we expected that some vegetation indicators and ecological states will be shared among subsets of the classes, thus limiting how much variation could be explained. Models that include both soil geomorphic units and climate parameters tend to be stronger than models with just SGUs or just the climate parameters (see Fig. 7 a −7c). Models with SGU by climate interaction terms explained the most variation (i.e., 2x_sgu_clim models, see Fig. 7 ). The ESG models, meant to capture the synergy of using both SGU and climate data in one grouping, show similar strengths to models with the SGU and climate parameters but fail to encompass all the variation explained by models with SGU-climate interaction terms (see Fig. 7 a −7c). When testing NPMANOVA models for the marginal influence of SGUs and climate parameters (sgu_clim models, see Fig. 7 d −7e), SGUs explain more variance than any climate variable (see Fig. 7 d-e; P < 0.001 for all SGU estimates). The most influential climate variable in all sgu_clim models only represented 2.4% of marginal variation (aridity index in cover model, P < 0.001), and no clear best climate variable emerged in models.

Discussion
We present a quantitative approach for developing a landclassification framework based on ecosystem potential and behavior for a 349,500-km 2 dryland area of the southwestern United States, building on qualitative a priori frameworks done previously (e.g., McAuliffe 1994 ;Bestelmeyer et al. 2011 ;Michaud et al. 2013 ;Duniway et al. 2016 ;Stringham et al. 2016 ). Our approach differs from previous efforts in the scale of the effort, systematic quantit ative supporting analysis and detailed 30 m mapping. Classes break up the landscape by productivity (ESD reference community production) and ecosystem behavior documented in observable ecological states that occur as a result of drivers of change (e.g., erosion, grazing, fire, land use disturbance, climate variation, management actions).
Results demonstrate that basic soil-geomorphological groupings are generally more important than climate variation in characterizing the ecosystems assessed, but the best classifications carefully consider interactions among all abiotic factors. Although we selected one final model for practical purposes, we stress that from an information theoretics standpoint, there is no single right model ( Anderson and Burnham 2002 ;Burnham and Anderson 2004 ), especially for a land-classification system designed for pragmatic land management that can include widely varying objectives. Considerable uncertainties exist in the expert-synthesized ESD data (no primary plot data available to evaluate) and soil survey data used in this analysis, warranting future testing with more quantitative primary data. We built the class definitions as transparently as possible to serve as a springboard for future refinements. There is a large need for robust land potential-based classification systems, as well as accurate maps of those classifications Duniway et al. 2016 ;Stringham et al. 2016 ;Maynard et al. 2019 ). Critical to both these needs is workflows (e.g., Fig. 1 ) that can merge existing data and expert knowledge into frameworks supported by quantitative analysis and mapping at scales consistent with management. Figure 6. Boxplots showing distribution of reference state mean ecological site description production ( a, c ) and line-point-intercept plot cover ( b, d ) for ecological site groups in the Semiarid Warm climate zone. Boxplots include a, total production, b, total foliar cover, c, tree production, and d, tree cover.
These four ESDs are in similar soils (Sand-Sandy Loams, usually a calcic, moderately deep to very deep, nonsaline); climate (generally 8 −11" precipitation with some slight difference in the low or high value and similar seasonality); and landform (generally structural benches, alluvial fans, mesas, alluvial flats, ∼1 −12% slopes), and they are mapped in the same geographic area (southeast Utah). The work here groups these similar ESDs into the Arid Warm Loamy Uplands/Sandy Uplands ESG, supporting a premise that Figure 7. a −c, Comparing nonparametric multivariate analysis of variance models for total explained fractional variation for differing sets of explanatory variables for the reference community production data (dark blue), state-and-transition model states (gray), and plot cover data (yellow). d −f, Plots of the marginal fractional variation explained by soil geomorphic units and climate variables in the sgu_clim models for each dataset. Terms: 'clim' = climate, 'sgu' = soil geomorphic units, 'esg' = ecological site groups, '2x_sgu_clim' = model with sgu and two-way interactions of sgus with climate variables, 'ai' = aridity index, 'pptrt' = monsoon/total precipitation ratio, 'maxt' = mean maximum temperature of hottest month, 'mint' = mean minimum temperature of the coldest month.
there are broader potential possibilities for reference state communities in these sites. This highlights the need for a more systematic, quantitative regional analysis as undertaken in this project to keep ESDs from proliferating into a myriad of overlapping concepts that are confusing for even the most adept users. The apparent tendency to create narrow definitions of ESDs that do not fully reflect the potential diversity of reference communities is perhaps a lingering bias of the rangeland community described well by Fuhlendorf et al. (2012) .

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 for-mation 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 ), 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 .
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  and S5) ( 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 and S5, available online at 10.50 6 6/P92OPRMV ) ( Draut et al., 2013 ;Ellwein et al. 2018 ;Miller et al. 2011 ;Thomas and Redsteer 2016 ).
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 ). 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 and S5). Considerable work has examined the degradation of these ESGs by grazing (e.g., Cole et al. 1997 ;Neff et al. 20 05 ;Belnap et al. 20 09 ) 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 and S5 ). 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. 2012aSchaepfer et al. , 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 effort s f or 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 ) 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). Cattle also preferentially use bottoms areas that are more productive and closer to water, especially in saline landscapes . 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 effort s 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 ).

Implications
As patterns of land use and climate change evolve quickly across broad scales, direct management efforts need frameworks that are consistent, systematic, and mechanistic ( Copeland et al. 2017 ). We systematized a wide array of rangeland and soil knowledge and data to quantitatively classify ecological site groups (ESGs) that represent systems of varied ecological behavior across regional extents. Additionally, 30-m resolution mapping of ESGs represents variation at a landscape scale that allows for more efficient field planning and appropriate ecological interpretation of the quickly expanding suite of high-resolution remotely sensed vegetation monitoring datasets (e.g., Zhou et al. 2020 ). This fills a shortcoming of conventional soil surveys that often generalize multiple contrasting soils within map unit polygons at scales too coarse for many remotely sensed data products and field operations.
Analysis revealed substantial redundancy and ambiguity in current ecological sites, highlighting the need for more systematic analysis of relevant ecological indicators to guide ecological site development. By separating more static soil and geomorphic properties from nonstationary climate parameters and building a quantitative workflow for optimizing class thresholds, the framework presented provides concise and transparent classes that can be easily updated to assimilate new climate conditions and expanding knowledge and data. The resulting ESGs also allow for assimilation of new data in combination with existing STMs to document ecological dynamics in generalized STMs. New ESGs have immediate utility for monitoring sample design, remote sensing analysis, and land management planning, particularly where ecological sites do not already exist. To facilitate use of the ESGs developed here, the ESG and SGU maps are available online ( Nauman 2021 ) and we provide a table linking ESGs to existing ESDs ( Table S6 , available online at 10.50 6 6/P92OPRMV ). Together, the maps, links to exiting ESDs, and ESG data provided here (see Tables  4, A .2, A .3, S5) can be used by land and resource managers to better understand land potential and ecological dynamics in the UCRB.

Declaration of Competing Interest
None.

Acknowledgments
We acknowledge the National Cooperative Soil Survey program and managing agency (USDA-NRCS) for the incredible resources they have made available in the soil survey program. Their databases are truly an invaluable resource. We also thank Dr. Skye Wills (USDA-NRCS), Dr. Henry Ferguson (formerly USDA-NRCS), and Dr. Jonathan Maynard (USDA-ARS) for assistance in acquiring and preparing NASIS field soil observations. Jeb Williamson (New Mexico State University, Jornada Experimental Range) kindly provided queries of ecological site description tables. We thank Dr. Shawn Salley (USDA-NRCS), journal editors, and two anonymous reviewers for their valuable feedback. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US government. Mapping data for this study are available online ( Nauman, 2021 : https://doi.org/10. 50 6 6/P92OPRMV ). Code developed for this study is also available at https://github.com/usgs/Predictive-Soil-Mapping/tree/master /ESGs . Very low productivity areas with vegetation in sparse areas and spread out in pockets or fissures. Outcrops with more fractured bedrock can support more vegetation.

Riparian
Variety of soils in floodplains or areas with perennial plant or tree available water tables or surface water.
Aridification, gullying, or channelization can cause these sites to irreversibly revert to uplands.

Saline Bottoms
Gently sloping, low-lying areas that receive excess moisture beyond ambient precipitation (run-on or subsurface moisture). Most have ephemeral washes and streams (not perennial). Soils are influenced by salts and have subsurface * soil electrical conductivity (ec) > 4 dS/m (sat. paste).
Higher productivity than other saline groups. Alkali sacaton and black greasewood are common dominant species. Generally more grass-dominated when in a reference state.
Gullying or channelization can lead to alternative states or cause these sites to irreversibly revert to uplands. Salts also make them less resilient to surface disturbance. Greasewood and other shrubs often increase in alternative states.

Sandy Bottoms
Other gently sloping, low-lying areas that receive excess moisture beyond ambient precipitation (run-on or subsurface moisture).
Most have ephemeral washes and streams (not perennial). Soils are sandier, averaging > 50% sand and < 27% clay in both surface * * and subsurface horizons.
Diverse shrubs and C4 grasses often dominate. These sites have higher productivities than upland counterparts. Can support basin big sage in semiarid climate zones (aridity index ∼> 0.144).
Often more prone to bare ground exposure and associated wind erosion. Can become an aeolian sand source for downwind dunes. Also highly prone to loss of perennial species.

Bottoms
Other gently sloping, low-lying areas that receive excess moisture beyond ambient precipitation (run-on or subsurface moisture).
Most have ephemeral washes and streams (not perennial).
Dominated by grasses and shrubs associated with run-in landscape settings (higher surface or ground water available). Basin big sage can dominate. These sites have higher productivities.
Gullying or channelization can cause state transition or even irreversible reversion to uplands. Woody encroachment is also commonly observed in these areas.

Gypsum
Upland * * * areas with soils averaging > 5% gypsum in surface or > 10% gypsum in subsurface, but with surface sodium adsorption ratio < 8. These areas are often hilly badlands but can be more gentle terrain as well.
Subshrublands with limited grasses dominated by C4 species and low overall productivities. Species composition determined by gypsum tolerance. Often have very high biological soil crust cover.
Favor biological soil crust development. Have the least number of documented alternative states, indicating a high resistance to state change. Limited annual and shrub invasions have been observed.

Saline Hills
Other upland areas that are highly salt limited (often sodic), erosion features common, often dissected badland hillslopes. These soils include surface sodium adsorption ratios > 7 and/or average ec > 4 dS/m in surface or average ec > 8 dS/m in subsurface.
Mat, Castle Valley, and Gardner's saltbush often dominate with associated salt-tolerant species. Low productivity with even grass and shrub production in reference communities.

Saline Uplands
Other uplands with moderate salt limitations including average surface ec > 1.5 dS/m or average subsurface ec > 2.
Less susceptible to herbaceous and woody invasion, as well as erosion and bare ground, than similar SGUs with more or less salinity.

Breaks
Other uplands on steep slopes ( > 35%) and rocky soils with > 40% (by vol.) rock content in surface soil horizons.
Very low productivity areas that favor woody species or resilient forbs. Vegetation is often sparse and limited by unstable slopes, poor water retention, and high rock content.
Particularly susceptible to cheatgrass and other annual invasions. Few other alternative states issues observed.

Very Shallow
Other uplands with soils < 30 cm depth until a bedrock contact. Sites are often rocky and rugged.
Generally low production with even mix of trees (above a certain aridity level), shrubs, and grasses. Blackbrush can dominate in drier areas.
Drought prone and susceptible to annual invasion. Can have bare ground states, erosion issues, and perennial loss of both grass and woody species.

Shallow
Other uplands with soils < 55 cm to a bedrock contact.
Commonly low to moderately productive pinyon-juniper woodlands but also support substantial grass and shrub components that vary in relative abundance by climate. Blackbrush can dominate in drier areas.
Drought prone and susceptible to annual invasion. More woody encroachment than Very Shallow. Bare ground, biocrust loss, eroded, and perennial loss states possible.

Deep Rocky
Other uplands with soils that average > 30% rock fragments by volume in either surface or subsurface horizons. Often rugged topography but can be high-energy alluvial deposits. These soils also tend to have high calcium carbonate contents.
Exhibit a wide variety of dominant grasses, shrubs and trees, including blackbrush, big sagebrush, and juniper at lower elevations. Generally moderate to moderately high production. Species composition is generally very mixed among species and functional groups.
High propensity for herbaceous invasion, moderate for woody encroachment. Resistant to erosional states, but moderately susceptible to bare ground and perennial loss states
Productive savannahs and grasslands often dominated by grasses more adapted to shrink-swell soils. However, big sage can dominate these areas in wetter climates.
Moderate water erosion, herbaceous invasion, and bare ground state risk. Common loss of perennial grasses and woody encroachment.

Sandy Uplands
Other uplands with very sandy eolian and alluvial deposits that average > 75% sand in both surface and subsurface horizons. These soils are generally quite young and low in carbonates ( < 10%, but usually < 5).
Productive savannas and grasslands with substantial shrub component (primarily four-wing saltbush, but with some sand sage, blackbrush, Ephedra spp., and big sage on wetter sites). Blackbrush can dominate this group as a long-term "state" but is less common than on shallow sites or sites with calcics and slightly finer textures. Dunes have been described as a reference state possibility for driest and most exposed areas. Big sage can also dominate in wetter climates.
Drought and disturbance can cause severe wind erosion and dune mobilization. Sites with water erosion issues have also been observed. High propensity for annual invasion, woody encroachment, and perennial species loss (particularly grass). Also moderate risk of bare ground states.

Loamy Uplands
Other uplands with surface soil textures of Sand, Loamy Sand, or Sandy Loam, but finer subsoil field textures or carbonate content higher than 10%.
Grasslands and savannah communities. Some areas have blackbrush communities that can dominate, but often in mosaic with grasslands as a long-term state. Big sage can also dominate.
Similar to Sandy Uplands, but with less risk of most alternative states and no erosional states related to water erosion documented.

Finer Uplands
Other uplands that tend to have finer loamy textures.
Savannas and shrublands with grasses, these are mostly dominated by Wyoming big sage at middle elevations, but include some sites dominated by winterfat.
High risk of herbaceous invasion, woody encroachment, perennial species loss, and bare ground states. Some documentation of eroded states.