Study area

Our study area was centered on the Uintah and Ouray Reservation and surroundings in northeastern Utah (Fig. 1), including three broad geographic areas: the higher elevation Uinta Mountains foothills in the north and the Tavaputs Plateau in the south, and the lower elevation Uinta Basin between these two. The Reservation is divided into 148 management units referred to as range units. The Reservation covers an elevation range of over 2 000 m, from 1 300 m to 3 350 m above sea level. Climate varies greatly, with an arid climate with cold winters and hot summers in the lowest elevations, while higher elevations have a subhumid climate with cold winters and short, cool summers. The lowest elevations in the Uinta Basin receive approximately 150 mm mean annual precipitation (MAP) and 8°C mean annual temperature (MAT) (1 540 m, Ft. Duchesne Station ID:USC00422996, Utah Climate Center 2019). The highest elevations receive approximately 885 mm MAP and 1°C MAT (3 231 m elevation, Brown Duck Station ID:USS0010J30S, Utah Climate Center 2019).

Geologically, the Uinta Basin is dominated by the Duchesne River and Uinta Formations, the Tavaputs Plateau is composed mainly of the Green River Formation, and the Uinta Mountains foothills include Mancos Shale and older Jurassic to Triassic Formations (Hintze et al. 2000). Soils include Aridisols in much of the lower-elevation Uinta Basin, Entisols in middle elevations of the Tavaputs Plateau and Uinta Basin, and Mollisols in higher-elevation, moister sites in the Tavaputs Plateau and Uinta Mountains (Boettinger 2009).

Many distinct vegetation communities occur throughout this region. Generally, the lowest elevations are dominated by saltbush (Atriplex L. spp.) and galleta grass (Pleuraphis jamesii Torr.), and middle elevations include pinyon-juniper (Pinus edulis Engelm. and Juniperus osteosperma Torr.), Wyoming big sagebrush (Artemisia tridentata Nutt. ssp. wyomingensis Beetle & Young), and bunchgrasses (Leymus Hochst. spp., Achnatherum P. Beauv. spp., and others). The highest elevations feature mountain big sagebrush (Artemisia tridentata ssp. vaseyana Rydb.) and other shrubs, lodgepole pine (Pinus contorta Dougl.), aspen (Populus tremuloides Michx.), and a variety of graminoids. Species names and authorities are from the USDA NRCS plants database (USDA–NRCS 2020).

Managed grazing by cattle and other livestock has been limited in recent years (A. Pingree, personal communication, June 21, 2019). However, unmanaged grazing still occurs on the Reservation; a bison herd introduced in 1986 has increased to several hundred individuals (Bates and Hersey 2016), and wild horses, elk, and other ungulates are common in some areas.

Field-sampled vegetation data

We compiled annual plant production data collected at 872 transects across the Uintah and Ouray Reservation, each surveyed once between 2010 and 2017 (see Fig. 1). Transect locations were stratified by soil map unit, with approximately one transect per 1 000 acres of each soil map unit within each range unit. Transect locations were randomly generated before sampling. At each transect, a 30.5-m (100-ft) transect tape was laid in a random orientation, and the current year's growth of grasses and forbs was sampled in ten 0.89-m² hoops located every 3.05 m along the transect. Biomass was sampled through either double-sampling or total harvest of current year's growth, depending on vegetation density (Natural Resources Conservation Service 2003; Herrick et al. 2005). Shrub biomass was sampled by weight unit estimation in two larger circular areas along the transect, each 40.72 m² in size.

Field sampling occurred from approximately May through August of each year. Sampling began at the lowest elevations and moved to progressively higher elevations throughout the season to capture as much of the year's total plant production as possible. The current stage of development for each species at each transect was also noted to estimate total annual production when development was not yet complete. For example, a forb flowering at the time of sampling was assumed to have produced only approximately 75% of its total annual production, and its observed biomass was multiplied by 1.33 to estimate its total annual production (Herrick et al. 2005; S. Green, personal communication, July 16, 2020). When wildlife grazing was evident, we also estimated what proportion of a species' total biomass had been grazed by weighing grazed and ungrazed individuals of similar development and increased observed measurements to estimate the total biomass produced before utilization (Natural Resources Conservation Service 2003; Herrick et al. 2005).

Small samples of all species were collected and oven-dried to calculate dry biomass proportions. Finally, annual ANPP of each transect was calculated by adding the growth measured for each species at the transect. Available forage was calculated by multiplying the total production of each species by a palatability factor for that species, representing the proportion of the species' total production that cattle are expected to use. Palatability factors had previously been established by the Bureau of Indian Affairs to quantify forage utilization by cattle in the region (see  Table S1 (mmc1.docx), available online at …). ANPP and forage values used in this analysis represent oven-dried biomass and are reported in kilograms per hectare. Full vegetation sampling methods are available in supplementary materials.

Calculating edaphic data

We used the NRCS Soil Data Viewer Tool (Natural Resources Conservation Service 2019) to calculate edaphic data in the study area. We summarized these data from the soil surface to a 20-cm depth, calculated as a weighted average of soil components > 5% abundance within each soil map unit. Data were rasterized at 30-m resolution.

Processing Landsat imagery

Landsat imagery available from 1984 to 2018 were accessed from Google Earth Engine (Gorelick et al. 2017). We accessed Thematic Mapper (Landsat 5) from 1984 to 2013, Enhanced Thematic Mapper (Landsat 7) from 1999 to 2018, and Operational Land Imager (Landsat 8) from 2013 through 2018. Landsat data were used to calculate Normalized Difference Vegetation Index (NDVI), formulated as:

where Red and NIR represent spectral reflectance in the red and near-infrared regions, respectively.

Only Tier 1 data, representing the highest-quality images, were selected for analysis. Images were processed to terrain-corrected surface reflectance and then further screened for clouds, snow, and shadows using the CFMask algorithm (Zhu et al. 2012). This process is often used to produce consistent time series such as Landsat Analysis Ready Data (Dwyer et al. 2018). After this imagery selection and processing, we determined the maximum NDVI value attained in each pixel in each yr from 1984 to 2018. NDVI values obtained from the three sensors differ slightly, by < 5% (Teillet et al. 2001; Li et al. 2014; Roy et al. 2016).

Calculating climatic data

Climatic data were obtained from Daymet (Thornton et al. 2018). Precipitation, maximum temperature, and minimum temperature were determined in multiple temporal windows. Vapor pressure deficit was calculated by subtracting actual vapor pressure from saturation vapor pressure (Zotarelli et al. 2010). Originally in 1-km resolution, climatic data were resampled by nearest neighbor to 30-m resolution to align with all other data.

Modeling forage availability and ANPP

We used a random forest model to predict forage availability and ANPP at transects from environmental and climatic covariates. Random forest fits a series of classification trees to random samples of input data (“in-bag observations”) and evaluates accuracy by predicting the values of data not included in a given tree (“out-of-bag observations”) (Breiman 2001; Cutler et al. 2007). Each tree uses only a random subset of predictor variables (“split variables”) per node, making trees more variable and reducing correlation between trees (Cutler et al. 2007). Random forest is appropriate for modeling moderate size datasets with many nonlinear, interacting variables.

We used the randomForest package in R (Liaw and Wiener 2002; R Core Team 2018), building 800 trees with 5 split variables per node. We used only two thirds of the input data for model training and reserved one third of data for model validation. The validation points were therefore never incorporated in model construction and are not to be confused with out-of-bag observations, which vary from tree to tree. Model performance was evaluated by calculating the mean absolute error (MAE) and root mean square error (RMSE) of the model's predictions of validation points not included in modeling. MAE evaluates the average absolute value of errors by the following equation, where represents the predicted value, y_t represents the observed value, and T represents the total number of points:

Table 1

Covariates included as predictors of ANPP and forage availability at transects and the source and processing of those covariates.

RMSE calculates the square of all errors, determines the average squared error, and then evaluates the square root of the average squared error:

Low MAE and RMSE indicate a close fit between predicted and observed values of points. All errors are weighted equally in MAE, whereas RMSE allows large errors to be identified by giving them greater weight. We prioritized the minimization of RMSE over MAE in model selection to reduce the occurrence of large errors.

There is some debate as to how variables should be selected or excluded from random forest models (Behnamian et al. 2017; Fox et al. 2017; Degenhardt et al. 2019). Fox et al. (2017) found, using a dataset similar to ours, that iteratively eliminating variables after attempting to identify unimportant variables did not significantly affect model performance. Therefore, we eliminated only two predictors that random forest variable importance measures strongly suggested were unimportant.

The final variables included in modeling forage and ANPP were all edaphic characteristics (soil available water capacity, cation exchange capacity, depth to restrictive layer, organic matter, pH, sodium adsorption ratio, proportion of sand, silt, and clay; estimated rangeland production), elevation, slope, northness, compound topographic index, estimated tree cover (from Landfire), maximum annual NDVI (from Landsat), and climatic data from Daymet (January–June precipitation sum, May–June mean vapor pressure deficit, May–June mean maximum temperature, and May–June mean minimum temperature). These covariates and their derivations are shown in Table 1.

Bias correction

Random forest predictions often have a systematic bias, where predictions are too high at very low and too low at very high observed values (Zhang and Lu 2012; Xu 2013). Therefore, prediction errors (observed values minus predicted values) tend to be negative at low predicted values and positive at high predicted values. We corrected this bias by employing a method from Xu (2013) and Zhang (2012), which applies a second random forest to the results of the first random forest, modeling the prediction error of the first forest as a function of the prediction. This modeled prediction bias is then subtracted from the first random forest to calculate a bias-corrected prediction.

Modeled forage availability layers

Layers of modeled annual forage availability and ANPP were calculated in R by associating the random forest model with data layers corresponding to all covariates included in the model and calculating pixel-by-pixel predictions annually. This required raster layers corresponding to all covariates included in modeling. These were projected to the same coordinate system (EPSG:4326) and resampled when necessary to 30-m resolution by nearest neighbor resampling.

We calculated layers of modeled forage availability and ANPP annually from 1984 to 2018. To correct for random forest bias, layers of predicted bias were also calculated and then subtracted from the modeled layers to calculate bias-corrected results. These resulting layers were masked where the underlying landscape corresponded to urban areas, agricultural land, or water.

Summarizing forage availability of range units

The Reservation is divided into 148 grazing allotments, referred to as range units (see Fig. 1). To summarize forage availability and ANPP in each unit, we calculated the mean forage availability and ANPP in each range unit in each year by masking the layers to include only pixels falling within a single range unit and calculating the mean and standard deviation of pixel values in each year. This summarized results at a more relevant scale and allowed comparison of forage availability in individual units. We also calculated other relevant attributes of forage availability in range units, including the median forage availability across years, and the maximum and minimum deviations from median forage availability. We calculated the forage coefficient of variation (CV), the ratio of standard deviation to mean, to assess the overall interannual variability of forage in range units. Lastly, we examined relationships between modeled forage availability and modeled ANPP at sampling locations and allotments as a whole.

Trends

We used the time series of mean forage availability in each range unit to test for trends in forage availability over time. For this analysis we used the Mann-Kendall test, a common method to assess nonparametric trends (Mann 1945; Kendall 1948). We also used this test to examine trends in the climatic variables included in modeling.

Random forest model fit and correlates of forage availability

Individual covariates were only weakly correlated with forage availability at transects (Table S2, available online at …), meaning each covariate alone only partially explained the variance in forage availability. The highest correlations were found for annual maximum NDVI (r² = 0.138), NRCS estimated rangeland production of soil map units (r² = 0.116), and May–June vapor pressure deficit (r² = 0.104). Lesser forage availability was associated with typical drought conditions—higher vapor pressure deficit, lower precipitation, and higher temperatures. Random forest variable importance metrics differed slightly from these raw variable correlations and showed tree cover, maximum NDVI, elevation, and precipitation had the greatest explanatory value in predicting forage availability (Fig. 2).

Figure 2.

Random forest model variable importance. Values show the percent increase in mean squared error among random forest trees without a given variable.

Though individual covariates were only weakly correlated with forage, the variance explained by the random forest model was high. After bias correction, r² = 0.72, RMSE = 99.11, and MAE = 64.54 among validation points not included in model training. Fit for the entire dataset including validation and training points was higher (r² = 0.86, RMSE = 73.37, MAE = 44.07 after bias correction) (Fig. 3). Model predictions had a slight bias before bias correction, with prediction errors slightly skewed negative at low modeled values and positive at high modeled values. Correction reduced this systematic bias and improved model performance, reducing RMSE by 27% and MAE by 33% among validation points.

Modeled forage availability

We used the random forest model to calculate layers of modeled forage availability and ANPP annually from 1984 to 2018 at 30-m resolution (Fig. 4a). Forage availability was highly variable spatially, ranging from under 100 kg • ha^–1 to over 800 kg • ha^–1. Generally, available forage was higher at higher elevations in the Uinta Mountains and Tavaputs Plateau than in the Uinta Basin.

Time series of the mean forage availability in two range units with highly contrasting forage dynamics are shown in Figs. 4b and 4c. These units are Steer Ridge, ranging from 2 500 to 2 800 m above sea level in the Tavaputs Plateau region, and Alger Draw, from 1 600 to 1 800 m above sea level in the Uinta Basin region. Vegetation composition varies greatly between these units, with Steer Ridge characterized by mountain big sagebrush, aspen, Gambel oak, and sedges, and Alger Draw dominated by greasewood, saltbushes, snakeweed, halogeton, and basin big sagebrush. These units are outlined in Fig. 4a.

In Steer Ridge, median forage availability was 485 kg • ha^–1 and most years fell within 10% of the median, between 436 and 506 kg • ha^–1 (see Fig. 4b). However, forage declined dramatically in a few years (2002, 2007, and 2014), as low as 314 kg • ha^–1, 35% below the median. In contrast to the drastic forage declines in some years, large increases above the median were not evident. Forage varied differently within Alger Draw, in the Uinta Basin region (see Fig. 4c). Forage availability was centered around a median of 162 kg • ha^–1 in this unit. Most years fell within 17% of the median, but forage increased dramatically in some yr (1995 and 1999), up to 43% above the median. However, large declines below median conditions were not evident in this unit.

Figure 5 shows histograms of these deviations from median conditions colored by geographic region. These histograms demonstrate that the dynamics mentioned earlier are typical for their respective regions, not unique to these range units. More frequently, units in the cooler, moister Uinta Mountains and Tavaputs Plateau had large decreases below median (Fig. 5a) and small increases above median conditions (Fig. 5b). Conversely, units in the warmer, drier Uinta Basin had more frequently small decreases below median but large increases above median. This matches the relationship shown in the Figure 4b and 4c time series results.

The coefficient of variation (CV) of forage, showing the overall interannual variability of available forage, ranged from 6% to 28% among all range units, with a mean of 16%. CV was negatively correlated with median forage availability, indicating that greater median forage availability was associated with lesser interannual variability. This relationship was weak but statistically significant (Fig. 6).

Relationship between forage availability and ANPP

We examined the relationship between forage availability and ANPP using both the plant production data collected at transects and the modeled forage and ANPP results. Forage and ANPP are clearly associated at transects, but with considerable variability (Fig. 7a), indicating that ANPP alone is not an adequate predictor of forage availability. Median forage availability and median ANPP at range units were more closely related (Fig. 7b).

The relationship between modeled forage availability and ANPP varied spatially (Fig. 8). Generally, lower-elevation areas in the Uinta Basin had lower forage to ANPP ratios, whereas higher elevations in the Uinta Mountains and Tavaputs Plateau had higher forage to ANPP ratios. This means available forage constituted a greater proportion of total ANPP in the wetter and cooler high-elevation units, indicative of greater forage quality. Forage to ANPP ratios ranged from 0.23 to 0.63 among all range units.

Trends

We found no range units with statistically significant (P < 0.05) trends in forage availability from 1984–2018, though some approached significant declines. In 21 of the 148 range units, primarily in the southern Tavaputs Plateau, we found declines approaching significance (P values from 0.08–0.25). Figure 9 shows the time series of annual forage availability in one range unit where declines approached statistical significance (P = 0.08). Analyzing climatic trends, we found significant increases in maximum and minimum temperature in approximately one-third of range units (P < 0.05) but found no units with significant changes in precipitation and only three with significant increases in vapor pressure deficit.

Figure 3.

Postcorrection modeled forage versus observed forage among validation points (a); and validation and training points (b). Line of best fit shown in solid blue line, 1:1 relationship shown in dashed line.

Modeled forage availability

We found large deviations between median forage availability and extremes in years with minimum and maximum forage availability (Figs. 4 and 5). On average across range units, minimum forage was 32% below median forage and up to 58% below median in some units. Conversely, maximum forage was on average 33% above median forage and up to 100% above median. Such variability has important management implications—high variability indicates typical conditions have limited significance and points to a greater need for adaptive management responsive to yearly conditions. We also found a negative relationship between median forage and coefficient of variation for forage, meaning variability is most significant in units with lesser forage availability (Fig. 6).

Our results are similar to those from a study that found that across publicly owned western US rangelands, ANPP from 1993 to 2017 varied from 23% below to 28% above mean ANPP (Robinson et al. 2019). However, variability throughout our predominantly shrubland study area is lower than what others have found in grasslands (Swemmer et al. 2007; Wehlage et al. 2016). Whereas the mean CV of forage in Reservation range units was 16%, CV of ANPP in grasslands may typically be closer to 30% (Knapp and Smith 2001) and even 50% in some shortgrass steppe regions (Reeves et al. 2020).

Our results showed that higher-elevation areas in the Uinta Mountains and Tavaputs Plateau experienced large forage availability declines below median conditions in some years but less dramatic increases above median (see Fig. 4b, 5). Conversely, units in the lower-elevation Uinta Basin experienced large spikes in forage availability above the median in some years but lesser declines below median conditions (see Fig. 4c, 5). Grazing management should be tailored to these local forage availability dynamics, anticipating the sharp forage availability declines in the Uinta Mountains and Tavaputs Plateau units during drought and the spikes during favorable conditions in the Uinta Basin. Unsurprisingly, increased forage availability was correlated with wetter, cooler conditions (see Table S2).

The relative dominance of perennial and annual vegetation in these regions could contribute to the differing forage availability dynamics (Schmidt and Karnieli 2000). For example, annual exotic cheatgrass (Bromus tectorum L.) is most prominent in the Uinta Basin and may capitalize on favorable conditions and produce additional biomass more readily than perennial species (Hardegree et al. 2018). In addition, many dominant plants in the Uinta Basin perform C₄ photosynthesis, whereas these species are nearly absent in the Uinta Mountains and Tavaputs Plateau. C₄ species typically tolerate drought better than C₃ species (Ward et al. 1999) and may grow more rapidly under more favorable conditions (Nippert et al. 2007), closely matching the differing forage availability dynamics we observed between these regions.

Figure 4.

Example surface of modeled forage availability in 1984 (a) and time series results of mean forage availability in the Steer Ridge (b) and Alger Draw (c) range units. Borders of the Reservation range units are shown in (a), with Steer Ridge and Alger Draw outlined in red in the south and center, respectively.

Figure 5.

Histogram of maximum forage decrease below median forage in range units (a) and maximum forage increase above median forage (b), colored by region.

Trends

We found no statistically significant (P < 0.05) trends in forage availability from 1984 to 2018, though we found trends approaching significance in some range units (Fig. 9). Others have also not found clear NPP trends in recent history across western US rangelands (Robinson et al. 2019). Shrublands, particularly in the US southwest, may have experienced particularly little NPP change on the time scales studied (Hicke et al. 2002). Though climate change is predicted to significantly impact rangeland productivity (Reeves et al. 2014), these changes may not yet be evident in the time frame we analyzed. We found significant temperature increases in our area but no precipitation trends, which would more drastically impact productivity.

Figure 6.

Coefficient of variation (CV) of forage availability in range units versus median forage availability of range units. Line of best fit shown in blue line.

Livestock and wildlife can also influence productivity (Menke & Bradford 1992; Fleischner 1994), but managed livestock grazing in our study area has been limited (A. Pingree, personal communication, June 21, 2019). Bison populations have grown in the Tavaputs Plateau (Bates and Hersey 2016), where many range units had weak forage availability declines (P < 0.25), but these factors are difficult to link. Regardless, trends in forage availability can continue to be monitored using the model, particularly as climate, grazing, or wildlife management may change. Though any model may have limited accuracy under novel conditions, random forest can successfully be extrapolated beyond training conditions (Carrasco et al. 2015). However, with continued climate and vegetation changes it is likely that the model will be improved by recalibration in the future based on updated field sampling. Nonetheless, this is still vastly more feasible than annually measuring forage over vast and varied landscapes.

Figure 7.

Forage availability versus annual net primary productivity (ANPP) as measured at individual transects (a) and median forage availability versus median ANPP as calculated from modeled results in each range unit (b). Line of best fit shown in blue.

Random forest model fit and correlates of forage availability

Individual covariates were only weakly associated with forage availability measured at transects, with the strongest association found for annual maximum NDVI (r² = 0.138). Unsurprisingly, lower forage availability was associated with higher vapor pressure deficit, lower precipitation, and higher temperatures, which are typical drought conditions that reduce productivity (Roby et al. 2020). Fit achieved by the random forest model was generally high, achieving r² = 0.72 among validation points and r² = 0.86 for training and validation points after bias correction (see Fig. 3). Since bias correction improved model fit, we recommend others follow random forest bias correction methods (Zhang and Lu 2012; Xu 2013) when appropriate.

Predicting forage availability from annual maximum NDVI is clearly not appropriate across our study area. Others have found stronger correlations between NPP and NDVI (e.g., Schloss et al. 1999; Paruelo et al. 2000), but NDVI is likely a weaker predictor of forage availability, especially across complex habitats covering multiple plant communities and soil types (Garroutte et al. 2016). Our results show that NDVI in conjunction with edaphic, topographic, and climatic covariates, which also impact productivity (Wessels et al. 2007; Popp et al. 2009; Mitchell 2010), can predict forage with acceptable accuracy in a model such as random forest.

Relationship between forage availability and ANPP

Across the Reservation, there was considerable variation in the relationship between forage availability and ANPP at the level of transects (Fig. 7a). Since ANPP does not fully explain forage availability, directly modeling forage with multiple covariates is preferable to considering ANPP as a proxy, as is sometimes done (Paruelo et al. 2000; Mitchell 2010). We were able to directly model forage availability because our field data identified vegetation to the species level, and the palatability of these species to cattle was previously determined for the region (see  Table S1 (mmc1.docx)).

Forage and ANPP differences were less pronounced when modeled ANPP and forage availability results were aggregated across range units (Fig. 7b). However, the ratio between the two did vary in a somewhat predictable pattern, with higher-elevation units in the Uinta Mountains and Tavaputs Plateau having a higher forage to ANPP ratio (Fig. 8). These forage-to-ANPP ratios were in line with previous literature, which reported a ratio of 0.37–0.54 in shrublands and higher in sites with more grass and less shrub cover (Mueggler & Stewart 1981; Hirata et al. 2005). In our study area, lower elevations are dominated by saltbush and sagebrush shrublands, while higher elevations include more meadow-like sites with less shrub cover and therefore higher forage to ANPP ratios.

Figure 8.

Forage availability–to–annual net primary productivity ratio of range units, from median modeled results in each range unit.