Rangelands provide numerous ecosystem services, including forage for livestock grazing. Effective grazing management requires measuring forage availability, which influences the level of grazing that can be sustained while maintaining healthy ecological conditions. However, spatiotemporal variability makes such determinations of forage quantity difficult, potentially hindering optimal management. These determinations are especially difficult across large, remote areas. To address this, we developed an approach using data from a one-time sampling of vegetation throughout the Uintah and Ouray Reservation in northeast Utah. By associating these data in a random forest model with environmental and climatic covariates that vary annually, we produced yearly predictions of forage availability on a pixel-by-pixel basis for the Reservation and surroundings from 1984 to 2018. This method addresses and quantifies the spatiotemporal variability of available forage. The model confirms that forage availability is highly variable throughout the area. On average, forage availability in Reservation management units declined as much as 32% below median availability in some years and increased 33% above median availability in others. Moreover, some management units experienced large increases in favorable years but less significant declines in unfavorable years, while the opposite was true in others. In comparison to determining a single “typical” forage availability of management units, recognizing this inherent variability and how it differs between units provides a fuller picture of the range of possible forage availability. This can improve grazing management by revealing how forage quantities vary from year to year and may help avoid forage overutilization during unfavorable years such as drought. The model can continue to be used into the future to monitor vegetation trends, though with ongoing climate and vegetation changes periodic recalibration may be necessary. In addition, the modeling method may be applicable to other similar study systems.
Rangelands, including grasslands, savannas, shrublands, and steppe, cover about 650 million acres in the western United States (Reeves and Mitchell 2011) and provide numerous ecological and social benefits (Reid et al. 2008). Though general methods to assess key indicators of rangeland health such as vegetation composition, bare ground cover, and annual production are well established (Galt et al. 2000; Herrick et al. 2005; Pellant et al. 2005), accurate and economical inventory and monitoring of rangelands still pose basic challenges (Booth and Tueller 2003).
Annual production and subsequent forage availability are among the most widely studied rangeland attributes, given their inherent link to stocking rates (Holechek 1988). Vegetation production is also an indicator of ecological condition, as declining productivity can signal land degradation, potentially caused by grazing (Bai et al. 2008; Wessels et al. 2008). Appropriate stocking rates help prevent overutilization and the biotic and abiotic degradation that can follow (Fuls 1992; Menke and Bradford 1992).
Accurately measuring forage availability within management units and assigning appropriate stocking rates is essential to grazing management. However, spatiotemporal variability in vegetation production throughout units makes determining forage availability difficult across large areas. Given this variability, any single stocking rate is a simplification (Stoddart 1960), but a “typical” stocking rate or grazing capacity is a useful starting point for management (Holechek 1988; Holechek and Pieper 1992; Galt et al. 2000). Spatial variability in vegetation production can be partially addressed by stratifying field sampling locations in a management unit by soil type or ecological sites (Karl and Herrick 2010). However, vegetation can vary within these units due to microclimatic and topographic variation. Temporal variability is difficult to account for without multiple years of data collection, and field vegetation sampling may inevitably occur in unusually wet or dry years. Therefore, a presumed typical production must be inferred from collected data (Holechek 1988), which can lead to spurious assumptions about typical conditions or trends through time.
Remote sensing is an attractive tool to address rangeland vegetation sampling challenges given its widespread availability through space and time, allowing this spatiotemporal variability to be measured. Remote sensing has been used for decades to determine indicators such as ground cover (Booth and Tueller 2003; Boswell et al. 2017; Ford et al. 2019), land degradation (Allbed and Kumar 2013), and total vegetation production (Hunt Jr. et al. 2003; Running et al. 2004; Del Grosso et al. 2008). However, a gap still exists in using remote sensing to determine forage availability in rangelands.
Most commonly, spectral indices like normalized difference vegetation index (NDVI) are assumed to represent available forage, given its association with total aboveground annual net primary productivity (ANPP) (Paruelo et al. 2000; Mitchell 2010; Borowik et al. 2013). However, the proportion of vegetation that is palatable and available to livestock as forage varies widely between vegetation communities (Miller and Krueger 1976; Mueggler and Stewart 1981), so ANPP may not be an appropriate proxy for forage availability across environments dominated by multiple vegetation types or complex vegetation mixtures. Therefore, directly determining forage availability is preferable to treating ANPP as a proxy. Furthermore, NDVI alone is an imperfect predictor of vegetation biomass, and its accuracy varies by soil, topography, and habitat type (Garroutte et al. 2016).
Here, we detail a method to predict annual forage availability and ANPP throughout more than 2 million ha (5 million acres) of the Uintah and Ouray Reservation and surroundings in northeastern Utah by fitting climatic, topographic, and edaphic covariates to plant production data collected in the field. This location is an ideal study area for this analysis given its wide elevation range, diversity of dominant vegetation types, high interannual variability in vegetation production, and the availability of a large plant production dataset collected in the field.
We aimed to enhance a typical rangeland inventory by producing a time series of high-resolution gridded layers of forage availability and ANPP. This was intended to address the following goals:
Assess the range of variability in forage availability in the study area, and analyze differences in variability between management units.
Examine whether trends in forage availability are apparent throughout the study period.
Interpret the significance of forage variability and trends for grazing management.
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-m2 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 m2 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, yt represents the observed value, and T represents the total number of points:
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.
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.
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 (r2 = 0.138), NRCS estimated rangeland production of soil map units (r2 = 0.116), and May–June vapor pressure deficit (r2 = 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).
Though individual covariates were only weakly correlated with forage, the variance explained by the random forest model was high. After bias correction, r2 = 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 (r2 = 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.
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.
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 C4 photosynthesis, whereas these species are nearly absent in the Uinta Mountains and Tavaputs Plateau. C4 species typically tolerate drought better than C3 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.
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.
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.
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 (r2 = 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 r2 = 0.72 among validation points and r2 = 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.
Average forage availability has limited value in grazing management—not planning for declines in drought years could promote overgrazing (Menke and Bradford 1992), as potentially less than half of median forage is available in these years. On the other hand, increasing stocking will allow the utilization of the additional available forage in favorable years. We found differing forage availability dynamics between the Uinta Basin region and the Uinta Mountains and Tavaputs Plateau regions, with the former more likely to experience forage availability spikes and the latter more susceptible to large declines. This spatial variation suggests management should plan for these local conditions, significantly reducing stocking during drought in the Uinta Mountains and Tavaputs Plateau and increasing stocking in favorable years in the Uinta Basin based on this known variability. Since large forage spikes in the Uinta Basin may represent an increase in cheatgrass productivity, increasing stocking in this area as a targeted-grazing method could also fulfill management goals by controlling cheatgrass populations and wildfire risk (Diamond et al. 2009). Advancements in long lead forecasting (e.g., National Weather Service Climate Prediction Center; https://www.cpc.ncep.noaa.gov/) or consideration of pregrowing season conditions (Raynor et al. 2020) could improve this fine-tuning of stocking rates year to year.
Directly modeling forage provided more reliable estimates of spatial and temporal patterns of forage availability than treating ANPP as a proxy for forage availability. However, we know of few others who have directly modeled available forage rather than ANPP (though see Hirata et al. 2005 and Smolko et al. 2018). Doing so should improve forage availability and stocking rate determinations by addressing variations in productivity and vegetation palatability throughout an area. Our modeling method was simple and incorporated only freely available data in conjunction with field observations, so it could be easily replicated at low cost wherever training data is available. Such methods can help minimize the expenditure of time and money required for sampling by monitoring change remotely.
Effective, economical inventory and monitoring is a basic challenge in rangeland environments. Due to high spatiotemporal variability and large spatial scales, it remains difficult to adequately sample and summarize indicators essential for informed management, such as forage availability. We used a random forest model to predict annual forage availability across the Uintah and Ouray Reservation and surroundings, trained on field-collected plant production data and freely available climatic and biophysical data. Producing results from 1984 to 2018, this model greatly enhances standard rangeland vegetation inventory in the area by explicitly addressing variability in forage availability. Understanding this variability and how it differs regionally may help improve rangeland management by informing how to adjust stocking rates in atypical years and avoid overgrazing or insufficient forage availability in drought years.
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.rama.2021.07.007.
We would like to thank everyone at the Bureau of Indian Affairs who helped make this work possible, including Patricia Wright, Antonio Pingree, Gary Dean, and Chris Secakuku. Thank you to Shane Green, Dean Stacy, and Jedd Bodily from the Natural Resources Conservation Service for their help planning this study, to Susan Durham for guidance with data analysis, and to all technicians who helped with data collection. Lastly, we would like to thank the reviewers whose input helped improve the manuscript.